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

Contents of /src/code/irrat.lisp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.47 - (hide annotations)
Wed Jul 19 02:54:30 2006 UTC (7 years, 9 months ago) by rtoy
Branch: MAIN
Changes since 1.46: +47 -1 lines
Add the trig argument reduction routines from Sun's fdlibm so we can
accurately reduce the arg and therefore compute the value of trig
functions accurately.

lisp/Config.linux_gencgc:
o Compile e_rem_pio2.c and k_rem_pio2.c

code/irrat.lisp:
o Disable %sin, %cos, %tan functions.
o Implement %sin, %cos, and %tan to call the fdlibm routine
  __ieee754_rem_pio2 to do argument reduction before calling the sin,
  cos, tan vops.

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

  ViewVC Help
Powered by ViewVC 1.1.5