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

Contents of /src/code/irrat.lisp

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.5