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

Contents of /src/code/irrat.lisp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.54 - (hide annotations)
Mon Jan 28 18:21:03 2008 UTC (6 years, 2 months ago) by rtoy
Branch: MAIN
Changes since 1.53: +37 -23 lines
Add new interface to ieee754_rem_pio2.  No longer need to pass in an
array. The new function returns 2 new output values.

code/irrat.lisp:
o Rename the original %ieee754-rem-pi/2 to %%ieee754-rem-pi-2.
o Define the new %ieee754-rem-pi/2 function.  This returns the output
  as two output values instead of an array.
o Use the new function.  We should have wrapped the original with
  without-gcing, but we don't have to anymore.

lisp/ppc-arch.c:
lisp/x86-arch.c:
o Implement the new C interface to __ieee754_rem_pio2

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

  ViewVC Help
Powered by ViewVC 1.1.5