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

Contents of /src/code/irrat.lisp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.46 - (hide annotations)
Fri Jun 30 18:41:22 2006 UTC (7 years, 9 months ago) by rtoy
Branch: MAIN
CVS Tags: snapshot-2006-07
Changes since 1.45: +136 -16 lines
This large checkin merges the double-double float support to HEAD.
The merge is from the tag "double-double-irrat-end".  The
double-double branch is now obsolete.

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

  ViewVC Help
Powered by ViewVC 1.1.5