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

Contents of /src/code/irrat.lisp

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.5