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

Contents of /src/code/irrat.lisp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.64 - (hide annotations)
Mon Feb 28 17:14:45 2011 UTC (3 years, 1 month ago) by rtoy
Branch: MAIN
CVS Tags: GIT-CONVERSION, snapshot-2011-09, snapshot-2011-06, snapshot-2011-07, snapshot-2011-04, snapshot-2011-03, HEAD
Changes since 1.63: +2 -2 lines
Need to define %sqrt for x86/sse2 so the compiler can constant-fold
SQRT.
1 wlott 1.1 ;;; -*- Mode: Lisp; Package: KERNEL; Log: code.log -*-
2     ;;;
3     ;;; **********************************************************************
4 ram 1.7 ;;; 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.
6     ;;;
7     (ext:file-comment
8 rtoy 1.64 "$Header: /tiger/var/lib/cvsroots/cmucl/src/code/irrat.lisp,v 1.64 2011/02/28 17:14:45 rtoy Exp $")
9 ram 1.7 ;;;
10 wlott 1.1 ;;; **********************************************************************
11     ;;;
12     ;;; This file contains all the irrational functions. Actually, most of the
13     ;;; work is done by calling out to C...
14     ;;;
15     ;;; Author: William Lott.
16     ;;;
17    
18     (in-package "KERNEL")
19 rtoy 1.61 (intl:textdomain "cmucl")
20 wlott 1.1
21    
22     ;;;; Random constants, utility functions, and macros.
23    
24     (defconstant pi 3.14159265358979323846264338327950288419716939937511L0)
25 wlott 1.2 ;(defconstant e 2.71828182845904523536028747135266249775724709369996L0)
26 wlott 1.1
27 ram 1.5 ;;; Make these INLINE, since the call to C is at least as compact as a Lisp
28     ;;; call, and saves number consing to boot.
29     ;;;
30 wlott 1.1 (defmacro def-math-rtn (name num-args)
31     (let ((function (intern (concatenate 'simple-string
32     "%"
33     (string-upcase name)))))
34 ram 1.4 `(progn
35 pw 1.31 (declaim (inline ,function))
36 ram 1.4 (export ',function)
37 wlott 1.10 (alien:def-alien-routine (,name ,function) double-float
38 ram 1.4 ,@(let ((results nil))
39     (dotimes (i num-args (nreverse results))
40     (push (list (intern (format nil "ARG-~D" i))
41     'double-float)
42     results)))))))
43 wlott 1.1
44     (eval-when (compile load eval)
45    
46     (defun handle-reals (function var)
47     `((((foreach fixnum single-float bignum ratio))
48     (coerce (,function (coerce ,var 'double-float)) 'single-float))
49     ((double-float)
50 rtoy 1.46 (,function ,var))
51     #+double-double
52     ((double-double-float)
53     (,(symbolicate "DD-" function) ,var))))
54 wlott 1.1
55     ); eval-when (compile load eval)
56    
57    
58     ;;;; Stubs for the Unix math library.
59    
60     ;;; Please refer to the Unix man pages for details about these routines.
61    
62     ;;; Trigonometric.
63 rtoy 1.58 #-(and x86 (not sse2))
64     (progn
65     ;; 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 wlott 1.1 (def-math-rtn "asin" 1)
74     (def-math-rtn "acos" 1)
75     (def-math-rtn "sinh" 1)
76     (def-math-rtn "cosh" 1)
77     (def-math-rtn "tanh" 1)
78 pw 1.18 (def-math-rtn "asinh" 1)
79     (def-math-rtn "acosh" 1)
80     (def-math-rtn "atanh" 1)
81 wlott 1.1
82     ;;; Exponential and Logarithmic.
83 rtoy 1.58 #-(and x86 (not sse2))
84     (progn
85     (def-math-rtn "exp" 1)
86     (def-math-rtn "log" 1)
87     (def-math-rtn "log10" 1))
88    
89 wlott 1.1 (def-math-rtn "pow" 2)
90 rtoy 1.58 #-(or x86 sparc-v7 sparc-v8 sparc-v9)
91     (def-math-rtn "sqrt" 1)
92 wlott 1.1 (def-math-rtn "hypot" 2)
93 ram 1.17
94 rtoy 1.58 ;; 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 ram 1.17 (progn
103 rtoy 1.47 #+nil
104 ram 1.17 (defun %sin (x)
105 dtc 1.19 (declare (double-float x)
106 ram 1.17 (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 rtoy 1.47 #+nil
113 ram 1.17 (defun %cos (x)
114 dtc 1.19 (declare (double-float x)
115 ram 1.17 (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 rtoy 1.47 #+nil
122 ram 1.17 (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 dtc 1.19 (declare (double-float x)
140 ram 1.17 (values double-float))
141     (%exp x))
142     (defun %log (x)
143 dtc 1.19 (declare (double-float x)
144 ram 1.17 (values double-float))
145     (%log x))
146     (defun %log10 (x)
147 dtc 1.19 (declare (double-float x)
148 ram 1.17 (values double-float))
149     (%log10 x))
150     #+nil ;; notyet
151     (defun %pow (x y)
152     (declare (type (double-float 0d0) x)
153 dtc 1.19 (double-float y)
154 ram 1.17 (values (double-float 0d0)))
155     (%pow x y))
156     (defun %sqrt (x)
157 dtc 1.19 (declare (double-float x)
158 ram 1.17 (values double-float))
159     (%sqrt x))
160     (defun %scalbn (f ex)
161 dtc 1.19 (declare (double-float f)
162 ram 1.17 (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 dtc 1.19 (declare (double-float x)
171     (values double-float))
172     (%logb x))
173 dtc 1.27 (defun %log1p (x)
174     (declare (double-float x)
175     (values double-float))
176     (%log1p x))
177 ram 1.17 ) ; progn
178 wlott 1.1
179 rtoy 1.44
180     ;; As above for x86. It also seems to be needed to handle
181     ;; constant-folding in the compiler.
182 rtoy 1.64 #+(or sparc (and x86 sse2))
183 rtoy 1.44 (progn
184     (defun %sqrt (x)
185     (declare (double-float x)
186     (values double-float))
187     (%sqrt x))
188     )
189    
190 rtoy 1.49 ;;; 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 rtoy 1.48 #+(or ppc x86)
200 rtoy 1.47 (progn
201 rtoy 1.54 (declaim (inline %%ieee754-rem-pi/2))
202 rtoy 1.49 ;; 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 rtoy 1.54 (alien:def-alien-routine ("__ieee754_rem_pio2" %%ieee754-rem-pi/2) c-call:int
208 rtoy 1.47 (x double-float)
209     (y (* double-float)))
210 rtoy 1.54
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 rtoy 1.48 )
222 rtoy 1.47
223 rtoy 1.58 #+(or ppc sse2)
224 rtoy 1.48 (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 rtoy 1.47 )
233    
234 rtoy 1.48 #+(or ppc x86)
235     (macrolet
236     ((frob (sin cos tan)
237 rtoy 1.54 `(progn
238 rtoy 1.55 ;; 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 rtoy 1.49 ;;
246 rtoy 1.55 ;; For x86 with an fsin/fcos/fptan instruction, the pi/4 is
247     ;; probably too restrictive.
248 rtoy 1.48 (defun %sin (x)
249     (declare (double-float x))
250     (if (< (abs x) (/ pi 4))
251     (,sin x)
252     ;; Argument reduction needed
253 rtoy 1.54 (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 rtoy 1.48 (defun %cos (x)
262     (declare (double-float x))
263     (if (< (abs x) (/ pi 4))
264     (,cos x)
265     ;; Argument reduction needed
266 rtoy 1.54 (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 rtoy 1.48 (defun %tan (x)
275     (declare (double-float x))
276     (if (< (abs x) (/ pi 4))
277     (,tan x)
278     ;; Argument reduction needed
279 rtoy 1.54 (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 rtoy 1.58 ;; Don't want %sin-quick and friends with sse2.
286     #+(and x86 (not sse2))
287 rtoy 1.48 (frob %sin-quick %cos-quick %tan-quick)
288 rtoy 1.58 #+(or ppc sse2)
289 rtoy 1.48 (frob %%sin %%cos %%tan))
290    
291    
292 wlott 1.1
293     ;;;; Power functions.
294    
295     (defun exp (number)
296 rtoy 1.62 "Return e raised to the power NUMBER."
297 wlott 1.1 (number-dispatch ((number number))
298     (handle-reals %exp number)
299     ((complex)
300     (* (exp (realpart number))
301     (cis (imagpart number))))))
302    
303     ;;; INTEXP -- Handle the rational base, integer power case.
304    
305     (defparameter *intexp-maximum-exponent* 10000)
306    
307 rtoy 1.57 (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 rtoy 1.63 (format stream (intl:gettext "The absolute value of ~S exceeds limit ~S.")
312 rtoy 1.57 (intexp-power condition)
313     *intexp-maximum-exponent*))))
314    
315 ram 1.6 ;;; 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
317     ;;; can be calculated more efficiently if power is a positive integer. Values
318     ;;; of power are calculated as positive integers, and inverted if negative.
319     ;;;
320 wlott 1.1 (defun intexp (base power)
321 rtoy 1.45 ;; 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 wlott 1.1 (when (> (abs power) *intexp-maximum-exponent*)
328 rtoy 1.57 ;; Allow user the option to continue with calculation, possibly
329     ;; increasing the limit to the given power.
330     (restart-case
331     (error 'intexp-limit-error
332     :base base
333     :power power)
334     (continue ()
335 rtoy 1.61 :report (lambda (stream)
336 rtoy 1.63 (write-string (intl:gettext "Continue with calculation") stream)))
337 rtoy 1.57 (new-limit ()
338 rtoy 1.61 :report (lambda (stream)
339 rtoy 1.63 (write-string (intl:gettext "Continue with calculation, update limit") stream))
340 rtoy 1.60 (setq *intexp-maximum-exponent* (abs power)))))
341 wlott 1.1 (cond ((minusp power)
342     (/ (intexp base (- power))))
343     ((eql base 2)
344     (ash 1 power))
345     (t
346     (do ((nextn (ash power -1) (ash power -1))
347     (total (if (oddp power) base 1)
348     (if (oddp power) (* base total) total)))
349     ((zerop nextn) total)
350     (setq base (* base base))
351     (setq power nextn)))))
352    
353    
354 ram 1.6 ;;; EXPT -- Public
355     ;;;
356     ;;; If an integer power of a rational, use INTEXP above. Otherwise, do
357     ;;; floating point stuff. If both args are real, we try %POW right off,
358     ;;; assuming it will return 0 if the result may be complex. If so, we call
359     ;;; COMPLEX-POW which directly computes the complex result. We also separate
360     ;;; the complex-real and real-complex cases from the general complex case.
361     ;;;
362 wlott 1.1 (defun expt (base power)
363 rtoy 1.62 "Returns BASE raised to the POWER."
364 wlott 1.1 (if (zerop power)
365 rtoy 1.40 ;; CLHS says that if the power is 0, the result is 1, subject to
366     ;; numeric contagion. But what happens if base is infinity or
367     ;; NaN? Do we silently return 1? For now, I think we should
368     ;; signal an error if the FP modes say so.
369     (let ((result (1+ (* base power))))
370     ;; If we get an NaN here, that means base*power above didn't
371     ;; produce 0 and FP traps were disabled, so we handle that
372     ;; here. Should this be a continuable restart?
373     (if (and (floatp result) (float-nan-p result))
374     (float 1 result)
375     result))
376 dtc 1.22 (labels (;; determine if the double float is an integer.
377     ;; 0 - not an integer
378     ;; 1 - an odd int
379     ;; 2 - an even int
380     (isint (ihi lo)
381     (declare (type (unsigned-byte 31) ihi)
382     (type (unsigned-byte 32) lo)
383     (optimize (speed 3) (safety 0)))
384     (let ((isint 0))
385     (declare (type fixnum isint))
386     (cond ((>= ihi #x43400000) ; exponent >= 53
387     (setq isint 2))
388     ((>= ihi #x3ff00000)
389     (let ((k (- (ash ihi -20) #x3ff))) ; exponent
390     (declare (type (mod 53) k))
391     (cond ((> k 20)
392     (let* ((shift (- 52 k))
393     (j (logand (ash lo (- shift))))
394     (j2 (ash j shift)))
395     (declare (type (mod 32) shift)
396     (type (unsigned-byte 32) j j2))
397     (when (= j2 lo)
398     (setq isint (- 2 (logand j 1))))))
399     ((= lo 0)
400     (let* ((shift (- 20 k))
401     (j (ash ihi (- shift)))
402     (j2 (ash j shift)))
403     (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 dtc 1.24 (complex
468 dtc 1.22 (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 dtc 1.24 (complex
490 dtc 1.22 (coerce (* pow (%cos y*pi)) rtype)
491     (coerce (* pow (%sin y*pi)) rtype)))))))))))))
492 rtoy 1.56 ;; 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 dtc 1.22 (number-dispatch ((base number) (power number))
500 rtoy 1.56 (((foreach fixnum (or bignum ratio) (complex rational))
501     integer)
502 dtc 1.22 (intexp base power))
503 rtoy 1.56 (((foreach single-float double-float)
504     rational)
505 dtc 1.22 (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 rtoy 1.46 #+double-double
515 rtoy 1.56 (((foreach fixnum (or bignum ratio) single-float double-float
516     double-double-float)
517 rtoy 1.46 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 rtoy 1.56 (((foreach (complex rational) (complex single-float) (complex double-float)
524     #+double-double (complex double-double-float))
525     rational)
526 dtc 1.22 (* (expt (abs base) power)
527     (cis (* power (phase base)))))
528 rtoy 1.56 #+double-double
529     ((double-double-float
530 dtc 1.22 complex)
531     (if (and (zerop base) (plusp (realpart power)))
532     (* base power)
533 rtoy 1.56 (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 dtc 1.22 (exp (* power (log base)))))
581 rtoy 1.56 #+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 dtc 1.22 (if (and (zerop base) (plusp (realpart power)))
591     (* base power)
592     (exp (* power (log base)))))))))
593 wlott 1.1
594 rtoy 1.52 ;; 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 rtoy 1.51 #.(log 2d0))
603 rtoy 1.52 #+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 rtoy 1.51 ;; 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 rtoy 1.52 (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 rtoy 1.51 bignum)))
644 rtoy 1.52 (values n (log (scale-float (float f float-type) (- exp))
645     (float 2 float-type)))))))))
646 rtoy 1.51 (etypecase x
647     (float
648 rtoy 1.52 (/ (log (float x float-type)) (log-of-2 float-type)))
649 rtoy 1.51 (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 rtoy 1.52 (/ (log1p (float (- x 1) float-type))
666     (log-of-2 float-type))
667 rtoy 1.51 (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 rtoy 1.43
678 wlott 1.1 (defun log (number &optional (base nil base-p))
679 rtoy 1.62 "Return the logarithm of NUMBER in the base BASE, which defaults to e."
680 wlott 1.1 (if base-p
681 toy 1.34 (cond ((zerop base)
682     ;; ANSI spec
683     base)
684 rtoy 1.43 ((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 rtoy 1.51 (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 rtoy 1.43 double-float)
709 rtoy 1.51 (/ (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 rtoy 1.43 (/ (log (coerce number 'double-float)) (log base)))
722 rtoy 1.46 #+double-double
723     ((double-double-float
724 rtoy 1.51 (foreach fixnum bignum ratio))
725 rtoy 1.52 (/ (log2 number 1w0) (log2 base 1w0)))
726 rtoy 1.51 #+double-double
727     ((double-double-float
728     (foreach double-double-float double-float single-float))
729 rtoy 1.46 (/ (log number) (log (coerce base 'double-double-float))))
730     #+double-double
731 rtoy 1.51 (((foreach fixnum bignum ratio)
732     double-double-float)
733 rtoy 1.52 (/ (log2 number 1w0) (log2 base 1w0)))
734 rtoy 1.51 #+double-double
735     (((foreach double-float single-float)
736 rtoy 1.46 double-double-float)
737     (/ (log (coerce number 'double-double-float)) (log base)))
738 rtoy 1.51 (((foreach single-float)
739     (foreach single-float))
740 rtoy 1.43 ;; 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 toy 1.34 (t
754 rtoy 1.43 ;; FIXME: This probably needs some work as well.
755 toy 1.34 (/ (log number) (log base))))
756 wlott 1.1 (number-dispatch ((number number))
757 toy 1.34 (((foreach fixnum bignum))
758 wlott 1.3 (if (minusp number)
759 toy 1.34 (complex (coerce (log (- number)) 'single-float)
760     (coerce pi 'single-float))
761 toy 1.37 (coerce (/ (log2 number) #.(log (exp 1d0) 2d0)) 'single-float)))
762 toy 1.34 ((ratio)
763     (if (minusp number)
764     (complex (coerce (log (- number)) 'single-float)
765     (coerce pi 'single-float))
766 toy 1.35 ;; 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 toy 1.37 ;; 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 toy 1.35 (if (= (integer-length top)
782 toy 1.36 (integer-length bot))
783 toy 1.35 (coerce (%log1p (coerce (- number 1) 'double-float))
784     'single-float)
785 toy 1.37 (coerce (/ (- (log2 top) (log2 bot))
786     #.(log (exp 1d0) 2d0))
787     'single-float)))))
788 ram 1.17 (((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 rtoy 1.46 #+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 ram 1.17 ((complex)
804     (complex-log number)))))
805 wlott 1.1
806     (defun sqrt (number)
807 rtoy 1.62 "Return the square root of NUMBER."
808 wlott 1.1 (number-dispatch ((number number))
809 ram 1.17 (((foreach fixnum bignum ratio))
810 wlott 1.3 (if (minusp number)
811 ram 1.17 (complex-sqrt number)
812 wlott 1.3 (coerce (%sqrt (coerce number 'double-float)) 'single-float)))
813 ram 1.17 (((foreach single-float double-float))
814 dtc 1.28 (if (minusp number)
815 ram 1.17 (complex-sqrt number)
816     (coerce (%sqrt (coerce number 'double-float))
817     '(dispatch-type number))))
818 rtoy 1.46 #+double-double
819     ((double-double-float)
820     (if (minusp number)
821     (dd-complex-sqrt number)
822     (multiple-value-bind (hi lo)
823     (c::sqrt-dd (kernel:double-double-hi number) (kernel:double-double-lo number))
824 rtoy 1.50 (kernel:%make-double-double-float hi lo))))
825 rtoy 1.46 ((complex)
826     (complex-sqrt number))))
827 wlott 1.1
828    
829     ;;;; Trigonometic and Related Functions
830    
831 wlott 1.2 (defun abs (number)
832 rtoy 1.62 "Returns the absolute value of the number."
833 wlott 1.2 (number-dispatch ((number number))
834 rtoy 1.46 (((foreach single-float double-float fixnum rational
835     #+double-double double-double-float))
836 wlott 1.2 (abs number))
837     ((complex)
838     (let ((rx (realpart number))
839     (ix (imagpart number)))
840     (etypecase rx
841     (rational
842     (sqrt (+ (* rx rx) (* ix ix))))
843     (single-float
844     (coerce (%hypot (coerce rx 'double-float)
845     (coerce ix 'double-float))
846     'single-float))
847     (double-float
848 rtoy 1.46 (%hypot rx ix))
849     #+double-double
850     (double-double-float
851 rtoy 1.50 (multiple-value-bind (abs^2 scale)
852     (dd-cssqs number)
853     (scale-float (sqrt abs^2) scale))))))))
854 wlott 1.1
855     (defun phase (number)
856 rtoy 1.62 "Returns the angle part of the polar representation of a complex number.
857 wlott 1.3 For complex numbers, this is (atan (imagpart number) (realpart number)).
858 wlott 1.1 For non-complex positive numbers, this is 0. For non-complex negative
859     numbers this is PI."
860 wlott 1.3 (etypecase number
861 ram 1.17 (rational
862 wlott 1.3 (if (minusp number)
863     (coerce pi 'single-float)
864     0.0f0))
865 ram 1.17 (single-float
866     (if (minusp (float-sign number))
867     (coerce pi 'single-float)
868     0.0f0))
869 wlott 1.3 (double-float
870 ram 1.17 (if (minusp (float-sign number))
871 wlott 1.3 (coerce pi 'double-float)
872     0.0d0))
873 rtoy 1.46 #+double-double
874     (double-double-float
875     (if (minusp (float-sign number))
876     dd-pi
877     0w0))
878 wlott 1.3 (complex
879     (atan (imagpart number) (realpart number)))))
880 wlott 1.1
881    
882     (defun sin (number)
883 rtoy 1.62 "Return the sine of NUMBER."
884 wlott 1.1 (number-dispatch ((number number))
885     (handle-reals %sin number)
886     ((complex)
887     (let ((x (realpart number))
888     (y (imagpart number)))
889 ram 1.17 (complex (* (sin x) (cosh y))
890     (* (cos x) (sinh y)))))))
891 wlott 1.1
892     (defun cos (number)
893 rtoy 1.62 "Return the cosine of NUMBER."
894 wlott 1.1 (number-dispatch ((number number))
895     (handle-reals %cos number)
896     ((complex)
897     (let ((x (realpart number))
898     (y (imagpart number)))
899 ram 1.17 (complex (* (cos x) (cosh y))
900     (- (* (sin x) (sinh y))))))))
901 wlott 1.1
902     (defun tan (number)
903 rtoy 1.62 "Return the tangent of NUMBER."
904 wlott 1.1 (number-dispatch ((number number))
905     (handle-reals %tan number)
906     ((complex)
907 ram 1.17 (complex-tan number))))
908 wlott 1.1
909 wlott 1.2 (defun cis (theta)
910 rtoy 1.62 "Return cos(Theta) + i sin(Theta), AKA exp(i Theta)."
911 wlott 1.2 (if (complexp theta)
912 rtoy 1.63 (error (intl:gettext "Argument to CIS is complex: ~S") theta)
913 wlott 1.2 (complex (cos theta) (sin theta))))
914 wlott 1.1
915     (defun asin (number)
916 rtoy 1.62 "Return the arc sine of NUMBER."
917 wlott 1.1 (number-dispatch ((number number))
918 wlott 1.3 ((rational)
919     (if (or (> number 1) (< number -1))
920     (complex-asin number)
921     (coerce (%asin (coerce number 'double-float)) 'single-float)))
922     (((foreach single-float double-float))
923 rtoy 1.42 (if (or (float-nan-p number)
924     (and (<= number (coerce 1 '(dispatch-type number)))
925     (>= number (coerce -1 '(dispatch-type number)))))
926 wlott 1.3 (coerce (%asin (coerce number 'double-float))
927 rtoy 1.42 '(dispatch-type number))
928 rtoy 1.46 (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 wlott 1.1 ((complex)
937 wlott 1.3 (complex-asin number))))
938 wlott 1.1
939     (defun acos (number)
940 rtoy 1.62 "Return the arc cosine of NUMBER."
941 wlott 1.1 (number-dispatch ((number number))
942 wlott 1.3 ((rational)
943     (if (or (> number 1) (< number -1))
944     (complex-acos number)
945     (coerce (%acos (coerce number 'double-float)) 'single-float)))
946     (((foreach single-float double-float))
947 rtoy 1.42 (if (or (float-nan-p number)
948     (and (<= number (coerce 1 '(dispatch-type number)))
949     (>= number (coerce -1 '(dispatch-type number)))))
950 wlott 1.3 (coerce (%acos (coerce number 'double-float))
951 rtoy 1.42 '(dispatch-type number))
952 rtoy 1.46 (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 wlott 1.1 ((complex)
961 wlott 1.3 (complex-acos number))))
962 wlott 1.1
963    
964     (defun atan (y &optional (x nil xp))
965 rtoy 1.62 "Return the arc tangent of Y if X is omitted or Y/X if X is supplied."
966 wlott 1.1 (if xp
967 wlott 1.12 (flet ((atan2 (y x)
968 wlott 1.13 (declare (type double-float y x)
969     (values double-float))
970     (if (zerop x)
971     (if (zerop y)
972     (if (plusp (float-sign x))
973     y
974     (float-sign y pi))
975     (float-sign y (/ pi 2)))
976     (%atan2 y x))))
977 toy 1.33 ;; If X is given, both X and Y must be real numbers.
978     (number-dispatch ((y real) (x real))
979 wlott 1.13 ((double-float
980     (foreach double-float single-float fixnum bignum ratio))
981     (atan2 y (coerce x 'double-float)))
982     (((foreach single-float fixnum bignum ratio)
983     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 rtoy 1.46 '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 wlott 1.1 (number-dispatch ((y number))
998     (handle-reals %atan y)
999     ((complex)
1000 ram 1.17 (complex-atan y)))))
1001 wlott 1.3
1002 ram 1.17 (defun sinh (number)
1003 rtoy 1.62 "Return the hyperbolic sine of NUMBER."
1004 ram 1.17 (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)
1013 rtoy 1.62 "Return the hyperbolic cosine of NUMBER."
1014 ram 1.17 (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 wlott 1.3 (defun tanh (number)
1023 rtoy 1.62 "Return the hyperbolic tangent of NUMBER."
1024 ram 1.17 (number-dispatch ((number number))
1025     (handle-reals %tanh number)
1026     ((complex)
1027     (complex-tanh number))))
1028 wlott 1.3
1029     (defun asinh (number)
1030 rtoy 1.62 "Return the hyperbolic arc sine of NUMBER."
1031 ram 1.17 (number-dispatch ((number number))
1032     (handle-reals %asinh number)
1033     ((complex)
1034     (complex-asinh number))))
1035 wlott 1.3
1036     (defun acosh (number)
1037 rtoy 1.62 "Return the hyperbolic arc cosine of NUMBER."
1038 ram 1.17 (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 rtoy 1.46 #+double-double
1050     ((double-double-float)
1051     (if (< number 1w0)
1052     (complex-acosh number)
1053     (dd-%acosh number)))
1054 ram 1.17 ((complex)
1055     (complex-acosh number))))
1056 wlott 1.3
1057     (defun atanh (number)
1058 rtoy 1.62 "Return the hyperbolic arc tangent of NUMBER."
1059 ram 1.17 (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 rtoy 1.46 #+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 ram 1.17 ((complex)
1078     (complex-atanh number))))
1079 wlott 1.14
1080 toy 1.32 ;;; 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 wlott 1.14
1084     #+hpux
1085 pw 1.18 (declaim (inline %log1p))
1086 wlott 1.14 #+hpux
1087 pw 1.18 (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 ram 1.17
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 toy 1.32
1144 ram 1.17 (declaim (inline square))
1145     (defun square (x)
1146 rtoy 1.46 (declare (float x))
1147 ram 1.17 (* 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 dtc 1.19 (declaim (inline scalb))
1154 ram 1.17 (defun scalb (x n)
1155 rtoy 1.62 "Compute 2^N * X without compute 2^N first (use properties of the
1156 ram 1.17 underlying floating-point format"
1157 rtoy 1.46 (declare (type float x)
1158 dtc 1.19 (type double-float-exponent n))
1159 ram 1.17 (scale-float x n))
1160    
1161 toy 1.32 (declaim (inline logb-finite))
1162     (defun logb-finite (x)
1163 rtoy 1.62 "Same as logb but X is not infinity and non-zero and not a NaN, so
1164 toy 1.32 that we can always return an integer"
1165 rtoy 1.46 (declare (type float x))
1166 toy 1.32 (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 ram 1.17 (defun logb (x)
1174 rtoy 1.62 "Compute an integer N such that 1 <= |2^(-N) * x| < 2.
1175 ram 1.17 For the special cases, the following values are used:
1176    
1177     x logb
1178     NaN NaN
1179     +/- infinity +infinity
1180     0 -infinity
1181     "
1182 rtoy 1.46 (declare (type float x))
1183 dtc 1.26 (cond ((float-nan-p x)
1184 ram 1.17 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 toy 1.32 ;; signal divide-by-zero, so do the actual division
1190 rtoy 1.46 (/ -1 x)
1191 ram 1.17 )
1192     (t
1193 toy 1.32 (logb-finite x))))
1194    
1195    
1196 ram 1.17
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 rtoy 1.62 "Create complex number with real part X and imaginary part Y such that
1203 ram 1.17 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 toy 1.32 (number z)
1207     (optimize (extensions:inhibit-warnings 3)))
1208     (if (typep (realpart z) 'double-float)
1209 ram 1.17 (complex x y)
1210     ;; Convert anything that's not a double-float to a single-float.
1211 toy 1.32 (complex (float x 1f0)
1212     (float y 1f0))))
1213 ram 1.17
1214     (defun cssqs (z)
1215 dtc 1.21 ;; 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 ram 1.17
1218     ;; Save all FP flags
1219 dtc 1.21 (let ((x (float (realpart z) 1d0))
1220 toy 1.32 (y (float (imagpart z) 1d0)))
1221 dtc 1.21 ;; 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 toy 1.32 (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 ram 1.17
1252     (defun complex-sqrt (z)
1253 rtoy 1.62 "Principle square root of Z
1254 ram 1.17
1255     Z may be any number, but the result is always a complex."
1256     (declare (number z))
1257 rtoy 1.46 #+double-double
1258     (when (typep z '(or double-double-float (complex double-double-float)))
1259     (return-from complex-sqrt (dd-complex-sqrt z)))
1260 ram 1.17 (multiple-value-bind (rho k)
1261     (cssqs z)
1262 toy 1.32 (declare (type (or (member 0d0) (double-float 0d0)) rho)
1263     (type fixnum k))
1264 ram 1.17 (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 toy 1.32 (locally
1271     ;; space 0 to get maybe-inline functions inlined.
1272     (declare (optimize (speed 3) (space 0)))
1273    
1274 rtoy 1.53 (if (not (locally (declare (optimize (inhibit-warnings 3)))
1275     (float-nan-p x)))
1276 toy 1.32 (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 ram 1.17
1297     (defun complex-log-scaled (z j)
1298 rtoy 1.62 "Compute log(2^j*z).
1299 ram 1.17
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 toy 1.32 (declare (optimize (speed 3)))
1317 ram 1.17 (let ((beta (max (abs x) (abs y)))
1318     (theta (min (abs x) (abs y))))
1319 toy 1.32 (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 ram 1.17
1332     (defun complex-log (z)
1333 rtoy 1.62 "Log of Z = log |Z| + i * arg Z
1334 ram 1.17
1335     Z may be any number, but the result is always a complex."
1336     (declare (number z))
1337 rtoy 1.46 #+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 ram 1.17 (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 rtoy 1.62 "Compute atanh z = (log(1+z) - log(1-z))/2"
1349 ram 1.17 (declare (number z))
1350 rtoy 1.46 #+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 rtoy 1.41 (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 ram 1.17
1411     (defun complex-tanh (z)
1412 rtoy 1.62 "Compute tanh z = sinh z / cosh z"
1413 ram 1.17 (declare (number z))
1414 rtoy 1.46 #+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 ram 1.17 (let ((x (float (realpart z) 1.0d0))
1419     (y (float (imagpart z) 1.0d0)))
1420 toy 1.32 (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 ram 1.17
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 toy 1.35 ;; and these two expressions are equal if and only if arg conj z =
1465     ;; -arg z, which is clearly true for all z.
1466 ram 1.17
1467 rtoy 1.38 ;; 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 ram 1.17 (defun complex-acos (z)
1497 rtoy 1.62 "Compute acos z = pi/2 - asin z
1498 ram 1.17
1499     Z may be any number, but the result is always a complex."
1500     (declare (number z))
1501 rtoy 1.46 #+double-double
1502     (when (typep z '(or double-double-float (complex double-double-float)))
1503     (return-from complex-acos (dd-complex-acos z)))
1504 rtoy 1.41 (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 ram 1.17
1515     (defun complex-acosh (z)
1516 rtoy 1.62 "Compute acosh z = 2 * log(sqrt((z+1)/2) + sqrt((z-1)/2))
1517 ram 1.17
1518     Z may be any number, but the result is always a complex."
1519     (declare (number z))
1520 rtoy 1.38 (let ((sqrt-z-1 (complex-sqrt (z-1 z)))
1521     (sqrt-z+1 (complex-sqrt (z+1 z))))
1522 dtc 1.21 (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 ram 1.17
1528    
1529     (defun complex-asin (z)
1530 rtoy 1.62 "Compute asin z = asinh(i*z)/i
1531 ram 1.17
1532     Z may be any number, but the result is always a complex."
1533     (declare (number z))
1534 rtoy 1.46 #+double-double
1535     (when (typep z '(or double-double-float (complex double-double-float)))
1536     (return-from complex-asin (dd-complex-asin z)))
1537 rtoy 1.41 (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 ram 1.17
1548     (defun complex-asinh (z)
1549 rtoy 1.62 "Compute asinh z = log(z + sqrt(1 + z*z))
1550 ram 1.17
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 rtoy 1.46 #+double-double
1555     (when (typep z '(or double-double-float (complex double-double-float)))
1556     (return-from complex-asinh (dd-complex-asinh z)))
1557 ram 1.17 (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 rtoy 1.62 "Compute atan z = atanh (i*z) / i
1564 ram 1.17
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 rtoy 1.46 #+double-double
1569     (when (typep z '(or double-double-float (complex double-double-float)))
1570     (return-from complex-atan (dd-complex-atan z)))
1571 ram 1.17 (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 rtoy 1.62 "Compute tan z = -i * tanh(i * z)
1578 ram 1.17
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 rtoy 1.46 #+double-double
1583     (when (typep z '(or double-double-float (complex double-double-float)))
1584     (return-from complex-tan (dd-complex-tan z)))
1585 ram 1.17 (let* ((iz (complex (- (imagpart z)) (realpart z)))
1586     (result (complex-tanh iz)))
1587     (complex (imagpart result)
1588     (- (realpart result)))))
1589 toy 1.32

  ViewVC Help
Powered by ViewVC 1.1.5