# Diff of /src/code/irrat.lisp

revision 1.5 by ram, Wed Oct 24 16:42:48 1990 UTC revision 1.6 by ram, Thu Jan 3 13:16:27 1991 UTC
# Line 96  Line 96
96
97  (defparameter *intexp-maximum-exponent* 10000)  (defparameter *intexp-maximum-exponent* 10000)
98
99    ;;; This function precisely calculates base raised to an integral power.  It
100    ;;; separates the cases by the sign of power, for efficiency reasons, as powers
101    ;;; can be calculated more efficiently if power is a positive integer.  Values
102    ;;; of power are calculated as positive integers, and inverted if negative.
103    ;;;
104  (defun intexp (base power)  (defun intexp (base power)
105    (when (> (abs power) *intexp-maximum-exponent*)    (when (> (abs power) *intexp-maximum-exponent*)
106      (cerror "Continue with calculation."      (cerror "Continue with calculation."
# Line 113  Line 118
118             (setq base (* base base))             (setq base (* base base))
119             (setq power nextn)))))             (setq power nextn)))))
120
;;; This function calculates x raised to the nth power.  It separates
;;; the  cases by the type of n, for efficiency reasons, as powers can
;;; be calculated more efficiently if n is a positive integer,  Therefore,
;;; All integer values of n are calculated as positive integers, and
;;; inverted if negative.
121
122    ;;; EXPT  --  Public
123    ;;;
124    ;;;    If an integer power of a rational, use INTEXP above.  Otherwise, do
125    ;;; floating point stuff.  If both args are real, we try %POW right off,
126    ;;; assuming it will return 0 if the result may be complex.  If so, we call
127    ;;; COMPLEX-POW which directly computes the complex result.  We also separate
128    ;;; the complex-real and real-complex cases from the general complex case.
129    ;;;
130  (defun expt (base power)  (defun expt (base power)
131    "Returns BASE raised to the POWER."    "Returns BASE raised to the POWER."
132    (if (zerop power)    (if (zerop power)
133        ;; This is wrong if power isn't an integer.        (1+ (* base power))
134        (typecase (realpart base)        (labels ((real-expt (base power rtype)
135          (single-float (coerce 1 'single-float))                   (let* ((fbase (coerce base 'double-float))
136          (double-float (coerce 1 'double-float))                          (fpower (coerce power 'double-float))
137          (t 1))                          (res (coerce (%pow fbase fpower) rtype)))
138        (number-dispatch ((base number) (power number))                     (if (and (zerop res) (minusp fbase))
139          (((foreach fixnum bignum ratio (complex rational)) integer)                         (multiple-value-bind (re im)
140           (intexp base power))                                              (complex-pow fbase fpower)
141          (((foreach single-float double-float) integer)                           (%make-complex (coerce re rtype) (coerce im rtype)))
142           (coerce (%pow (coerce base 'double-float)                         res)))
143                         (coerce power 'double-float))                 (complex-pow (fbase fpower)
144                   '(dispatch-type base)))                   (let ((pow (%pow (- fbase) fpower))
145          (((foreach fixnum bignum ratio single-float)                         (fpower*pi (* fpower pi)))
146            (foreach ratio single-float))                     (values (* pow (%cos fpower*pi))
147           (coerce (%pow (coerce base 'double-float)                             (* pow (%sin fpower*pi))))))
148                         (coerce power 'double-float))          (declare (inline real-expt))
149                   'single-float))          (number-dispatch ((base number) (power number))
150          (((foreach fixnum bignum ratio single-float double-float) double-float)            (((foreach fixnum (or bignum ratio) (complex rational)) integer)
151           (%pow (coerce base 'double-float) (coerce power 'double-float)))             (intexp base power))
152          (((complex rational) ratio)            (((foreach single-float double-float) integer)
153           (* (expt (abs base) power)             (real-expt base power '(dispatch-type base)))
154              (cis (* power (phase base)))))            (((foreach fixnum (or bignum ratio) single-float)
155          (((complex float) (foreach integer ratio))              (foreach ratio single-float))
156           (* (expt (abs base) power)             (real-expt base power 'single-float))
157              (cis (* power (phase base)))))            (((foreach fixnum (or bignum ratio) single-float double-float)
158          (((foreach fixnum bignum ratio single-float double-float) complex)              double-float)
159           (if (minusp base)             (real-expt base power 'double-float))
160               (/ (exp (* power (log (- base)))))            ((double-float single-float)
161               (exp (* power (log base)))))             (real-expt base power 'double-float))
162          (((foreach (complex float) (complex rational)) complex)            (((foreach (complex rational) (complex float)) rational)
163           (exp (* power (log base)))))))             (* (expt (abs base) power)
164                  (cis (* power (phase base)))))
165              (((foreach fixnum (or bignum ratio) single-float double-float)
166                complex)
167               (if (minusp base)
168                   (/ (exp (* power (truly-the float (log (- base))))))
169                   (exp (* power (truly-the float (log base))))))
170              (((foreach (complex float) (complex rational)) complex)
171               (exp (* power (log base))))))))
172
173  (defun log (number &optional (base nil base-p))  (defun log (number &optional (base nil base-p))
174    "Return the logarithm of NUMBER in the base BASE, which defaults to e."    "Return the logarithm of NUMBER in the base BASE, which defaults to e."

Legend:
 Removed from v.1.5 changed lines Added in v.1.6