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

Contents of /src/code/irrat.lisp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.29 - (hide annotations)
Wed Jun 23 14:05:17 1999 UTC (14 years, 9 months ago) by dtc
Branch: MAIN
Changes since 1.28: +2 -2 lines
Use the fsqrt VOP on the SPARC V7 and up.
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 dtc 1.29 "$Header: /tiger/var/lib/cvsroots/cmucl/src/code/irrat.lisp,v 1.29 1999/06/23 14:05:17 dtc 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 ram 1.5 (proclaim '(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.29 #-(or x86 sparc-v7) (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 ;; It seems that everyone has a C version of sinh, cosh, and
528     ;; tanh. Let's use these for reals because the original
529     ;; implementations based on the definitions lose big in round-off
530     ;; error. These bad definitions also mean that sin and cos for
531     ;; complex numbers can also lose big.
532 wlott 1.3
533 ram 1.17 #+nil
534 wlott 1.3 (defun sinh (number)
535     "Return the hyperbolic sine of NUMBER."
536     (/ (- (exp number) (exp (- number))) 2))
537    
538 ram 1.17 (defun sinh (number)
539     "Return the hyperbolic sine of NUMBER."
540     (number-dispatch ((number number))
541     (handle-reals %sinh number)
542     ((complex)
543     (let ((x (realpart number))
544     (y (imagpart number)))
545     (complex (* (sinh x) (cos y))
546     (* (cosh x) (sin y)))))))
547    
548     #+nil
549 wlott 1.3 (defun cosh (number)
550     "Return the hyperbolic cosine of NUMBER."
551     (/ (+ (exp number) (exp (- number))) 2))
552    
553 ram 1.17 (defun cosh (number)
554     "Return the hyperbolic cosine of NUMBER."
555     (number-dispatch ((number number))
556     (handle-reals %cosh number)
557     ((complex)
558     (let ((x (realpart number))
559     (y (imagpart number)))
560     (complex (* (cosh x) (cos y))
561     (* (sinh x) (sin y)))))))
562    
563 wlott 1.3 (defun tanh (number)
564     "Return the hyperbolic tangent of NUMBER."
565 ram 1.17 (number-dispatch ((number number))
566     (handle-reals %tanh number)
567     ((complex)
568     (complex-tanh number))))
569 wlott 1.3
570     (defun asinh (number)
571     "Return the hyperbolic arc sine of NUMBER."
572 ram 1.17 (number-dispatch ((number number))
573     (handle-reals %asinh number)
574     ((complex)
575     (complex-asinh number))))
576 wlott 1.3
577     (defun acosh (number)
578     "Return the hyperbolic arc cosine of NUMBER."
579 ram 1.17 (number-dispatch ((number number))
580     ((rational)
581     ;; acosh is complex if number < 1
582     (if (< number 1)
583     (complex-acosh number)
584     (coerce (%acosh (coerce number 'double-float)) 'single-float)))
585     (((foreach single-float double-float))
586     (if (< number (coerce 1 '(dispatch-type number)))
587     (complex-acosh number)
588     (coerce (%acosh (coerce number 'double-float))
589     '(dispatch-type number))))
590     ((complex)
591     (complex-acosh number))))
592 wlott 1.3
593     (defun atanh (number)
594     "Return the hyperbolic arc tangent of NUMBER."
595 ram 1.17 (number-dispatch ((number number))
596     ((rational)
597     ;; atanh is complex if |number| > 1
598     (if (or (> number 1) (< number -1))
599     (complex-atanh number)
600     (coerce (%atanh (coerce number 'double-float)) 'single-float)))
601     (((foreach single-float double-float))
602     (if (or (> number (coerce 1 '(dispatch-type number)))
603     (< number (coerce -1 '(dispatch-type number))))
604     (complex-atanh number)
605     (coerce (%atanh (coerce number 'double-float))
606     '(dispatch-type number))))
607     ((complex)
608     (complex-atanh number))))
609 wlott 1.14
610 pw 1.18 ;;; HP-UX does not supply a C version of log1p, so
611     ;;; use the definition.
612 wlott 1.14
613     #+hpux
614 pw 1.18 (declaim (inline %log1p))
615 wlott 1.14 #+hpux
616 pw 1.18 (defun %log1p (number)
617     (declare (double-float number)
618     (optimize (speed 3) (safety 0)))
619     (the double-float (log (the (double-float 0d0) (+ number 1d0)))))
620 ram 1.17
621    
622     #+old-elfun
623     (progn
624     ;;; Here are the old definitions of the special functions, for
625     ;;; complex-valued arguments. Some of these functions suffer from
626     ;;; severe round-off error or unnecessary overflow.
627    
628     (proclaim '(inline mult-by-i))
629     (defun mult-by-i (number)
630     (complex (- (imagpart number))
631     (realpart number)))
632    
633     (defun complex-sqrt (x)
634     (exp (/ (log x) 2)))
635    
636     (defun complex-log (x)
637     (complex (log (abs x))
638     (phase x)))
639    
640     (defun complex-atanh (number)
641 ram 1.16 (/ (- (log (1+ number)) (log (- 1 number))) 2))
642 ram 1.17
643     (defun complex-tanh (number)
644     (/ (- (exp number) (exp (- number)))
645     (+ (exp number) (exp (- number)))))
646    
647     (defun complex-acos (number)
648     (* -2 (mult-by-i (log (+ (sqrt (/ (1+ number) 2))
649     (mult-by-i (sqrt (/ (- 1 number) 2))))))))
650    
651     (defun complex-acosh (number)
652     (* 2 (log (+ (sqrt (/ (1+ number) 2)) (sqrt (/ (1- number) 2))))))
653    
654     (defun complex-asin (number)
655     (- (mult-by-i (log (+ (mult-by-i number) (sqrt (- 1 (* number number))))))))
656    
657     (defun complex-asinh (number)
658     (log (+ number (sqrt (1+ (* number number))))))
659    
660     (defun complex-atan (y)
661     (let ((im (imagpart y))
662     (re (realpart y)))
663     (/ (- (log (complex (- 1 im) re))
664     (log (complex (+ 1 im) (- re))))
665     (complex 0 2))))
666    
667     (defun complex-tan (number)
668     (let* ((num (sin number))
669     (denom (cos number)))
670     (if (zerop denom) (error "~S undefined tangent." number)
671     (/ num denom))))
672     )
673    
674     #-old-specfun
675     (progn
676     ;;;;
677     ;;;; This is a set of routines that implement many elementary
678     ;;;; transcendental functions as specified by ANSI Common Lisp. The
679     ;;;; implementation is based on Kahan's paper.
680     ;;;;
681     ;;;; I believe I have accurately implemented the routines and are
682     ;;;; correct, but you may want to check for your self.
683     ;;;;
684     ;;;; These functions are written for CMU Lisp and take advantage of
685     ;;;; some of the features available there. It may be possible,
686     ;;;; however, to port this to other Lisps.
687     ;;;;
688     ;;;; Some functions are significantly more accurate than the original
689     ;;;; definitions in CMU Lisp. In fact, some functions in CMU Lisp
690     ;;;; give the wrong answer like (acos #c(-2.0 0.0)), where the true
691     ;;;; answer is pi + i*log(2-sqrt(3)).
692     ;;;;
693     ;;;; All of the implemented functions will take any number for an
694     ;;;; input, but the result will always be a either a complex
695     ;;;; single-float or a complex double-float.
696     ;;;;
697     ;;;; General functions
698     ;;;; complex-sqrt
699     ;;;; complex-log
700     ;;;; complex-atanh
701     ;;;; complex-tanh
702     ;;;; complex-acos
703     ;;;; complex-acosh
704     ;;;; complex-asin
705     ;;;; complex-asinh
706     ;;;; complex-atan
707     ;;;; complex-tan
708     ;;;;
709     ;;;; Utility functions:
710     ;;;; scalb logb
711     ;;;;
712     ;;;; Internal functions:
713     ;;;; square coerce-to-complex-type cssqs complex-log-scaled
714     ;;;;
715     ;;;;
716     ;;;; Please send any bug reports, comments, or improvements to Raymond
717     ;;;; Toy at toy@rtp.ericsson.se.
718     ;;;;
719     ;;;; References
720     ;;;;
721     ;;;; Kahan, W. "Branch Cuts for Complex Elementary Functions, or Much
722     ;;;; Ado About Nothing's Sign Bit" in Iserles and Powell (eds.) "The
723     ;;;; State of the Art in Numerical Analysis", pp. 165-211, Clarendon
724     ;;;; Press, 1987
725     ;;;;
726     (declaim (inline square))
727     (declaim (ftype (function (double-float) (double-float 0d0)) square))
728     (defun square (x)
729     (declare (double-float x)
730     (values (double-float 0d0)))
731     (* x x))
732    
733     ;; If you have these functions in libm, perhaps they should be used
734     ;; instead of these Lisp versions. These versions are probably good
735     ;; enough, especially since they are portable.
736    
737 dtc 1.19 (declaim (inline scalb))
738 ram 1.17 (defun scalb (x n)
739     "Compute 2^N * X without compute 2^N first (use properties of the
740     underlying floating-point format"
741     (declare (type double-float x)
742 dtc 1.19 (type double-float-exponent n))
743 ram 1.17 (scale-float x n))
744    
745     (defun logb (x)
746     "Compute an integer N such that 1 <= |2^N * x| < 2.
747     For the special cases, the following values are used:
748    
749     x logb
750     NaN NaN
751     +/- infinity +infinity
752     0 -infinity
753     "
754     (declare (type double-float x))
755 dtc 1.26 (cond ((float-nan-p x)
756 ram 1.17 x)
757     ((float-infinity-p x)
758     #.ext:double-float-positive-infinity)
759     ((zerop x)
760     ;; The answer is negative infinity, but we are supposed to
761     ;; signal divide-by-zero.
762     ;; (error 'division-by-zero :operation 'logb :operands (list x))
763     (/ -1.0d0 x)
764     )
765     (t
766     (multiple-value-bind (signif expon sign)
767     (decode-float x)
768     (declare (ignore signif sign))
769     ;; decode-float is almost right, except that the exponent
770     ;; is off by one
771     (1- expon)))))
772    
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     (number z))
783     (if (subtypep (type-of (realpart z)) 'double-float)
784     (complex x y)
785     ;; Convert anything that's not a double-float to a single-float.
786     (complex (float x 1.0)
787     (float y 1.0))))
788    
789     (defun cssqs (z)
790 dtc 1.21 ;; Compute |(x+i*y)/2^k|^2 scaled to avoid over/underflow. The
791     ;; result is r + i*k, where k is an integer.
792 ram 1.17
793     ;; Save all FP flags
794 dtc 1.21 (let ((x (float (realpart z) 1d0))
795     (y (float (imagpart z) 1d0))
796     (k 0)
797     (rho 0d0))
798 ram 1.17 (declare (double-float x y)
799     (type (double-float 0d0) rho)
800     (fixnum k))
801 dtc 1.21 ;; Would this be better handled using an exception handler to
802     ;; catch the overflow or underflow signal? For now, we turn all
803     ;; traps off and look at the accrued exceptions to see if any
804     ;; signal would have been raised.
805     (with-float-traps-masked (:underflow :overflow)
806     (setf rho (+ (square x) (square y)))
807 dtc 1.26 (cond ((and (or (float-nan-p rho)
808 dtc 1.21 (float-infinity-p rho))
809     (or (float-infinity-p (abs x))
810     (float-infinity-p (abs y))))
811     (setf rho #.ext:double-float-positive-infinity))
812     ((let ((threshold #.(/ least-positive-double-float
813     double-float-epsilon))
814 dtc 1.23 (traps (ldb vm::float-sticky-bits
815 dtc 1.21 (vm:floating-point-modes))))
816     ;; Overflow raised or (underflow raised and rho <
817     ;; lambda/eps)
818     (or (not (zerop (logand vm:float-overflow-trap-bit traps)))
819     (and (not (zerop (logand vm:float-underflow-trap-bit traps)))
820     (< rho threshold))))
821     (setf k (logb (max (abs x) (abs y))))
822     (setf rho (+ (square (scalb x (- k)))
823     (square (scalb y (- k))))))))
824     (values rho k)))
825 ram 1.17
826     (defun complex-sqrt (z)
827     "Principle square root of Z
828    
829     Z may be any number, but the result is always a complex."
830     (declare (number z))
831     (multiple-value-bind (rho k)
832     (cssqs z)
833     (declare (type (double-float 0d0) rho)
834     (fixnum k))
835     (let ((x (float (realpart z) 1.0d0))
836     (y (float (imagpart z) 1.0d0))
837     (eta 0d0)
838     (nu 0d0))
839     (declare (double-float x y eta nu))
840    
841 dtc 1.26 (if (not (float-nan-p x))
842 ram 1.17 (setf rho (+ (scalb (abs x) (- k)) (sqrt rho))))
843    
844     (cond ((oddp k)
845     (setf k (ash k -1)))
846     (t
847     (setf k (1- (ash k -1)))
848     (setf rho (+ rho rho))))
849    
850     (setf rho (scalb (sqrt rho) k))
851    
852     (setf eta rho)
853     (setf nu y)
854    
855     (when (/= rho 0d0)
856     (when (not (float-infinity-p (abs nu)))
857     (setf nu (/ (/ nu rho) 2d0)))
858     (when (< x 0d0)
859     (setf eta (abs nu))
860     (setf nu (float-sign y rho))))
861     (coerce-to-complex-type eta nu z))))
862    
863     (defun complex-log-scaled (z j)
864     "Compute log(2^j*z).
865    
866     This is for use with J /= 0 only when |z| is huge."
867     (declare (number z)
868     (fixnum j))
869     ;; The constants t0, t1, t2 should be evaluated to machine
870     ;; precision. In addition, Kahan says the accuracy of log1p
871     ;; influences the choices of these constants but doesn't say how to
872     ;; choose them. We'll just assume his choices matches our
873     ;; implementation of log1p.
874     (let ((t0 #.(/ 1 (sqrt 2.0d0)))
875     (t1 1.2d0)
876     (t2 3d0)
877     (ln2 #.(log 2d0))
878     (x (float (realpart z) 1.0d0))
879     (y (float (imagpart z) 1.0d0)))
880     (multiple-value-bind (rho k)
881     (cssqs z)
882     (declare (type (double-float 0d0) rho)
883     (fixnum k))
884     (let ((beta (max (abs x) (abs y)))
885     (theta (min (abs x) (abs y))))
886     (declare (type (double-float 0d0) beta theta))
887     (if (and (zerop k)
888     (< t0 beta)
889     (or (<= beta t1)
890     (< rho t2)))
891     (setf rho (/ (%log1p (+ (* (- beta 1.0d0)
892     (+ beta 1.0d0))
893     (* theta theta)))
894     2d0))
895     (setf rho (+ (/ (log rho) 2d0)
896     (* (+ k j) ln2))))
897     (setf theta (atan y x))
898     (coerce-to-complex-type rho theta z)))))
899    
900     (defun complex-log (z)
901     "Log of Z = log |Z| + i * arg Z
902    
903     Z may be any number, but the result is always a complex."
904     (declare (number z))
905     (complex-log-scaled z 0))
906    
907     ;; Let us note the following "strange" behavior. atanh 1.0d0 is
908     ;; +infinity, but the following code returns approx 176 + i*pi/4. The
909     ;; reason for the imaginary part is caused by the fact that arg i*y is
910     ;; never 0 since we have positive and negative zeroes.
911    
912     (defun complex-atanh (z)
913     "Compute atanh z = (log(1+z) - log(1-z))/2"
914     (declare (number z))
915     (let* (;; Constants
916     (theta #.(/ (sqrt most-positive-double-float) 4.0d0))
917     (rho #.(/ 4.0d0 (sqrt most-positive-double-float)))
918     (half-pi #.(/ pi 2.0d0))
919     (rp (float (realpart z) 1.0d0))
920     (beta (float-sign rp 1.0d0))
921     (x (* beta rp))
922     (y (* beta (- (float (imagpart z) 1.0d0))))
923     (eta 0.0d0)
924     (nu 0.0d0))
925     (declare (double-float theta rho half-pi rp beta y eta nu)
926     (type (double-float 0d0) x))
927     (cond ((or (> x theta)
928     (> (abs y) theta))
929     ;; To avoid overflow...
930     (setf eta (float-sign y half-pi))
931     ;; nu is real part of 1/(x + iy). This is x/(x^2+y^2),
932     ;; which can cause overflow. Arrange this computation so
933     ;; that it won't overflow.
934     (setf nu (let* ((x-bigger (> x (abs y)))
935     (r (if x-bigger (/ y x) (/ x y)))
936     (d (+ 1.0d0 (* r r))))
937     (declare (double-float r d))
938     (if x-bigger
939     (/ (/ x) d)
940     (/ (/ r y) d)))))
941     ((= x 1.0d0)
942     ;; Should this be changed so that if y is zero, eta is set
943     ;; to +infinity instead of approx 176? In any case
944     ;; tanh(176) is 1.0d0 within working precision.
945     (let ((t1 (+ 4d0 (square y)))
946     (t2 (+ (abs y) rho)))
947     (declare (type (double-float 0d0) t1 t2))
948     #+nil
949     (setf eta (log (/ (sqrt (sqrt t1)))
950     (sqrt t2)))
951     (setf eta (* 0.5d0 (log (the (double-float 0.0d0)
952     (/ (sqrt t1) t2)))))
953     (setf nu (* 0.5d0
954     (float-sign y
955     (+ half-pi (atan (* 0.5d0 t2))))))))
956     (t
957     (let ((t1 (+ (abs y) rho)))
958     (declare (double-float t1))
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    
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     (declare (double-float x y))
979     (cond ((> (abs x)
980 pw 1.18 #-(or linux hpux) #.(/ (%asinh most-positive-double-float) 4d0)
981 ram 1.17 ;; This is more accurate under linux.
982 pw 1.18 #+(or linux hpux) #.(/ (+ (%log 2.0d0)
983     (%log most-positive-double-float)) 4d0))
984 ram 1.17 (complex (float-sign x)
985     (float-sign y 0.0d0)))
986     (t
987 dtc 1.20 (let* ((tv (%tan y))
988 ram 1.17 (beta (+ 1.0d0 (* tv tv)))
989     (s (sinh x))
990     (rho (sqrt (+ 1.0d0 (* s s)))))
991     (declare (double-float tv s)
992     (type (double-float 0.0d0) beta rho))
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    
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     ;;.and these two expressions are equal if and only if arg conj z =
1023     ;;-arg z, which is clearly true for all z.
1024    
1025     (defun complex-acos (z)
1026     "Compute acos z = pi/2 - asin z
1027    
1028     Z may be any number, but the result is always a complex."
1029     (declare (number z))
1030     (let ((sqrt-1+z (complex-sqrt (+ 1 z)))
1031     (sqrt-1-z (complex-sqrt (- 1 z))))
1032 dtc 1.21 (with-float-traps-masked (:divide-by-zero)
1033     (complex (* 2 (atan (/ (realpart sqrt-1-z)
1034     (realpart sqrt-1+z))))
1035     (asinh (imagpart (* (conjugate sqrt-1+z)
1036     sqrt-1-z)))))))
1037 ram 1.17
1038     (defun complex-acosh (z)
1039     "Compute acosh z = 2 * log(sqrt((z+1)/2) + sqrt((z-1)/2))
1040    
1041     Z may be any number, but the result is always a complex."
1042     (declare (number z))
1043     (let ((sqrt-z-1 (complex-sqrt (- z 1)))
1044     (sqrt-z+1 (complex-sqrt (+ z 1))))
1045 dtc 1.21 (with-float-traps-masked (:divide-by-zero)
1046     (complex (asinh (realpart (* (conjugate sqrt-z-1)
1047     sqrt-z+1)))
1048     (* 2 (atan (/ (imagpart sqrt-z-1)
1049     (realpart sqrt-z+1))))))))
1050 ram 1.17
1051    
1052     (defun complex-asin (z)
1053     "Compute asin z = asinh(i*z)/i
1054    
1055     Z may be any number, but the result is always a complex."
1056     (declare (number z))
1057     (let ((sqrt-1-z (complex-sqrt (- 1 z)))
1058     (sqrt-1+z (complex-sqrt (+ 1 z))))
1059 dtc 1.21 (with-float-traps-masked (:divide-by-zero)
1060     (complex (atan (/ (realpart z)
1061     (realpart (* sqrt-1-z sqrt-1+z))))
1062     (asinh (imagpart (* (conjugate sqrt-1-z)
1063     sqrt-1+z)))))))
1064 ram 1.17
1065     (defun complex-asinh (z)
1066     "Compute asinh z = log(z + sqrt(1 + z*z))
1067    
1068     Z may be any number, but the result is always a complex."
1069     (declare (number z))
1070     ;; asinh z = -i * asin (i*z)
1071     (let* ((iz (complex (- (imagpart z)) (realpart z)))
1072     (result (complex-asin iz)))
1073     (complex (imagpart result)
1074     (- (realpart result)))))
1075    
1076     (defun complex-atan (z)
1077     "Compute atan z = atanh (i*z) / i
1078    
1079     Z may be any number, but the result is always a complex."
1080     (declare (number z))
1081     ;; atan z = -i * atanh (i*z)
1082     (let* ((iz (complex (- (imagpart z)) (realpart z)))
1083     (result (complex-atanh iz)))
1084     (complex (imagpart result)
1085     (- (realpart result)))))
1086    
1087     (defun complex-tan (z)
1088     "Compute tan z = -i * tanh(i * z)
1089    
1090     Z may be any number, but the result is always a complex."
1091     (declare (number z))
1092     ;; tan z = -i * tanh(i*z)
1093     (let* ((iz (complex (- (imagpart z)) (realpart z)))
1094     (result (complex-tanh iz)))
1095     (complex (imagpart result)
1096     (- (realpart result)))))
1097     )

  ViewVC Help
Powered by ViewVC 1.1.5