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

Contents of /src/code/irrat.lisp

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.5