share/numeric/expense.lisp

   1 ;;; -*- Mode: Lisp; Package: Macsyma -*-                                 ;;;
   2 ;;;    (c) Copyright 1984 the Regents of the University of California.   ;;;
   3 ;;;        All Rights Reserved.                                          ;;;
   4 ;;;        This work was produced under the sponsorship of the           ;;;
   5 ;;;        U.S. Department of Energy.  The Government retains            ;;;
   6 ;;;        certain rights therein.                                       ;;;
   7 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   8
   9 (macsyma-module expens)
  10
  11 (defmvar $cost_reciprocal 4
  12          "The expense of computing a floating point reciprocal in terms of
  13   scalar floating point additions on the CRAY-1(Note that this can be redefined
  14   for vector mode on the CRAY-1, another computer, or put in terms of absolute
  15   machine cycles.  However, all COST_-type variables would need to be
  16   consistently redefined.  Note further that EXPENSE would probably need to
  17   be completely rethought for a multiprocessor or data-flow machine)."
  18          fixnum
  19          modified-commands '$expense)
  20
  21 (defmvar $cost_divide 5
  22          "The expense of computing a floating point divide in terms of
  23   scalar floating point additions on the CRAY-1(For further discussion do:
  24   DESCRIBE(COST_RECIPROCAL) )."
  25          fixnum
  26          modified-commands '$expense)
  27
  28 (defmvar $cost_sqrt 29.
  29          "The expense of computing a floating point square root in terms of
  30   scalar floating point additions on the CRAY-1(For further discussion do:
  31   DESCRIBE(COST_RECIPROCAL) )."
  32          fixnum
  33          modified-commands '$expense)
  34
  35 (defmvar $cost_exp 30.
  36          "The expense of computing a floating point exponentiation in terms
  37   of scalar floating point additions on the CRAY-1(For further discussion do:
  38   DESCRIBE(COST_RECIPROCAL) )."
  39          fixnum
  40          modified-commands '$expense)
  41
  42 (defmvar $cost_sin_cos_log 32.
  43          "The expense of computing a floating point SIN, COS, or LOG in
  44   terms of scalar floating point additions on the CRAY-1.  Note that this
  45   is a mean of sorts for the three operations(For further discussion do:
  46   DESCRIBE(COST_RECIPROCAL) )."
  47          fixnum
  48          modified-commands '$expense)
  49
  50 (defmvar $cost_float_power (+ $cost_exp $cost_sin_cos_log)
  51          "The expense of computing a floating point power in terms of scalar
  52   floating point additions on the CRAY-1(For further discussion do:
  53   DESCRIBE(COST_RECIPROCAL) )."
  54          fixnum
  55          modified-commands '($expense $gather_exponents))
  56
  57 (defmvar $cost_hyper_trig 35.
  58          "The expense of computing a floating point ARCSIN, ARCCOS, ARCTAN,
  59   SINH, COSH, TANH, or TAN in terms of scalar floating point additions on the
  60   CRAY-1.  Note that this is a mean of sorts for these 7 different operations.
  61   (For further discussion do: DESCRIBE(COST_RECIPROCAL) )."
  62          fixnum
  63          modified-commands '$expense)
  64
  65 (defmvar $merge_ops '((mlist simp) $cvmgp $cvmgt)
  66          "A MACSYMA list of currently known CRAY-1 vector merge operations."
  67          modified-commands '($block_optimize $expense))
  68
  69 (defun multiplies-in-nth-power (nth)
  70    (cond ((< nth 2) 0)
  71          (t
  72           (let ((slow (bignump nth)))
  73             (do ((exin nth (cond (slow (- exin (* pw2 rem)))
  74                                  (t (- exin (* pw2 rem)))))
  75                  (rem 0)
  76                  (in-cut -2 (+ 1 in-cut rem))
  77                  (pw2 1 (cond (slow (+ pw2 pw2))
  78                               (t (+ pw2 pw2)))))
  79                 ((or (zerop exin) (> in-cut $cost_float_power))
  80                  (cond ((< in-cut $cost_float_power) in-cut)
  81                        (t $cost_float_power)))
  82               (declare (fixnum exin rem in-cut pw2))
  83               (setq rem (cond (slow (rem (quotient exin pw2) 2))
  84                               (t (/ (truncate exin pw2) 2)))))))))
  85
  86 ;;; the following macro is courtesy of gjc.
  87
  88 (defmacro eval&reduce (oper eval list
  89                             &aux (temp (gensym))
  90                             (val (gensym)))
  91           `(let ((,temp ,list))
  92              (do ((,val (funcall ,eval (pop ,temp))
  93                         (funcall ,oper ,val (funcall ,eval (pop ,temp)))))
  94                  ((null ,temp) ,val))))
  95
  96 (defun $expense (x)
  97   (cond (($mapatom x) 0)
  98         (t (let ((opr (caar x)))
  99              (cond ((member opr '(mplus mtimes) :test #'eq)
 100                     (let ((cl (+ (- (length x) 2)
 101                                  (eval&reduce '+ '$expense (cdr x)))))
 102                       (cond ((and (eq opr 'mtimes) (equal -1 (cadr x))) (1- cl))
 103                             (t cl))))
 104                    ((eq opr 'mexpt)
 105                     (let ((expon (caddr x)))
 106                       (+ ($expense (cadr x))
 107                          (cond ((integerp expon)
 108                                 (cond ((< expon 0)
 109                                        (+ $cost_reciprocal
 110                                           (multiplies-in-nth-power (- expon))))
 111                                       (t (multiplies-in-nth-power expon))))
 112                                (t (+ ($expense expon)
 113                                      $cost_exp
 114                                      (cond ((eq (cadr x) '$%e) 0)
 115                                            (t $cost_sin_cos_log))))))))
 116                    ((eq opr 'mminus) ($expense (cadr x)))
 117                    ((eq opr '%sqrt) (+ $cost_sqrt ($expense (cadr x))))
 118                    ((member opr $merge_ops :test #'eq) (+ 4
 119                                              ($expense (cadr x))
 120                                              ($expense (caddr x))
 121                                              ($expense (cadddr x))))
 122                    ((eq opr 'mquotient)
 123                     (cond ((member (cadr x) '(1 -1))
 124                            (+ $cost_reciprocal ($expense (caddr x))))
 125                           (t (+ (* $cost_divide (- (length x) 2))
 126                                 (eval&reduce '+ '$expense (cdr x))))))
 127                    ((member opr '(%acos %asin %atan %cosh %sinh %tan %tanh) :test #'eq)
 128                     (+ $cost_hyper_trig ($expense (cadr x))))
 129                    ((member opr '(%cos %log %sin) :test #'eq)
 130                     (+ $cost_sin_cos_log ($expense (cadr x))))
 131                    ((eq opr '$atan2)
 132                     (+ $cost_hyper_trig ($expense (cadr x)) ($expense (caddr x))))
 133                    (t
 134                     (mformat nil "Beware: ~M is not in function base of EXPENSE~%" opr)
 135                     (eval&reduce '+ '$expense (cdr x))))))))