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

Contents of /src/code/irrat.lisp

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.5