/[cmucl]/src/code/irrat.lisp
ViewVC logotype

Diff of /src/code/irrat.lisp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.9.1.1 by wlott, Fri Jan 24 08:47:04 1992 UTC revision 1.64 by rtoy, Mon Feb 28 17:14:45 2011 UTC
# Line 3  Line 3 
3  ;;; **********************************************************************  ;;; **********************************************************************
4  ;;; This code was written as part of the CMU Common Lisp project at  ;;; This code was written as part of the CMU Common Lisp project at
5  ;;; Carnegie Mellon University, and has been placed in the public domain.  ;;; Carnegie Mellon University, and has been placed in the public domain.
 ;;; If you want to use this code or any part of CMU Common Lisp, please contact  
 ;;; Scott Fahlman or slisp-group@cs.cmu.edu.  
6  ;;;  ;;;
7  (ext:file-comment  (ext:file-comment
8    "$Header$")    "$Header$")
# Line 18  Line 16 
16  ;;;  ;;;
17    
18  (in-package "KERNEL")  (in-package "KERNEL")
19    (intl:textdomain "cmucl")
20    
21    
22  ;;;; Random constants, utility functions, and macros.  ;;;; Random constants, utility functions, and macros.
# Line 33  Line 32 
32                                         "%"                                         "%"
33                                         (string-upcase name)))))                                         (string-upcase name)))))
34      `(progn      `(progn
35         (proclaim '(inline ,function))         (declaim (inline ,function))
36         (export ',function)         (export ',function)
37         (alien:def-alien-routine (,name ,function) double-float         (alien:def-alien-routine (,name ,function) double-float
38           ,@(let ((results nil))           ,@(let ((results nil))
# Line 48  Line 47 
47    `((((foreach fixnum single-float bignum ratio))    `((((foreach fixnum single-float bignum ratio))
48       (coerce (,function (coerce ,var 'double-float)) 'single-float))       (coerce (,function (coerce ,var 'double-float)) 'single-float))
49      ((double-float)      ((double-float)
50       (,function ,var))))       (,function ,var))
51        #+double-double
52        ((double-double-float)
53         (,(symbolicate "DD-" function) ,var))))
54    
55  ); eval-when (compile load eval)  ); eval-when (compile load eval)
56    
# Line 58  Line 60 
60  ;;; Please refer to the Unix man pages for details about these routines.  ;;; Please refer to the Unix man pages for details about these routines.
61    
62  ;;; Trigonometric.  ;;; Trigonometric.
63  (def-math-rtn "sin" 1)  #-(and x86 (not sse2))
64  (def-math-rtn "cos" 1)  (progn
65  (def-math-rtn "tan" 1)    ;; For x86 (without sse2), we can use x87 instructions to implement
66      ;; these.  With sse2, we don't currently support that, so these
67      ;; should be disabled.
68      (def-math-rtn "sin" 1)
69      (def-math-rtn "cos" 1)
70      (def-math-rtn "tan" 1)
71      (def-math-rtn "atan" 1)
72      (def-math-rtn "atan2" 2))
73  (def-math-rtn "asin" 1)  (def-math-rtn "asin" 1)
74  (def-math-rtn "acos" 1)  (def-math-rtn "acos" 1)
 (def-math-rtn "atan" 1)  
 (def-math-rtn "atan2" 2)  
75  (def-math-rtn "sinh" 1)  (def-math-rtn "sinh" 1)
76  (def-math-rtn "cosh" 1)  (def-math-rtn "cosh" 1)
77  (def-math-rtn "tanh" 1)  (def-math-rtn "tanh" 1)
# Line 73  Line 80 
80  (def-math-rtn "atanh" 1)  (def-math-rtn "atanh" 1)
81    
82  ;;; Exponential and Logarithmic.  ;;; Exponential and Logarithmic.
83  (def-math-rtn "exp" 1)  #-(and x86 (not sse2))
84  (def-math-rtn "expm1" 1)  (progn
85  (def-math-rtn "log" 1)    (def-math-rtn "exp" 1)
86  (def-math-rtn "log10" 1)    (def-math-rtn "log" 1)
87  (def-math-rtn "log1p" 1)    (def-math-rtn "log10" 1))
88    
89  (def-math-rtn "pow" 2)  (def-math-rtn "pow" 2)
90  (def-math-rtn "cbrt" 1)  #-(or x86 sparc-v7 sparc-v8 sparc-v9)
91  (def-math-rtn "sqrt" 1)  (def-math-rtn "sqrt" 1)
92  (def-math-rtn "hypot" 2)  (def-math-rtn "hypot" 2)
93    
94    ;; Don't want log1p to use the x87 instruction.
95    #-(or hpux (and x86 (not sse2)))
96    (def-math-rtn "log1p" 1)
97    
98    ;; These are needed for use by byte-compiled files.  But don't use
99    ;; these with sse2 since we don't support using the x87 instructions
100    ;; here.
101    #+(and x86 (not sse2))
102    (progn
103      #+nil
104      (defun %sin (x)
105        (declare (double-float x)
106                 (values double-float))
107        (%sin x))
108      (defun %sin-quick (x)
109        (declare (double-float x)
110                 (values double-float))
111        (%sin-quick x))
112      #+nil
113      (defun %cos (x)
114        (declare (double-float x)
115                 (values double-float))
116        (%cos x))
117      (defun %cos-quick (x)
118        (declare (double-float x)
119                 (values double-float))
120        (%cos-quick x))
121      #+nil
122      (defun %tan (x)
123        (declare (double-float x)
124                 (values double-float))
125        (%tan x))
126      (defun %tan-quick (x)
127        (declare (double-float x)
128                 (values double-float))
129        (%tan-quick x))
130      (defun %atan (x)
131        (declare (double-float x)
132                 (values double-float))
133        (%atan x))
134      (defun %atan2 (x y)
135        (declare (double-float x y)
136                 (values double-float))
137        (%atan2 x y))
138      (defun %exp (x)
139        (declare (double-float x)
140                 (values double-float))
141        (%exp x))
142      (defun %log (x)
143        (declare (double-float x)
144                 (values double-float))
145        (%log x))
146      (defun %log10 (x)
147        (declare (double-float x)
148                 (values double-float))
149        (%log10 x))
150      #+nil ;; notyet
151      (defun %pow (x y)
152        (declare (type (double-float 0d0) x)
153                 (double-float y)
154                 (values (double-float 0d0)))
155        (%pow x y))
156      (defun %sqrt (x)
157        (declare (double-float x)
158                 (values double-float))
159        (%sqrt x))
160      (defun %scalbn (f ex)
161        (declare (double-float f)
162                 (type (signed-byte 32) ex)
163                 (values double-float))
164        (%scalbn f ex))
165      (defun %scalb (f ex)
166        (declare (double-float f ex)
167                 (values double-float))
168        (%scalb f ex))
169      (defun %logb (x)
170        (declare (double-float x)
171                 (values double-float))
172        (%logb x))
173      (defun %log1p (x)
174        (declare (double-float x)
175                 (values double-float))
176        (%log1p x))
177      ) ; progn
178    
179    
180    ;; As above for x86.  It also seems to be needed to handle
181    ;; constant-folding in the compiler.
182    #+(or sparc (and x86 sse2))
183    (progn
184      (defun %sqrt (x)
185        (declare (double-float x)
186                 (values double-float))
187        (%sqrt x))
188      )
189    
190    ;;; The standard libm routines for sin, cos, and tan on x86 (Linux)
191    ;;; and ppc are not very accurate for large arguments when compared to
192    ;;; sparc (and maxima).  This is basically caused by the fact that
193    ;;; those libraries do not do an accurate argument reduction.  The
194    ;;; following functions use some routines Sun's free fdlibm library to
195    ;;; do accurate reduction.  Then we call the standard C functions (or
196    ;;; vops for x86) on the reduced argument.  This produces much more
197    ;;; accurate values.
198    
199    #+(or ppc x86)
200    (progn
201    (declaim (inline %%ieee754-rem-pi/2))
202    ;; Basic argument reduction routine.  It returns two values: n and y
203    ;; such that (n + 8*k)*pi/2+y = x where |y|<pi/4 and n indicates in
204    ;; which octant the arg lies.  Y is actually computed in two parts,
205    ;; y[0] and y[1] such that the sum is y, for accuracy.
206    
207    (alien:def-alien-routine ("__ieee754_rem_pio2" %%ieee754-rem-pi/2) c-call:int
208      (x double-float)
209      (y (* double-float)))
210    
211    ;; Same as above, but instead of needing to pass an array in, the
212    ;; output array is broken up into two output values instead.  This is
213    ;; easier for the user, and we don't have to wrap calls with
214    ;; without-gcing.
215    (declaim (inline %ieee754-rem-pi/2))
216    (alien:def-alien-routine ("ieee754_rem_pio2" %ieee754-rem-pi/2) c-call:int
217      (x double-float)
218      (y0 double-float :out)
219      (y1 double-float :out))
220    
221    )
222    
223    #+(or ppc sse2)
224    (progn
225    (declaim (inline %%sin %%cos %%tan))
226    (macrolet ((frob (alien-name lisp-name)
227                 `(alien:def-alien-routine (,alien-name ,lisp-name) double-float
228                    (x double-float))))
229      (frob "sin" %%sin)
230      (frob "cos" %%cos)
231      (frob "tan" %%tan))
232    )
233    
234    #+(or ppc x86)
235    (macrolet
236        ((frob (sin cos tan)
237           `(progn
238              ;; In all of the routines below, we just compute the sum of
239              ;; y0 and y1 and use that as the (reduced) argument for the
240              ;; trig functions.  This is slightly less accurate than what
241              ;; fdlibm does, which calls special functions using y0 and
242              ;; y1 separately, for greater accuracy.  This isn't
243              ;; implemented, and some spot checks indicate that what we
244              ;; have here is accurate.
245              ;;
246              ;; For x86 with an fsin/fcos/fptan instruction, the pi/4 is
247              ;; probably too restrictive.
248              (defun %sin (x)
249                (declare (double-float x))
250                (if (< (abs x) (/ pi 4))
251                    (,sin x)
252                    ;; Argument reduction needed
253                    (multiple-value-bind (n y0 y1)
254                        (%ieee754-rem-pi/2 x)
255                      (let ((reduced (+ y0 y1)))
256                        (case (logand n 3)
257                          (0 (,sin reduced))
258                          (1 (,cos reduced))
259                          (2 (- (,sin reduced)))
260                          (3 (- (,cos reduced))))))))
261              (defun %cos (x)
262                (declare (double-float x))
263                (if (< (abs x) (/ pi 4))
264                    (,cos x)
265                    ;; Argument reduction needed
266                    (multiple-value-bind (n y0 y1)
267                        (%ieee754-rem-pi/2 x)
268                      (let ((reduced (+ y0 y1)))
269                        (case (logand n 3)
270                          (0 (,cos reduced))
271                          (1 (- (,sin reduced)))
272                          (2 (- (,cos reduced)))
273                          (3 (,sin reduced)))))))
274              (defun %tan (x)
275                (declare (double-float x))
276                (if (< (abs x) (/ pi 4))
277                    (,tan x)
278                    ;; Argument reduction needed
279                    (multiple-value-bind (n y0 y1)
280                        (%ieee754-rem-pi/2 x)
281                      (let ((reduced (+ y0 y1)))
282                        (if (evenp n)
283                            (,tan reduced)
284                            (- (/ (,tan reduced)))))))))))
285      ;; Don't want %sin-quick and friends with sse2.
286      #+(and x86 (not sse2))
287      (frob %sin-quick %cos-quick %tan-quick)
288      #+(or ppc sse2)
289      (frob %%sin %%cos %%tan))
290    
291    
292    
293  ;;;; Power functions.  ;;;; Power functions.
294    
# Line 98  Line 304 
304    
305  (defparameter *intexp-maximum-exponent* 10000)  (defparameter *intexp-maximum-exponent* 10000)
306    
307    (define-condition intexp-limit-error (error)
308      ((base :initarg :base :reader intexp-base)
309       (power :initarg :power :reader intexp-power))
310      (:report (lambda (condition stream)
311                 (format stream (intl:gettext "The absolute value of ~S exceeds limit ~S.")
312                         (intexp-power condition)
313                         *intexp-maximum-exponent*))))
314    
315  ;;; This function precisely calculates base raised to an integral power.  It  ;;; This function precisely calculates base raised to an integral power.  It
316  ;;; separates the cases by the sign of power, for efficiency reasons, as powers  ;;; separates the cases by the sign of power, for efficiency reasons, as powers
317  ;;; can be calculated more efficiently if power is a positive integer.  Values  ;;; can be calculated more efficiently if power is a positive integer.  Values
318  ;;; of power are calculated as positive integers, and inverted if negative.  ;;; of power are calculated as positive integers, and inverted if negative.
319  ;;;  ;;;
320  (defun intexp (base power)  (defun intexp (base power)
321      ;; Handle the special case of 1^power.  Maxima sometimes does this,
322      ;; and there's no need to cause a continuable error in this case.
323      ;; Should we also handle (-1)^power?
324      (when (eql base 1)
325        (return-from intexp base))
326    
327    (when (> (abs power) *intexp-maximum-exponent*)    (when (> (abs power) *intexp-maximum-exponent*)
328      (cerror "Continue with calculation."      ;; Allow user the option to continue with calculation, possibly
329              "The absolute value of ~S exceeds ~S."      ;; increasing the limit to the given power.
330              power '*intexp-maximum-exponent* base power))      (restart-case
331            (error 'intexp-limit-error
332                   :base base
333                   :power power)
334          (continue ()
335            :report (lambda (stream)
336                      (write-string (intl:gettext "Continue with calculation") stream)))
337          (new-limit ()
338            :report (lambda (stream)
339                      (write-string (intl:gettext "Continue with calculation, update limit") stream))
340            (setq *intexp-maximum-exponent* (abs power)))))
341    (cond ((minusp power)    (cond ((minusp power)
342           (/ (intexp base (- power))))           (/ (intexp base (- power))))
343          ((eql base 2)          ((eql base 2)
# Line 132  Line 362 
362  (defun expt (base power)  (defun expt (base power)
363    "Returns BASE raised to the POWER."    "Returns BASE raised to the POWER."
364    (if (zerop power)    (if (zerop power)
365        (1+ (* base power))        ;; CLHS says that if the power is 0, the result is 1, subject to
366        (labels ((real-expt (base power rtype)        ;; numeric contagion.  But what happens if base is infinity or
367                   (let* ((fbase (coerce base 'double-float))        ;; NaN?  Do we silently return 1?  For now, I think we should
368                          (fpower (coerce power 'double-float))        ;; signal an error if the FP modes say so.
369                          (res (coerce (%pow fbase fpower) rtype)))        (let ((result (1+ (* base power))))
370                     (if (and (zerop res) (minusp fbase))          ;; If we get an NaN here, that means base*power above didn't
371                         (multiple-value-bind (re im)          ;; produce 0 and FP traps were disabled, so we handle that
372                                              (complex-pow fbase fpower)          ;; here.  Should this be a continuable restart?
373                           (%make-complex (coerce re rtype) (coerce im rtype)))          (if (and (floatp result) (float-nan-p result))
374                         res)))              (float 1 result)
375                 (complex-pow (fbase fpower)              result))
376                   (let ((pow (%pow (- fbase) fpower))      (labels (;; determine if the double float is an integer.
377                         (fpower*pi (* fpower pi)))               ;;  0 - not an integer
378                     (values (* pow (%cos fpower*pi))               ;;  1 - an odd int
379                             (* pow (%sin fpower*pi))))))               ;;  2 - an even int
380          (declare (inline real-expt))               (isint (ihi lo)
381          (number-dispatch ((base number) (power number))                 (declare (type (unsigned-byte 31) ihi)
382            (((foreach fixnum (or bignum ratio) (complex rational)) integer)                          (type (unsigned-byte 32) lo)
383             (intexp base power))                          (optimize (speed 3) (safety 0)))
384            (((foreach single-float double-float) rational)                 (let ((isint 0))
385             (real-expt base power '(dispatch-type base)))                   (declare (type fixnum isint))
386            (((foreach fixnum (or bignum ratio) single-float)                   (cond ((>= ihi #x43400000)     ; exponent >= 53
387              (foreach ratio single-float))                          (setq isint 2))
388             (real-expt base power 'single-float))                         ((>= ihi #x3ff00000)
389            (((foreach fixnum (or bignum ratio) single-float double-float)                          (let ((k (- (ash ihi -20) #x3ff)))      ; exponent
390              double-float)                            (declare (type (mod 53) k))
391             (real-expt base power 'double-float))                            (cond ((> k 20)
392            ((double-float single-float)                                   (let* ((shift (- 52 k))
393             (real-expt base power 'double-float))                                          (j (logand (ash lo (- shift))))
394            (((foreach (complex rational) (complex float)) rational)                                          (j2 (ash j shift)))
395             (* (expt (abs base) power)                                     (declare (type (mod 32) shift)
396                (cis (* power (phase base)))))                                              (type (unsigned-byte 32) j j2))
397            (((foreach fixnum (or bignum ratio) single-float double-float)                                     (when (= j2 lo)
398              complex)                                       (setq isint (- 2 (logand j 1))))))
399             (if (minusp base)                                  ((= lo 0)
400                 (/ (exp (* power (truly-the float (log (- base))))))                                   (let* ((shift (- 20 k))
401                 (exp (* power (truly-the float (log base))))))                                          (j (ash ihi (- shift)))
402            (((foreach (complex float) (complex rational)) complex)                                          (j2 (ash j shift)))
403             (exp (* power (log base))))))))                                     (declare (type (mod 32) shift)
404                                                (type (unsigned-byte 31) j j2))
405                                       (when (= j2 ihi)
406                                         (setq isint (- 2 (logand j 1))))))))))
407                     isint))
408                 (real-expt (x y rtype)
409                   (let ((x (coerce x 'double-float))
410                         (y (coerce y 'double-float)))
411                     (declare (double-float x y))
412                     (let* ((x-hi (kernel:double-float-high-bits x))
413                            (x-lo (kernel:double-float-low-bits x))
414                            (x-ihi (logand x-hi #x7fffffff))
415                            (y-hi (kernel:double-float-high-bits y))
416                            (y-lo (kernel:double-float-low-bits y))
417                            (y-ihi (logand y-hi #x7fffffff)))
418                       (declare (type (signed-byte 32) x-hi y-hi)
419                                (type (unsigned-byte 31) x-ihi y-ihi)
420                                (type (unsigned-byte 32) x-lo y-lo))
421                       ;; y==zero: x**0 = 1
422                       (when (zerop (logior y-ihi y-lo))
423                         (return-from real-expt (coerce 1d0 rtype)))
424                       ;; +-NaN return x+y
425                       (when (or (> x-ihi #x7ff00000)
426                                 (and (= x-ihi #x7ff00000) (/= x-lo 0))
427                                 (> y-ihi #x7ff00000)
428                                 (and (= y-ihi #x7ff00000) (/= y-lo 0)))
429                         (return-from real-expt (coerce (+ x y) rtype)))
430                       (let ((yisint (if (< x-hi 0) (isint y-ihi y-lo) 0)))
431                         (declare (type fixnum yisint))
432                         ;; special value of y
433                         (when (and (zerop y-lo) (= y-ihi #x7ff00000))
434                           ;; y is +-inf
435                           (return-from real-expt
436                             (cond ((and (= x-ihi #x3ff00000) (zerop x-lo))
437                                    ;; +-1**inf is NaN
438                                    (coerce (- y y) rtype))
439                                   ((>= x-ihi #x3ff00000)
440                                    ;; (|x|>1)**+-inf = inf,0
441                                    (if (>= y-hi 0)
442                                        (coerce y rtype)
443                                        (coerce 0 rtype)))
444                                   (t
445                                    ;; (|x|<1)**-,+inf = inf,0
446                                    (if (< y-hi 0)
447                                        (coerce (- y) rtype)
448                                        (coerce 0 rtype))))))
449    
450                         (let ((abs-x (abs x)))
451                           (declare (double-float abs-x))
452                           ;; special value of x
453                           (when (and (zerop x-lo)
454                                      (or (= x-ihi #x7ff00000) (zerop x-ihi)
455                                          (= x-ihi #x3ff00000)))
456                             ;; x is +-0,+-inf,+-1
457                             (let ((z (if (< y-hi 0)
458                                          (/ 1 abs-x)       ; z = (1/|x|)
459                                          abs-x)))
460                               (declare (double-float z))
461                               (when (< x-hi 0)
462                                 (cond ((and (= x-ihi #x3ff00000) (zerop yisint))
463                                        ;; (-1)**non-int
464                                        (let ((y*pi (* y pi)))
465                                          (declare (double-float y*pi))
466                                          (return-from real-expt
467                                            (complex
468                                             (coerce (%cos y*pi) rtype)
469                                             (coerce (%sin y*pi) rtype)))))
470                                       ((= yisint 1)
471                                        ;; (x<0)**odd = -(|x|**odd)
472                                        (setq z (- z)))))
473                               (return-from real-expt (coerce z rtype))))
474    
475                           (if (>= x-hi 0)
476                               ;; x>0
477                               (coerce (kernel::%pow x y) rtype)
478                               ;; x<0
479                               (let ((pow (kernel::%pow abs-x y)))
480                                 (declare (double-float pow))
481                                 (case yisint
482                                   (1 ; Odd
483                                    (coerce (* -1d0 pow) rtype))
484                                   (2 ; Even
485                                    (coerce pow rtype))
486                                   (t ; Non-integer
487                                    (let ((y*pi (* y pi)))
488                                      (declare (double-float y*pi))
489                                      (complex
490                                       (coerce (* pow (%cos y*pi)) rtype)
491                                       (coerce (* pow (%sin y*pi)) rtype)))))))))))))
492          ;; This is really messy and should be cleaned up.  The easiest
493          ;; way to see if we're doing what we should is the macroexpand
494          ;; the number-dispatch and check each branch.
495          ;;
496          ;; We try to apply the rule of float precision contagion (CLHS
497          ;; 12.1.4.4): the result has the same precision has the most
498          ;; precise argument.
499          (number-dispatch ((base number) (power number))
500            (((foreach fixnum (or bignum ratio) (complex rational))
501              integer)
502             (intexp base power))
503            (((foreach single-float double-float)
504              rational)
505             (real-expt base power '(dispatch-type base)))
506            (((foreach fixnum (or bignum ratio) single-float)
507              (foreach ratio single-float))
508             (real-expt base power 'single-float))
509            (((foreach fixnum (or bignum ratio) single-float double-float)
510              double-float)
511             (real-expt base power 'double-float))
512            ((double-float single-float)
513             (real-expt base power 'double-float))
514            #+double-double
515            (((foreach fixnum (or bignum ratio) single-float double-float
516                       double-double-float)
517              double-double-float)
518             (dd-%pow (coerce base 'double-double-float) power))
519            #+double-double
520            ((double-double-float
521              (foreach fixnum (or bignum ratio) single-float double-float))
522             (dd-%pow base (coerce power 'double-double-float)))
523            (((foreach (complex rational) (complex single-float) (complex double-float)
524                       #+double-double (complex double-double-float))
525              rational)
526             (* (expt (abs base) power)
527                (cis (* power (phase base)))))
528            #+double-double
529            ((double-double-float
530              complex)
531             (if (and (zerop base) (plusp (realpart power)))
532                 (* base power)
533                 (exp (* power (* (log2 base 1w0) (log 2w0))))))
534            (((foreach fixnum (or bignum ratio) single-float double-float)
535              (foreach (complex double-float)))
536             ;; Result should have double-float accuracy.  Use log2 in
537             ;; case the base won't fit in a double-float.
538             (if (and (zerop base) (plusp (realpart power)))
539                 (* base power)
540                 (exp (* power (* (log2 base) (log 2d0))))))
541            ((double-float
542              (foreach (complex rational) (complex single-float)))
543             (if (and (zerop base) (plusp (realpart power)))
544                 (* base power)
545                 (exp (* power (log base)))))
546            #+double-double
547            (((foreach fixnum (or bignum ratio) single-float double-float)
548              (foreach (complex double-double-float)))
549             ;; Result should have double-double-float accuracy.  Use log2
550             ;; in case the base won't fit in a double-float.
551             (if (and (zerop base) (plusp (realpart power)))
552                 (* base power)
553                 (exp (* power (* (log2 base 1w0) (log 2w0))))))
554            (((foreach fixnum (or bignum ratio) single-float)
555              (foreach (complex single-float)))
556             (if (and (zerop base) (plusp (realpart power)))
557                 (* base power)
558                 (exp (* power (log base)))))
559            (((foreach (complex rational) (complex single-float))
560              (foreach single-float (complex single-float)))
561             (if (and (zerop base) (plusp (realpart power)))
562                 (* base power)
563                 (exp (* power (log base)))))
564            (((foreach (complex rational) (complex single-float))
565              (foreach double-float (complex double-float)))
566             (if (and (zerop base) (plusp (realpart power)))
567                 (* base power)
568                 (exp (* power (log (coerce base '(complex double-float)))))))
569            #+double-double
570            (((foreach (complex rational) (complex single-float))
571              (foreach double-double-float (complex double-double-float)))
572             (if (and (zerop base) (plusp (realpart power)))
573                 (* base power)
574                 (exp (* power (log (coerce base '(complex double-double-float)))))))
575            (((foreach (complex double-float))
576              (foreach single-float double-float (complex single-float)
577                       (complex double-float)))
578             (if (and (zerop base) (plusp (realpart power)))
579                 (* base power)
580                 (exp (* power (log base)))))
581            #+double-double
582            (((foreach (complex double-float))
583              (foreach double-double-float (complex double-double-float)))
584             (if (and (zerop base) (plusp (realpart power)))
585                 (* base power)
586                 (exp (* power (log (coerce base '(complex double-double-float)))))))
587            #+double-double
588            (((foreach (complex double-double-float))
589              (foreach float (complex float)))
590             (if (and (zerop base) (plusp (realpart power)))
591                 (* base power)
592                 (exp (* power (log base)))))))))
593    
594    ;; Log base 2 of a real number.  The result is a either a double-float
595    ;; or double-double-float number (real or complex, as appropriate),
596    ;; depending on the type of FLOAT-TYPE.
597    (defun log2 (x &optional (float-type 1d0))
598      (labels ((log-of-2 (f)
599                 ;; log(2), with the precision specified by the type of F
600                 (number-dispatch ((f real))
601                   ((double-float)
602                    #.(log 2d0))
603                   #+double-double
604                   ((double-double-float)
605                    #.(log 2w0))))
606               (log-2-pi (f)
607                 ;; log(pi), with the precision specified by the type of F
608                 (number-dispatch ((f real))
609                   ((double-float)
610                    #.(/ pi (log 2d0)))
611                   #+double-double
612                   ((double-double-float)
613                    #.(/ dd-pi (log 2w0)))))
614               (log1p (x)
615                 ;; log(1+x), with the precision specified by the type of
616                 ;; X
617                 (number-dispatch ((x real))
618                   (((foreach single-float double-float))
619                    (%log1p (float x 1d0)))
620                   #+double-double
621                   ((double-double-float)
622                    (dd-%log1p x))))
623               (log2-bignum (bignum)
624                 ;; Write x = 2^n*f where 1/2 < f <= 1.  Then log2(x) = n
625                 ;; + log2(f).
626                 ;;
627                 ;; So we grab the top few bits of x and scale that
628                 ;; appropriately, take the log of it and add it to n.
629                 ;;
630                 ;; Return n and log2(f) separately.
631                 (if (minusp bignum)
632                     (multiple-value-bind (n frac)
633                         (log2-bignum (abs bignum))
634                       (values n (complex frac (log-2-pi float-type))))
635                     (let ((n (integer-length bignum))
636                           (float-bits (float-digits float-type)))
637                       (if (< n float-bits)
638                           (values 0 (log (float bignum float-type)
639                                          (float 2 float-type)))
640                           (let ((exp (min float-bits n))
641                                 (f (ldb (byte float-bits
642                                               (max 0 (- n float-bits)))
643                                         bignum)))
644                             (values n (log (scale-float (float f float-type) (- exp))
645                                            (float 2 float-type)))))))))
646        (etypecase x
647          (float
648           (/ (log (float x float-type)) (log-of-2 float-type)))
649          (ratio
650           (let ((top (numerator x))
651                 (bot (denominator x)))
652             ;; If the number of bits in the numerator and
653             ;; denominator are different, just use the fact
654             ;; log(x/y) = log(x) - log(y).  But to preserve
655             ;; accuracy, we actually do
656             ;; (log2(x)-log2(y))/log2(e)).
657             ;;
658             ;; However, if the numerator and denominator have the
659             ;; same number of bits, implying the quotient is near
660             ;; one, we use log1p(x) = log(1+x). Since the number is
661             ;; rational, we don't lose precision subtracting 1 from
662             ;; it, and converting it to double-float is accurate.
663             (if (= (integer-length top)
664                    (integer-length bot))
665                 (/ (log1p (float (- x 1) float-type))
666                    (log-of-2 float-type))
667                 (multiple-value-bind (top-n top-frac)
668                     (log2-bignum top)
669                   (multiple-value-bind (bot-n bot-frac)
670                       (log2-bignum bot)
671                     (+ (- top-n bot-n)
672                        (- top-frac bot-frac)))))))
673          (integer
674           (multiple-value-bind (n frac)
675               (log2-bignum x)
676             (+ n frac))))))
677    
678  (defun log (number &optional (base nil base-p))  (defun log (number &optional (base nil base-p))
679    "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."
680    (if base-p    (if base-p
681        (/ (log number) (log base))        (cond ((zerop base)
682                 ;; ANSI spec
683                 base)
684                ((and (realp number) (realp base))
685                 ;; CLHS 12.1.4.1 says
686                 ;;
687                 ;;   When rationals and floats are combined by a
688                 ;;   numerical function, the rational is first converted
689                 ;;   to a float of the same format.
690                 ;;
691                 ;; So assume this applies to floats as well convert all
692                 ;; numbers to the largest float format before computing
693                 ;; the log.
694                 ;;
695                 ;; This makes (log 17 10.0) = (log 17.0 10) and so on.
696                 (number-dispatch ((number real) (base real))
697                   ((double-float
698                     (foreach double-float single-float))
699                    (/ (log2 number) (log2 base)))
700                   (((foreach fixnum bignum ratio)
701                     (foreach fixnum bignum ratio single-float))
702                    (let* ((result (/ (log2 number) (log2 base))))
703                      ;; Figure out the right result type
704                      (if (realp result)
705                          (coerce result 'single-float)
706                          (coerce result '(complex single-float)))))
707                   (((foreach fixnum bignum ratio)
708                     double-float)
709                    (/ (log2 number) (log2 base)))
710                   ((single-float
711                     (foreach fixnum bignum ratio))
712                    (let* ((result (/ (log2 number) (log2 base))))
713                      ;; Figure out the right result type
714                      (if (realp result)
715                          (coerce result 'single-float)
716                          (coerce result '(complex single-float)))))
717                   ((double-float
718                     (foreach fixnum bignum ratio))
719                    (/ (log2 number) (log2 base)))
720                   ((single-float double-float)
721                    (/ (log (coerce number 'double-float)) (log base)))
722                   #+double-double
723                   ((double-double-float
724                     (foreach fixnum bignum ratio))
725                    (/ (log2 number 1w0) (log2 base 1w0)))
726                   #+double-double
727                   ((double-double-float
728                     (foreach double-double-float double-float single-float))
729                    (/ (log number) (log (coerce base 'double-double-float))))
730                   #+double-double
731                   (((foreach fixnum bignum ratio)
732                     double-double-float)
733                    (/ (log2 number 1w0) (log2 base 1w0)))
734                   #+double-double
735                   (((foreach double-float single-float)
736                     double-double-float)
737                    (/ (log (coerce number 'double-double-float)) (log base)))
738                   (((foreach single-float)
739                     (foreach single-float))
740                    ;; Converting everything to double-float helps the
741                    ;; cases like (log 17 10) = (/ (log 17) (log 10)).
742                    ;; This is usually handled above, but if we compute (/
743                    ;; (log 17) (log 10)), we get a slightly different
744                    ;; answer due to roundoff.  This makes it a bit more
745                    ;; consistent.
746                    ;;
747                    ;; FIXME: This probably needs more work.
748                    (let ((result (/ (log (float number 1d0))
749                                     (log (float base 1d0)))))
750                      (if (realp result)
751                          (coerce result 'single-float)
752                          (coerce result '(complex single-float)))))))
753                (t
754                 ;; FIXME:  This probably needs some work as well.
755                 (/ (log number) (log base))))
756        (number-dispatch ((number number))        (number-dispatch ((number number))
757          (((foreach fixnum bignum ratio single-float))          (((foreach fixnum bignum))
758           (if (minusp number)           (if (minusp number)
759               (complex (log (- number)) (coerce pi 'single-float))               (complex (coerce (log (- number)) 'single-float)
760               (coerce (%log (coerce number 'double-float)) 'single-float)))                        (coerce pi 'single-float))
761          ((double-float)               (coerce (/ (log2 number) #.(log (exp 1d0) 2d0)) 'single-float)))
762            ((ratio)
763           (if (minusp number)           (if (minusp number)
764               (complex (log (- number)) (coerce pi 'double-float))               (complex (coerce (log (- number)) 'single-float)
765               (%log number)))                        (coerce pi 'single-float))
766          ((complex) (complex (log (abs number)) (phase number))))))               ;; What happens when the ratio is close to 1?  We need to
767                 ;; be careful to preserve accuracy.
768                 (let ((top (numerator number))
769                       (bot (denominator number)))
770                   ;; If the number of bits in the numerator and
771                   ;; denominator are different, just use the fact
772                   ;; log(x/y) = log(x) - log(y).  But to preserve
773                   ;; accuracy, we actually do
774                   ;; (log2(x)-log2(y))/log2(e)).
775                   ;;
776                   ;; However, if the numerator and denominator have the
777                   ;; same number of bits, implying the quotient is near
778                   ;; one, we use log1p(x) = log(1+x). Since the number is
779                   ;; rational, we don't lose precision subtracting 1 from
780                   ;; it, and converting it to double-float is accurate.
781                   (if (= (integer-length top)
782                          (integer-length bot))
783                       (coerce (%log1p (coerce (- number 1) 'double-float))
784                               'single-float)
785                       (coerce (/ (- (log2 top) (log2 bot))
786                                  #.(log (exp 1d0) 2d0))
787                               'single-float)))))
788            (((foreach single-float double-float))
789             ;; Is (log -0) -infinity (libm.a) or -infinity + i*pi (Kahan)?
790             ;; Since this doesn't seem to be an implementation issue
791             ;; I (pw) take the Kahan result.
792             (if (< (float-sign number)
793                    (coerce 0 '(dispatch-type number)))
794                 (complex (log (- number)) (coerce pi '(dispatch-type number)))
795                 (coerce (%log (coerce number 'double-float))
796                         '(dispatch-type number))))
797            #+double-double
798            ((double-double-float)
799             (let ((hi (kernel:double-double-hi number)))
800               (if (< (float-sign hi) 0d0)
801                   (complex (dd-%log (- number)) dd-pi)
802                   (dd-%log number))))
803            ((complex)
804             (complex-log number)))))
805    
806  (defun sqrt (number)  (defun sqrt (number)
807    "Return the square root of NUMBER."    "Return the square root of NUMBER."
808    (number-dispatch ((number number))    (number-dispatch ((number number))
809      (((foreach fixnum bignum ratio single-float))      (((foreach fixnum bignum ratio))
810       (if (minusp number)       (if (minusp number)
811           (exp (/ (log number) 2))           (complex-sqrt number)
812           (coerce (%sqrt (coerce number 'double-float)) 'single-float)))           (coerce (%sqrt (coerce number 'double-float)) 'single-float)))
813      ((double-float)      (((foreach single-float double-float))
814       (if (minusp number)       (if (minusp number)
815           (exp (/ (log number) 2))           (complex-sqrt number)
816           (%sqrt number)))           (coerce (%sqrt (coerce number 'double-float))
817      ((complex) (exp (/ (log number) 2)))))                   '(dispatch-type number))))
818        #+double-double
819  ;;; ISQRT:  Integer square root - isqrt(n)**2 <= n      ((double-double-float)
820  ;;; Upper and lower bounds on the result are estimated using integer-length.       (if (minusp number)
821  ;;; On each iteration, one of the bounds is replaced by their mean.           (dd-complex-sqrt number)
822  ;;; The lower bound is returned when the bounds meet or differ by only 1.           (multiple-value-bind (hi lo)
823  ;;; Initial bounds guarantee that lg(sqrt(n)) = lg(n)/2 iterations suffice.               (c::sqrt-dd (kernel:double-double-hi number) (kernel:double-double-lo number))
824               (kernel:%make-double-double-float hi lo))))
825  (defun isqrt (n)      ((complex)
826    "Returns the root of the nearest integer less than       (complex-sqrt number))))
    n which is a perfect square."  
   (if (and (integerp n) (not (minusp n)))  
       (do* ((lg (integer-length n))  
             (lo (ash 1 (ash (1- lg) -1)))  
             (hi (+ lo (ash lo (if (oddp lg) -1 0))))) ;tighten by 3/4 if possible.  
            ((<= (1- hi) lo) lo)  
         (let ((mid (ash (+ lo hi) -1)))  
           (if (<= (* mid mid) n) (setq lo mid) (setq hi mid))))  
       (error "Isqrt: ~S argument must be a nonnegative integer" n)))  
   
827    
828    
829  ;;;; Trigonometic and Related Functions  ;;;; Trigonometic and Related Functions
# Line 225  Line 831 
831  (defun abs (number)  (defun abs (number)
832    "Returns the absolute value of the number."    "Returns the absolute value of the number."
833    (number-dispatch ((number number))    (number-dispatch ((number number))
834      (((foreach single-float double-float fixnum rational))      (((foreach single-float double-float fixnum rational
835                   #+double-double double-double-float))
836       (abs number))       (abs number))
837      ((complex)      ((complex)
838       (let ((rx (realpart number))       (let ((rx (realpart number))
# Line 238  Line 845 
845                            (coerce ix 'double-float))                            (coerce ix 'double-float))
846                    'single-float))                    'single-float))
847           (double-float           (double-float
848            (%hypot rx ix)))))))            (%hypot rx ix))
849             #+double-double
850             (double-double-float
851              (multiple-value-bind (abs^2 scale)
852                  (dd-cssqs number)
853                (scale-float (sqrt abs^2) scale))))))))
854    
855  (defun phase (number)  (defun phase (number)
856    "Returns the angle part of the polar representation of a complex number.    "Returns the angle part of the polar representation of a complex number.
# Line 246  Line 858 
858    For non-complex positive numbers, this is 0.  For non-complex negative    For non-complex positive numbers, this is 0.  For non-complex negative
859    numbers this is PI."    numbers this is PI."
860    (etypecase number    (etypecase number
861      ((or rational single-float)      (rational
862       (if (minusp number)       (if (minusp number)
863           (coerce pi 'single-float)           (coerce pi 'single-float)
864           0.0f0))           0.0f0))
865        (single-float
866         (if (minusp (float-sign number))
867             (coerce pi 'single-float)
868             0.0f0))
869      (double-float      (double-float
870       (if (minusp number)       (if (minusp (float-sign number))
871           (coerce pi 'double-float)           (coerce pi 'double-float)
872           0.0d0))           0.0d0))
873        #+double-double
874        (double-double-float
875         (if (minusp (float-sign number))
876             dd-pi
877             0w0))
878      (complex      (complex
879       (atan (imagpart number) (realpart number)))))       (atan (imagpart number) (realpart number)))))
880    
# Line 265  Line 886 
886      ((complex)      ((complex)
887       (let ((x (realpart number))       (let ((x (realpart number))
888             (y (imagpart number)))             (y (imagpart number)))
889         (complex (* (sin x) (cosh y)) (* (cos x) (sinh y)))))))         (complex (* (sin x) (cosh y))
890                    (* (cos x) (sinh y)))))))
891    
892  (defun cos (number)  (defun cos (number)
893    "Return the cosine of NUMBER."    "Return the cosine of NUMBER."
# Line 274  Line 896 
896      ((complex)      ((complex)
897       (let ((x (realpart number))       (let ((x (realpart number))
898             (y (imagpart number)))             (y (imagpart number)))
899         (complex (* (cos x) (cosh y)) (- (* (sin x) (sinh y))))))))         (complex (* (cos x) (cosh y))
900                    (- (* (sin x) (sinh y))))))))
901    
902  (defun tan (number)  (defun tan (number)
903    "Return the tangent of NUMBER."    "Return the tangent of NUMBER."
904    (number-dispatch ((number number))    (number-dispatch ((number number))
905      (handle-reals %tan number)      (handle-reals %tan number)
906      ((complex)      ((complex)
907       (let* ((num (sin number))       (complex-tan number))))
             (denom (cos number)))  
        (if (zerop denom) (error "~S undefined tangent." number)  
            (/ num denom))))))  
908    
909  (defun cis (theta)  (defun cis (theta)
910    "Return cos(Theta) + i sin(Theta), AKA exp(i Theta)."    "Return cos(Theta) + i sin(Theta), AKA exp(i Theta)."
911    (if (complexp theta)    (if (complexp theta)
912        (error "Argument to CIS is complex: ~S" theta)        (error (intl:gettext "Argument to CIS is complex: ~S") theta)
913        (complex (cos theta) (sin theta))))        (complex (cos theta) (sin theta))))
914    
 (proclaim '(inline mult-by-i))  
 (defun mult-by-i (number)  
   (complex (imagpart number)  
            (- (realpart number))))  
   
 (defun complex-asin (number)  
   (- (mult-by-i (log (+ (mult-by-i number) (sqrt (- 1 (* number number))))))))  
   
915  (defun asin (number)  (defun asin (number)
916    "Return the arc sine of NUMBER."    "Return the arc sine of NUMBER."
917    (number-dispatch ((number number))    (number-dispatch ((number number))
# Line 308  Line 920 
920           (complex-asin number)           (complex-asin number)
921           (coerce (%asin (coerce number 'double-float)) 'single-float)))           (coerce (%asin (coerce number 'double-float)) 'single-float)))
922      (((foreach single-float double-float))      (((foreach single-float double-float))
923       (if (or (> number (coerce 1 '(dispatch-type number)))       (if (or (float-nan-p number)
924               (< number (coerce -1 '(dispatch-type number))))               (and (<= number (coerce 1 '(dispatch-type number)))
925           (complex-asin number)                    (>= number (coerce -1 '(dispatch-type number)))))
926           (coerce (%asin (coerce number 'double-float))           (coerce (%asin (coerce number 'double-float))
927                   '(dispatch-type number))))                   '(dispatch-type number))
928             (complex-asin number)))
929        #+double-double
930        ((double-double-float)
931         (if (or (float-nan-p number)
932                 (and (<= number 1w0)
933                      (>= number -1w0)))
934             (dd-%asin number)
935             (dd-complex-asin number)))
936      ((complex)      ((complex)
937       (complex-asin number))))       (complex-asin number))))
938    
 (defun complex-acos (number)  
   (- (mult-by-i (log (+ number (mult-by-i (sqrt (- (* number number)))))))))  
   
939  (defun acos (number)  (defun acos (number)
940    "Return the arc cosine of NUMBER."    "Return the arc cosine of NUMBER."
941    (number-dispatch ((number number))    (number-dispatch ((number number))
# Line 327  Line 944 
944           (complex-acos number)           (complex-acos number)
945           (coerce (%acos (coerce number 'double-float)) 'single-float)))           (coerce (%acos (coerce number 'double-float)) 'single-float)))
946      (((foreach single-float double-float))      (((foreach single-float double-float))
947       (if (or (> number (coerce 1 '(dispatch-type number)))       (if (or (float-nan-p number)
948               (< number (coerce -1 '(dispatch-type number))))               (and (<= number (coerce 1 '(dispatch-type number)))
949           (complex-acos number)                    (>= number (coerce -1 '(dispatch-type number)))))
950           (coerce (%acos (coerce number 'double-float))           (coerce (%acos (coerce number 'double-float))
951                   '(dispatch-type number))))                   '(dispatch-type number))
952             (complex-acos number)))
953        #+double-double
954        ((double-double-float)
955         (if (or (float-nan-p number)
956                 (and (<= number 1w0)
957                      (>= number -1w0)))
958             (dd-%acos number)
959             (complex-acos number)))
960      ((complex)      ((complex)
961       (complex-acos number))))       (complex-acos number))))
962    
# Line 339  Line 964 
964  (defun atan (y &optional (x nil xp))  (defun atan (y &optional (x nil xp))
965    "Return the arc tangent of Y if X is omitted or Y/X if X is supplied."    "Return the arc tangent of Y if X is omitted or Y/X if X is supplied."
966    (if xp    (if xp
967        (if (and (zerop x) (zerop y))        (flet ((atan2 (y x)
968            (if (plusp (float-sign (float x)))                 (declare (type double-float y x)
969                y                          (values double-float))
970                (if (minusp (float-sign (float y)))                 (if (zerop x)
971                    (- pi)                     (if (zerop y)
972                    pi))                         (if (plusp (float-sign x))
973            (number-dispatch ((y real) (x real))                             y
974              (((foreach fixnum bignum ratio single-float)                             (float-sign y pi))
975                (foreach fixnum bignum ratio single-float))                         (float-sign y (/ pi 2)))
976               (coerce (%atan2 (coerce y 'double-float)                     (%atan2 y x))))
977                               (coerce x 'double-float))          ;; If X is given, both X and Y must be real numbers.
978                       'single-float))          (number-dispatch ((y real) (x real))
979              ((double-float (foreach fixnum bignum ratio single-float))            ((double-float
980               (%atan2 y (coerce x 'double-float)))              (foreach double-float single-float fixnum bignum ratio))
981              (((foreach fixnum bignum ratio single-float double-float)             (atan2 y (coerce x 'double-float)))
982                double-float)            (((foreach single-float fixnum bignum ratio)
983               (%atan2 (coerce y 'double-float) x))))              double-float)
984               (atan2 (coerce y 'double-float) x))
985              (((foreach single-float fixnum bignum ratio)
986                (foreach single-float fixnum bignum ratio))
987               (coerce (atan2 (coerce y 'double-float) (coerce x 'double-float))
988                       'single-float))
989              #+double-double
990              ((double-double-float
991                (foreach double-double-float double-float single-float fixnum bignum ratio))
992               (dd-%atan2 y (coerce x 'double-double-float)))
993              #+double-double
994              (((foreach double-float single-float fixnum bignum ratio)
995                double-double-float)
996               (dd-%atan2 (coerce y 'double-double-float) x))))
997        (number-dispatch ((y number))        (number-dispatch ((y number))
998          (handle-reals %atan y)          (handle-reals %atan y)
999          ((complex)          ((complex)
1000           (let ((im (imagpart y))           (complex-atan y)))))
                (re (realpart y)))  
            (/ (- (log (complex (- 1 im) re))  
                  (log (complex (+ 1 im) (- re))))  
               (complex 0 2)))))))  
   
1001    
1002  (defun sinh (number)  (defun sinh (number)
1003    "Return the hyperbolic sine of NUMBER."    "Return the hyperbolic sine of NUMBER."
1004    (/ (- (exp number) (exp (- number))) 2))    (number-dispatch ((number number))
1005        (handle-reals %sinh number)
1006        ((complex)
1007         (let ((x (realpart number))
1008               (y (imagpart number)))
1009           (complex (* (sinh x) (cos y))
1010                    (* (cosh x) (sin y)))))))
1011    
1012  (defun cosh (number)  (defun cosh (number)
1013    "Return the hyperbolic cosine of NUMBER."    "Return the hyperbolic cosine of NUMBER."
1014    (/ (+ (exp number) (exp (- number))) 2))    (number-dispatch ((number number))
1015        (handle-reals %cosh number)
1016        ((complex)
1017         (let ((x (realpart number))
1018               (y (imagpart number)))
1019           (complex (* (cosh x) (cos y))
1020                    (* (sinh x) (sin y)))))))
1021    
1022  (defun tanh (number)  (defun tanh (number)
1023    "Return the hyperbolic tangent of NUMBER."    "Return the hyperbolic tangent of NUMBER."
1024    (/ (- (exp number) (exp (- number)))    (number-dispatch ((number number))
1025       (+ (exp number) (exp (- number)))))      (handle-reals %tanh number)
1026        ((complex)
1027         (complex-tanh number))))
1028    
1029  (defun asinh (number)  (defun asinh (number)
1030    "Return the hyperbolic arc sine of NUMBER."    "Return the hyperbolic arc sine of NUMBER."
1031    (log (+ number (sqrt (1+ (* number number))))))    (number-dispatch ((number number))
1032        (handle-reals %asinh number)
1033        ((complex)
1034         (complex-asinh number))))
1035    
1036  (defun acosh (number)  (defun acosh (number)
1037    "Return the hyperbolic arc cosine of NUMBER."    "Return the hyperbolic arc cosine of NUMBER."
1038    (log (+ number (* (1+ number) (sqrt (/ (1- number) (1+ number)))))))    (number-dispatch ((number number))
1039        ((rational)
1040         ;; acosh is complex if number < 1
1041         (if (< number 1)
1042             (complex-acosh number)
1043             (coerce (%acosh (coerce number 'double-float)) 'single-float)))
1044        (((foreach single-float double-float))
1045         (if (< number (coerce 1 '(dispatch-type number)))
1046             (complex-acosh number)
1047             (coerce (%acosh (coerce number 'double-float))
1048                     '(dispatch-type number))))
1049        #+double-double
1050        ((double-double-float)
1051         (if (< number 1w0)
1052             (complex-acosh number)
1053             (dd-%acosh number)))
1054        ((complex)
1055         (complex-acosh number))))
1056    
1057  (defun atanh (number)  (defun atanh (number)
1058    "Return the hyperbolic arc tangent of NUMBER."    "Return the hyperbolic arc tangent of NUMBER."
1059    (log (* (1+ number) (sqrt (/ (- 1 (* number number)))))))    (number-dispatch ((number number))
1060        ((rational)
1061         ;; atanh is complex if |number| > 1
1062         (if (or (> number 1) (< number -1))
1063             (complex-atanh number)
1064             (coerce (%atanh (coerce number 'double-float)) 'single-float)))
1065        (((foreach single-float double-float))
1066         (if (or (> number (coerce 1 '(dispatch-type number)))
1067                 (< number (coerce -1 '(dispatch-type number))))
1068             (complex-atanh number)
1069             (coerce (%atanh (coerce number 'double-float))
1070                     '(dispatch-type number))))
1071        #+double-double
1072        ((double-double-float)
1073         (if (or (> number 1w0)
1074                 (< number -1w0))
1075             (complex-atanh number)
1076             (dd-%atanh (coerce number 'double-double-float))))
1077        ((complex)
1078         (complex-atanh number))))
1079    
1080    ;;; HP-UX does not supply a C version of log1p, so use the definition.
1081    ;;; We really need to fix this.  The definition really loses big-time
1082    ;;; in roundoff as x gets small.
1083    
1084    #+hpux
1085    (declaim (inline %log1p))
1086    #+hpux
1087    (defun %log1p (number)
1088      (declare (double-float number)
1089               (optimize (speed 3) (safety 0)))
1090      (the double-float (log (the (double-float 0d0) (+ number 1d0)))))
1091    
1092    
1093    ;;;;
1094    ;;;; This is a set of routines that implement many elementary
1095    ;;;; transcendental functions as specified by ANSI Common Lisp.  The
1096    ;;;; implementation is based on Kahan's paper.
1097    ;;;;
1098    ;;;; I believe I have accurately implemented the routines and are
1099    ;;;; correct, but you may want to check for your self.
1100    ;;;;
1101    ;;;; These functions are written for CMU Lisp and take advantage of
1102    ;;;; some of the features available there.  It may be possible,
1103    ;;;; however, to port this to other Lisps.
1104    ;;;;
1105    ;;;; Some functions are significantly more accurate than the original
1106    ;;;; definitions in CMU Lisp.  In fact, some functions in CMU Lisp
1107    ;;;; give the wrong answer like (acos #c(-2.0 0.0)), where the true
1108    ;;;; answer is pi + i*log(2-sqrt(3)).
1109    ;;;;
1110    ;;;; All of the implemented functions will take any number for an
1111    ;;;; input, but the result will always be a either a complex
1112    ;;;; single-float or a complex double-float.
1113    ;;;;
1114    ;;;; General functions
1115    ;;;;   complex-sqrt
1116    ;;;;   complex-log
1117    ;;;;   complex-atanh
1118    ;;;;   complex-tanh
1119    ;;;;   complex-acos
1120    ;;;;   complex-acosh
1121    ;;;;   complex-asin
1122    ;;;;   complex-asinh
1123    ;;;;   complex-atan
1124    ;;;;   complex-tan
1125    ;;;;
1126    ;;;; Utility functions:
1127    ;;;;   scalb logb
1128    ;;;;
1129    ;;;; Internal functions:
1130    ;;;;    square coerce-to-complex-type cssqs complex-log-scaled
1131    ;;;;
1132    ;;;;
1133    ;;;; Please send any bug reports, comments, or improvements to Raymond
1134    ;;;; Toy at toy@rtp.ericsson.se.
1135    ;;;;
1136    ;;;; References
1137    ;;;;
1138    ;;;; Kahan, W. "Branch Cuts for Complex Elementary Functions, or Much
1139    ;;;; Ado About Nothing's Sign Bit" in Iserles and Powell (eds.) "The
1140    ;;;; State of the Art in Numerical Analysis", pp. 165-211, Clarendon
1141    ;;;; Press, 1987
1142    ;;;;
1143    
1144    (declaim (inline square))
1145    (defun square (x)
1146      (declare (float x))
1147      (* x x))
1148    
1149    ;; If you have these functions in libm, perhaps they should be used
1150    ;; instead of these Lisp versions.  These versions are probably good
1151    ;; enough, especially since they are portable.
1152    
1153    (declaim (inline scalb))
1154    (defun scalb (x n)
1155      "Compute 2^N * X without compute 2^N first (use properties of the
1156    underlying floating-point format"
1157      (declare (type float x)
1158               (type double-float-exponent n))
1159      (scale-float x n))
1160    
1161    (declaim (inline logb-finite))
1162    (defun logb-finite (x)
1163      "Same as logb but X is not infinity and non-zero and not a NaN, so
1164    that we can always return an integer"
1165      (declare (type float x))
1166      (multiple-value-bind (signif expon sign)
1167          (decode-float x)
1168        (declare (ignore signif sign))
1169        ;; decode-float is almost right, except that the exponent
1170        ;; is off by one
1171        (1- expon)))
1172    
1173    (defun logb (x)
1174      "Compute an integer N such that 1 <= |2^(-N) * x| < 2.
1175    For the special cases, the following values are used:
1176    
1177        x             logb
1178       NaN            NaN
1179       +/- infinity   +infinity
1180       0              -infinity
1181    "
1182      (declare (type float x))
1183      (cond ((float-nan-p x)
1184             x)
1185            ((float-infinity-p x)
1186             #.ext:double-float-positive-infinity)
1187            ((zerop x)
1188             ;; The answer is negative infinity, but we are supposed to
1189             ;; signal divide-by-zero, so do the actual division
1190             (/ -1 x)
1191             )
1192            (t
1193             (logb-finite x))))
1194    
1195    
1196    
1197    ;; This function is used to create a complex number of the appropriate
1198    ;; type.
1199    
1200    (declaim (inline coerce-to-complex-type))
1201    (defun coerce-to-complex-type (x y z)
1202      "Create complex number with real part X and imaginary part Y such that
1203    it has the same type as Z.  If Z has type (complex rational), the X
1204    and Y are coerced to single-float."
1205      (declare (double-float x y)
1206               (number z)
1207               (optimize (extensions:inhibit-warnings 3)))
1208      (if (typep (realpart z) 'double-float)
1209          (complex x y)
1210          ;; Convert anything that's not a double-float to a single-float.
1211          (complex (float x 1f0)
1212                   (float y 1f0))))
1213    
1214    (defun cssqs (z)
1215      ;; Compute |(x+i*y)/2^k|^2 scaled to avoid over/underflow. The
1216      ;; result is r + i*k, where k is an integer.
1217    
1218      ;; Save all FP flags
1219      (let ((x (float (realpart z) 1d0))
1220            (y (float (imagpart z) 1d0)))
1221        ;; Would this be better handled using an exception handler to
1222        ;; catch the overflow or underflow signal?  For now, we turn all
1223        ;; traps off and look at the accrued exceptions to see if any
1224        ;; signal would have been raised.
1225        (with-float-traps-masked (:underflow :overflow)
1226          (let ((rho (+ (square x) (square y))))
1227            (declare (optimize (speed 3) (space 0)))
1228            (cond ((and (or (float-nan-p rho)
1229                            (float-infinity-p rho))
1230                        (or (float-infinity-p (abs x))
1231                            (float-infinity-p (abs y))))
1232                   (values ext:double-float-positive-infinity 0))
1233                  ((let ((threshold #.(/ least-positive-double-float
1234                                         double-float-epsilon))
1235                         (traps (ldb vm::float-sticky-bits
1236                                     (vm:floating-point-modes))))
1237                     ;; Overflow raised or (underflow raised and rho <
1238                     ;; lambda/eps)
1239                     (or (not (zerop (logand vm:float-overflow-trap-bit traps)))
1240                         (and (not (zerop (logand vm:float-underflow-trap-bit traps)))
1241                              (< rho threshold))))
1242                   ;; If we're here, neither x nor y are infinity and at
1243                   ;; least one is non-zero.. Thus logb returns a nice
1244                   ;; integer.
1245                   (let ((k (- (logb-finite (max (abs x) (abs y))))))
1246                     (values (+ (square (scalb x k))
1247                                (square (scalb y k)))
1248                             (- k))))
1249                  (t
1250                   (values rho 0)))))))
1251    
1252    (defun complex-sqrt (z)
1253      "Principle square root of Z
1254    
1255    Z may be any number, but the result is always a complex."
1256      (declare (number z))
1257      #+double-double
1258      (when (typep z '(or double-double-float (complex double-double-float)))
1259        (return-from complex-sqrt (dd-complex-sqrt z)))
1260      (multiple-value-bind (rho k)
1261          (cssqs z)
1262        (declare (type (or (member 0d0) (double-float 0d0)) rho)
1263                 (type fixnum k))
1264        (let ((x (float (realpart z) 1.0d0))
1265              (y (float (imagpart z) 1.0d0))
1266              (eta 0d0)
1267              (nu 0d0))
1268          (declare (double-float x y eta nu))
1269    
1270          (locally
1271              ;; space 0 to get maybe-inline functions inlined.
1272              (declare (optimize (speed 3) (space 0)))
1273    
1274            (if (not (locally (declare (optimize (inhibit-warnings 3)))
1275                       (float-nan-p x)))
1276                (setf rho (+ (scalb (abs x) (- k)) (sqrt rho))))
1277    
1278            (cond ((oddp k)
1279                   (setf k (ash k -1)))
1280                  (t
1281                   (setf k (1- (ash k -1)))
1282                   (setf rho (+ rho rho))))
1283    
1284            (setf rho (scalb (sqrt rho) k))
1285    
1286            (setf eta rho)
1287            (setf nu y)
1288    
1289            (when (/= rho 0d0)
1290              (when (not (float-infinity-p (abs nu)))
1291                (setf nu (/ (/ nu rho) 2d0)))
1292              (when (< x 0d0)
1293                (setf eta (abs nu))
1294                (setf nu (float-sign y rho))))
1295            (coerce-to-complex-type eta nu z)))))
1296    
1297    (defun complex-log-scaled (z j)
1298      "Compute log(2^j*z).
1299    
1300    This is for use with J /= 0 only when |z| is huge."
1301      (declare (number z)
1302               (fixnum j))
1303      ;; The constants t0, t1, t2 should be evaluated to machine
1304      ;; precision.  In addition, Kahan says the accuracy of log1p
1305      ;; influences the choices of these constants but doesn't say how to
1306      ;; choose them.  We'll just assume his choices matches our
1307      ;; implementation of log1p.
1308      (let ((t0 #.(/ 1 (sqrt 2.0d0)))
1309            (t1 1.2d0)
1310            (t2 3d0)
1311            (ln2 #.(log 2d0))
1312            (x (float (realpart z) 1.0d0))
1313            (y (float (imagpart z) 1.0d0)))
1314        (multiple-value-bind (rho k)
1315            (cssqs z)
1316          (declare (optimize (speed 3)))
1317          (let ((beta (max (abs x) (abs y)))
1318                (theta (min (abs x) (abs y))))
1319            (coerce-to-complex-type (if (and (zerop k)
1320                                             (< t0 beta)
1321                                             (or (<= beta t1)
1322                                                 (< rho t2)))
1323                                        (/ (%log1p (+ (* (- beta 1.0d0)
1324                                                         (+ beta 1.0d0))
1325                                                      (* theta theta)))
1326                                           2d0)
1327                                        (+ (/ (log rho) 2d0)
1328                                           (* (+ k j) ln2)))
1329                                    (atan y x)
1330                                    z)))))
1331    
1332    (defun complex-log (z)
1333      "Log of Z = log |Z| + i * arg Z
1334    
1335    Z may be any number, but the result is always a complex."
1336      (declare (number z))
1337      #+double-double
1338      (when (typep z '(or double-double-float (complex double-double-float)))
1339        (return-from complex-log (dd-complex-log-scaled z 0)))
1340      (complex-log-scaled z 0))
1341    
1342    ;; Let us note the following "strange" behavior.  atanh 1.0d0 is
1343    ;; +infinity, but the following code returns approx 176 + i*pi/4. The
1344    ;; reason for the imaginary part is caused by the fact that arg i*y is
1345    ;; never 0 since we have positive and negative zeroes.
1346    
1347    (defun complex-atanh (z)
1348      "Compute atanh z = (log(1+z) - log(1-z))/2"
1349      (declare (number z))
1350      #+double-double
1351      (when (typep z '(or double-double-float (complex double-double-float)))
1352        (return-from complex-atanh (dd-complex-atanh z)))
1353    
1354      (if (and (realp z) (< z -1))
1355          ;; atanh is continuous in quadrant III in this case.
1356          (complex-atanh (complex z -0f0))
1357          (let* ( ;; Constants
1358                 (theta (/ (sqrt most-positive-double-float) 4.0d0))
1359                 (rho (/ 4.0d0 (sqrt most-positive-double-float)))
1360                 (half-pi (/ pi 2.0d0))
1361                 (rp (float (realpart z) 1.0d0))
1362                 (beta (float-sign rp 1.0d0))
1363                 (x (* beta rp))
1364                 (y (* beta (- (float (imagpart z) 1.0d0))))
1365                 (eta 0.0d0)
1366                 (nu 0.0d0))
1367            ;; Shouldn't need this declare.
1368            (declare (double-float x y))
1369            (locally
1370                (declare (optimize (speed 3)))
1371              (cond ((or (> x theta)
1372                         (> (abs y) theta))
1373                     ;; To avoid overflow...
1374                     (setf nu (float-sign y half-pi))
1375                     ;; eta is real part of 1/(x + iy).  This is x/(x^2+y^2),
1376                     ;; which can cause overflow.  Arrange this computation so
1377                     ;; that it won't overflow.
1378                     (setf eta (let* ((x-bigger (> x (abs y)))
1379                                      (r (if x-bigger (/ y x) (/ x y)))
1380                                      (d (+ 1.0d0 (* r r))))
1381                                 (if x-bigger
1382                                     (/ (/ x) d)
1383                                     (/ (/ r y) d)))))
1384                    ((= x 1.0d0)
1385                     ;; Should this be changed so that if y is zero, eta is set
1386                     ;; to +infinity instead of approx 176?  In any case
1387                     ;; tanh(176) is 1.0d0 within working precision.
1388                     (let ((t1 (+ 4d0 (square y)))
1389                           (t2 (+ (abs y) rho)))
1390                       (setf eta (log (/ (sqrt (sqrt t1))
1391                                         (sqrt t2))))
1392                       (setf nu (* 0.5d0
1393                                   (float-sign y
1394                                               (+ half-pi (atan (* 0.5d0 t2))))))))
1395                    (t
1396                     (let ((t1 (+ (abs y) rho)))
1397                       ;; Normal case using log1p(x) = log(1 + x)
1398                       (setf eta (* 0.25d0
1399                                    (%log1p (/ (* 4.0d0 x)
1400                                               (+ (square (- 1.0d0 x))
1401                                                  (square t1))))))
1402                       (setf nu (* 0.5d0
1403                                   (atan (* 2.0d0 y)
1404                                         (- (* (- 1.0d0 x)
1405                                               (+ 1.0d0 x))
1406                                            (square t1))))))))
1407              (coerce-to-complex-type (* beta eta)
1408                                      (- (* beta nu))
1409                                      z)))))
1410    
1411    (defun complex-tanh (z)
1412      "Compute tanh z = sinh z / cosh z"
1413      (declare (number z))
1414      #+double-double
1415      (when (typep z '(or double-double-float (complex double-double-float)))
1416        (return-from complex-tanh (dd-complex-tanh z)))
1417    
1418      (let ((x (float (realpart z) 1.0d0))
1419            (y (float (imagpart z) 1.0d0)))
1420        (locally
1421            ;; space 0 to get maybe-inline functions inlined
1422            (declare (optimize (speed 3) (space 0)))
1423          (cond ((> (abs x)
1424                    #-(or linux hpux) #.(/ (%asinh most-positive-double-float) 4d0)
1425                    ;; This is more accurate under linux.
1426                    #+(or linux hpux) #.(/ (+ (%log 2.0d0)
1427                                              (%log most-positive-double-float)) 4d0))
1428                 (coerce-to-complex-type (float-sign x)
1429                                         (float-sign y) z))
1430                (t
1431                 (let* ((tv (%tan y))
1432                        (beta (+ 1.0d0 (* tv tv)))
1433                        (s (sinh x))
1434                        (rho (sqrt (+ 1.0d0 (* s s)))))
1435                   (if (float-infinity-p (abs tv))
1436                       (coerce-to-complex-type (/ rho s)
1437                                               (/ tv)
1438                                               z)
1439                       (let ((den (+ 1.0d0 (* beta s s))))
1440                         (coerce-to-complex-type (/ (* beta rho s)
1441                                                    den)
1442                                                 (/ tv den)
1443                                                 z)))))))))
1444    
1445    ;; Kahan says we should only compute the parts needed.  Thus, the
1446    ;; realpart's below should only compute the real part, not the whole
1447    ;; complex expression.  Doing this can be important because we may get
1448    ;; spurious signals that occur in the part that we are not using.
1449    ;;
1450    ;; However, we take a pragmatic approach and just use the whole
1451    ;; expression.
1452    
1453    ;; NOTE: The formula given by Kahan is somewhat ambiguous in whether
1454    ;; it's the conjugate of the square root or the square root of the
1455    ;; conjugate.  This needs to be checked.
1456    
1457    ;; I checked.  It doesn't matter because (conjugate (sqrt z)) is the
1458    ;; same as (sqrt (conjugate z)) for all z.  This follows because
1459    ;;
1460    ;; (conjugate (sqrt z)) = exp(0.5*log |z|)*exp(-0.5*j*arg z).
1461    ;;
1462    ;; (sqrt (conjugate z)) = exp(0.5*log|z|)*exp(0.5*j*arg conj z)
1463    ;;
1464    ;; and these two expressions are equal if and only if arg conj z =
1465    ;; -arg z, which is clearly true for all z.
1466    
1467    ;; NOTE: The rules of Common Lisp says that if you mix a real with a
1468    ;; complex, the real is converted to a complex before performing the
1469    ;; operation.  However, Kahan says in this paper (pg 176):
1470    ;;
1471    ;; (iii) Careless handling can turn infinity or the sign of zero into
1472    ;;       misinformation that subsequently disappears leaving behind
1473    ;;       only a plausible but incorrect result.  That is why compilers
1474    ;;       must not transform z-1 into z-(1+i*0), as we have seen above,
1475    ;;       nor -(-x-x^2) into (x+x^2), as we shall see below, lest a
1476    ;;       subsequent logarithm or square root produce a non-zero
1477    ;;       imaginary part whose sign is opposite to what was intended.
1478    ;;
1479    ;; The interesting examples are too long and complicated to reproduce
1480    ;; here.  We refer the reader to his paper.
1481    ;;
1482    ;; The functions below are intended to handle the cases where a real
1483    ;; is mixed with a complex and we don't want CL complex contagion to
1484    ;; occur..
1485    
1486    (declaim (inline 1+z 1-z z-1 z+1))
1487    (defun 1+z (z)
1488      (complex (+ 1 (realpart z)) (imagpart z)))
1489    (defun 1-z (z)
1490      (complex (- 1 (realpart z)) (- (imagpart z))))
1491    (defun z-1 (z)
1492      (complex (- (realpart z) 1) (imagpart z)))
1493    (defun z+1 (z)
1494      (complex (+ (realpart z) 1) (imagpart z)))
1495    
1496    (defun complex-acos (z)
1497      "Compute acos z = pi/2 - asin z
1498    
1499    Z may be any number, but the result is always a complex."
1500      (declare (number z))
1501      #+double-double
1502      (when (typep z '(or double-double-float (complex double-double-float)))
1503        (return-from complex-acos (dd-complex-acos z)))
1504      (if (and (realp z) (> z 1))
1505          ;; acos is continuous in quadrant IV in this case.
1506          (complex-acos (complex z -0f0))
1507          (let ((sqrt-1+z (complex-sqrt (1+z z)))
1508                (sqrt-1-z (complex-sqrt (1-z z))))
1509            (with-float-traps-masked (:divide-by-zero)
1510              (complex (* 2 (atan (/ (realpart sqrt-1-z)
1511                                     (realpart sqrt-1+z))))
1512                       (asinh (imagpart (* (conjugate sqrt-1+z)
1513                                           sqrt-1-z))))))))
1514    
1515    (defun complex-acosh (z)
1516      "Compute acosh z = 2 * log(sqrt((z+1)/2) + sqrt((z-1)/2))
1517    
1518    Z may be any number, but the result is always a complex."
1519      (declare (number z))
1520      (let ((sqrt-z-1 (complex-sqrt (z-1 z)))
1521            (sqrt-z+1 (complex-sqrt (z+1 z))))
1522        (with-float-traps-masked (:divide-by-zero)
1523          (complex (asinh (realpart (* (conjugate sqrt-z-1)
1524                                       sqrt-z+1)))
1525                   (* 2 (atan (/ (imagpart sqrt-z-1)
1526                                 (realpart sqrt-z+1))))))))
1527    
1528    
1529    (defun complex-asin (z)
1530      "Compute asin z = asinh(i*z)/i
1531    
1532    Z may be any number, but the result is always a complex."
1533      (declare (number z))
1534      #+double-double
1535      (when (typep z '(or double-double-float (complex double-double-float)))
1536        (return-from complex-asin (dd-complex-asin z)))
1537      (if (and (realp z) (> z 1))
1538          ;; asin is continuous in quadrant IV in this case.
1539          (complex-asin (complex z -0f0))
1540          (let ((sqrt-1-z (complex-sqrt (1-z z)))
1541                (sqrt-1+z (complex-sqrt (1+z z))))
1542            (with-float-traps-masked (:divide-by-zero)
1543              (complex (atan (/ (realpart z)
1544                                (realpart (* sqrt-1-z sqrt-1+z))))
1545                       (asinh (imagpart (* (conjugate sqrt-1-z)
1546                                           sqrt-1+z))))))))
1547    
1548    (defun complex-asinh (z)
1549      "Compute asinh z = log(z + sqrt(1 + z*z))
1550    
1551    Z may be any number, but the result is always a complex."
1552      (declare (number z))
1553      ;; asinh z = -i * asin (i*z)
1554      #+double-double
1555      (when (typep z '(or double-double-float (complex double-double-float)))
1556        (return-from complex-asinh (dd-complex-asinh z)))
1557      (let* ((iz (complex (- (imagpart z)) (realpart z)))
1558             (result (complex-asin iz)))
1559        (complex (imagpart result)
1560                 (- (realpart result)))))
1561    
1562    (defun complex-atan (z)
1563      "Compute atan z = atanh (i*z) / i
1564    
1565    Z may be any number, but the result is always a complex."
1566      (declare (number z))
1567      ;; atan z = -i * atanh (i*z)
1568      #+double-double
1569      (when (typep z '(or double-double-float (complex double-double-float)))
1570        (return-from complex-atan (dd-complex-atan z)))
1571      (let* ((iz (complex (- (imagpart z)) (realpart z)))
1572             (result (complex-atanh iz)))
1573        (complex (imagpart result)
1574                 (- (realpart result)))))
1575    
1576    (defun complex-tan (z)
1577      "Compute tan z = -i * tanh(i * z)
1578    
1579    Z may be any number, but the result is always a complex."
1580      (declare (number z))
1581      ;; tan z = -i * tanh(i*z)
1582      #+double-double
1583      (when (typep z '(or double-double-float (complex double-double-float)))
1584        (return-from complex-tan (dd-complex-tan z)))
1585      (let* ((iz (complex (- (imagpart z)) (realpart z)))
1586             (result (complex-tanh iz)))
1587        (complex (imagpart result)
1588                 (- (realpart result)))))
1589    

Legend:
Removed from v.1.9.1.1  
changed lines
  Added in v.1.64

  ViewVC Help
Powered by ViewVC 1.1.5