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

Contents of /src/code/irrat.lisp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.18 - (hide annotations)
Wed Feb 5 16:15:51 1997 UTC (17 years, 2 months ago) by pw
Branch: MAIN
Changes since 1.17: +18 -31 lines
initial post 1.3.7 merge
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 pw 1.18 "$Header: /tiger/var/lib/cvsroots/cmucl/src/code/irrat.lisp,v 1.18 1997/02/05 16:15:51 pw 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 ram 1.17 #-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 ram 1.17 #-x86 (def-math-rtn "sqrt" 1)
79 wlott 1.1 (def-math-rtn "hypot" 2)
80 pw 1.18 #-hpux
81 ram 1.17 (def-math-rtn "log1p" 1)
82    
83     #+x86 ;; These are needed for use by byte-compiled files.
84     (progn
85     (defun %sin (x)
86     (declare (type double-float x)
87     (values double-float))
88     (%sin x))
89     (defun %sin-quick (x)
90     (declare (double-float x)
91     (values double-float))
92     (%sin-quick x))
93     (defun %cos (x)
94     (declare (type double-float x)
95     (values double-float))
96     (%cos x))
97     (defun %cos-quick (x)
98     (declare (double-float x)
99     (values double-float))
100     (%cos-quick x))
101     (defun %tan (x)
102     (declare (double-float x)
103     (values double-float))
104     (%tan x))
105     (defun %tan-quick (x)
106     (declare (double-float x)
107     (values double-float))
108     (%tan-quick x))
109     (defun %atan (x)
110     (declare (double-float x)
111     (values double-float))
112     (%atan x))
113     (defun %atan2 (x y)
114     (declare (double-float x y)
115     (values double-float))
116     (%atan2 x y))
117     (defun %exp (x)
118     (declare (type double-float x)
119     (values double-float))
120     (%exp x))
121     (defun %log (x)
122     (declare (type double-float x)
123     (values double-float))
124     (%log x))
125     (defun %log10 (x)
126     (declare (type double-float x)
127     (values double-float))
128     (%log10 x))
129    
130     #+nil ;; notyet
131     (defun %pow (x y)
132     (declare (type (double-float 0d0) x)
133     (type double-float y)
134     (values (double-float 0d0)))
135     (%pow x y))
136     (defun %sqrt (x)
137     (declare (type double-float x)
138     (values double-float))
139     (%sqrt x))
140     (defun %scalbn (f ex)
141     (declare (type double-float f)
142     (type (signed-byte 32) ex)
143     (values double-float))
144     (%scalbn f ex))
145     (defun %scalb (f ex)
146     (declare (double-float f ex)
147     (values double-float))
148     (%scalb f ex))
149     (defun %logb (x)
150     (declare (double-float x)
151     (values double-float))
152     (%logb x))
153     ) ; progn
154 wlott 1.1
155    
156     ;;;; Power functions.
157    
158     (defun exp (number)
159     "Return e raised to the power NUMBER."
160     (number-dispatch ((number number))
161     (handle-reals %exp number)
162     ((complex)
163     (* (exp (realpart number))
164     (cis (imagpart number))))))
165    
166     ;;; INTEXP -- Handle the rational base, integer power case.
167    
168     (defparameter *intexp-maximum-exponent* 10000)
169    
170 ram 1.6 ;;; This function precisely calculates base raised to an integral power. It
171     ;;; separates the cases by the sign of power, for efficiency reasons, as powers
172     ;;; can be calculated more efficiently if power is a positive integer. Values
173     ;;; of power are calculated as positive integers, and inverted if negative.
174     ;;;
175 wlott 1.1 (defun intexp (base power)
176     (when (> (abs power) *intexp-maximum-exponent*)
177     (cerror "Continue with calculation."
178     "The absolute value of ~S exceeds ~S."
179     power '*intexp-maximum-exponent* base power))
180     (cond ((minusp power)
181     (/ (intexp base (- power))))
182     ((eql base 2)
183     (ash 1 power))
184     (t
185     (do ((nextn (ash power -1) (ash power -1))
186     (total (if (oddp power) base 1)
187     (if (oddp power) (* base total) total)))
188     ((zerop nextn) total)
189     (setq base (* base base))
190     (setq power nextn)))))
191    
192    
193 ram 1.6 ;;; EXPT -- Public
194     ;;;
195     ;;; If an integer power of a rational, use INTEXP above. Otherwise, do
196     ;;; floating point stuff. If both args are real, we try %POW right off,
197     ;;; assuming it will return 0 if the result may be complex. If so, we call
198     ;;; COMPLEX-POW which directly computes the complex result. We also separate
199     ;;; the complex-real and real-complex cases from the general complex case.
200     ;;;
201 wlott 1.1 (defun expt (base power)
202 wlott 1.3 "Returns BASE raised to the POWER."
203 wlott 1.1 (if (zerop power)
204 ram 1.6 (1+ (* base power))
205     (labels ((real-expt (base power rtype)
206     (let* ((fbase (coerce base 'double-float))
207     (fpower (coerce power 'double-float))
208     (res (coerce (%pow fbase fpower) rtype)))
209     (if (and (zerop res) (minusp fbase))
210     (multiple-value-bind (re im)
211     (complex-pow fbase fpower)
212     (%make-complex (coerce re rtype) (coerce im rtype)))
213     res)))
214     (complex-pow (fbase fpower)
215     (let ((pow (%pow (- fbase) fpower))
216     (fpower*pi (* fpower pi)))
217     (values (* pow (%cos fpower*pi))
218     (* pow (%sin fpower*pi))))))
219     (declare (inline real-expt))
220     (number-dispatch ((base number) (power number))
221     (((foreach fixnum (or bignum ratio) (complex rational)) integer)
222     (intexp base power))
223 ram 1.8 (((foreach single-float double-float) rational)
224 ram 1.6 (real-expt base power '(dispatch-type base)))
225 ram 1.9 (((foreach fixnum (or bignum ratio) single-float)
226 ram 1.6 (foreach ratio single-float))
227     (real-expt base power 'single-float))
228     (((foreach fixnum (or bignum ratio) single-float double-float)
229     double-float)
230     (real-expt base power 'double-float))
231     ((double-float single-float)
232     (real-expt base power 'double-float))
233     (((foreach (complex rational) (complex float)) rational)
234     (* (expt (abs base) power)
235     (cis (* power (phase base)))))
236     (((foreach fixnum (or bignum ratio) single-float double-float)
237     complex)
238 ram 1.17 (exp (* power (log base))))
239     (((foreach (complex float) (complex rational)) number)
240 ram 1.6 (exp (* power (log base))))))))
241 wlott 1.1
242     (defun log (number &optional (base nil base-p))
243 wlott 1.3 "Return the logarithm of NUMBER in the base BASE, which defaults to e."
244 wlott 1.1 (if base-p
245 pw 1.18 (if (zerop base)
246     base ; ANSI spec
247     (/ (log number) (log base)))
248 wlott 1.1 (number-dispatch ((number number))
249 ram 1.17 (((foreach fixnum bignum ratio))
250 wlott 1.3 (if (minusp number)
251     (complex (log (- number)) (coerce pi 'single-float))
252     (coerce (%log (coerce number 'double-float)) 'single-float)))
253 ram 1.17 (((foreach single-float double-float))
254     ;; Is (log -0) -infinity (libm.a) or -infinity + i*pi (Kahan)?
255     ;; Since this doesn't seem to be an implementation issue
256     ;; I (pw) take the Kahan result.
257     (if (< (float-sign number)
258     (coerce 0 '(dispatch-type number)))
259     (complex (log (- number)) (coerce pi '(dispatch-type number)))
260     (coerce (%log (coerce number 'double-float))
261     '(dispatch-type number))))
262     ((complex)
263     (complex-log number)))))
264 wlott 1.1
265     (defun sqrt (number)
266     "Return the square root of NUMBER."
267     (number-dispatch ((number number))
268 ram 1.17 (((foreach fixnum bignum ratio))
269 wlott 1.3 (if (minusp number)
270 ram 1.17 (complex-sqrt number)
271 wlott 1.3 (coerce (%sqrt (coerce number 'double-float)) 'single-float)))
272 ram 1.17 (((foreach single-float double-float))
273     ;; NOTE there is a problem with (at least x86 NPX) of what result
274     ;; should be returned for (sqrt -0.0). The x86 hardware FSQRT
275     ;; instruction returns -0d0. The result is that Python will perhaps
276     ;; note the following test in generic sqrt, non-negatively constrained
277     ;; float types will be passed to FSQRT (or libm on other boxes).
278     ;; So, in the interest of consistency of compiled and interpreted
279     ;; codes, the following test is disabled for now. Maybe the float-sign
280     ;; test could be moved to the optimization codes.
281     (if (< (#+nil float-sign #-nil identity number)
282     (coerce 0 '(dispatch-type number)))
283     (complex-sqrt number)
284     (coerce (%sqrt (coerce number 'double-float))
285     '(dispatch-type number))))
286     ((complex)
287     (complex-sqrt number))))
288 wlott 1.1
289    
290     ;;;; Trigonometic and Related Functions
291    
292 wlott 1.2 (defun abs (number)
293     "Returns the absolute value of the number."
294     (number-dispatch ((number number))
295     (((foreach single-float double-float fixnum rational))
296     (abs number))
297     ((complex)
298     (let ((rx (realpart number))
299     (ix (imagpart number)))
300     (etypecase rx
301     (rational
302     (sqrt (+ (* rx rx) (* ix ix))))
303     (single-float
304     (coerce (%hypot (coerce rx 'double-float)
305     (coerce ix 'double-float))
306     'single-float))
307     (double-float
308     (%hypot rx ix)))))))
309 wlott 1.1
310     (defun phase (number)
311     "Returns the angle part of the polar representation of a complex number.
312 wlott 1.3 For complex numbers, this is (atan (imagpart number) (realpart number)).
313 wlott 1.1 For non-complex positive numbers, this is 0. For non-complex negative
314     numbers this is PI."
315 wlott 1.3 (etypecase number
316 ram 1.17 (rational
317 wlott 1.3 (if (minusp number)
318     (coerce pi 'single-float)
319     0.0f0))
320 ram 1.17 (single-float
321     (if (minusp (float-sign number))
322     (coerce pi 'single-float)
323     0.0f0))
324 wlott 1.3 (double-float
325 ram 1.17 (if (minusp (float-sign number))
326 wlott 1.3 (coerce pi 'double-float)
327     0.0d0))
328     (complex
329     (atan (imagpart number) (realpart number)))))
330 wlott 1.1
331    
332     (defun sin (number)
333 wlott 1.3 "Return the sine of NUMBER."
334 wlott 1.1 (number-dispatch ((number number))
335     (handle-reals %sin number)
336     ((complex)
337     (let ((x (realpart number))
338     (y (imagpart number)))
339 ram 1.17 (complex (* (sin x) (cosh y))
340     (* (cos x) (sinh y)))))))
341 wlott 1.1
342     (defun cos (number)
343 wlott 1.3 "Return the cosine of NUMBER."
344 wlott 1.1 (number-dispatch ((number number))
345     (handle-reals %cos number)
346     ((complex)
347     (let ((x (realpart number))
348     (y (imagpart number)))
349 ram 1.17 (complex (* (cos x) (cosh y))
350     (- (* (sin x) (sinh y))))))))
351 wlott 1.1
352     (defun tan (number)
353 wlott 1.3 "Return the tangent of NUMBER."
354 wlott 1.1 (number-dispatch ((number number))
355     (handle-reals %tan number)
356     ((complex)
357 ram 1.17 (complex-tan number))))
358 wlott 1.1
359 wlott 1.2 (defun cis (theta)
360 wlott 1.3 "Return cos(Theta) + i sin(Theta), AKA exp(i Theta)."
361 wlott 1.2 (if (complexp theta)
362     (error "Argument to CIS is complex: ~S" theta)
363     (complex (cos theta) (sin theta))))
364 wlott 1.1
365     (defun asin (number)
366 wlott 1.3 "Return the arc sine of NUMBER."
367 wlott 1.1 (number-dispatch ((number number))
368 wlott 1.3 ((rational)
369     (if (or (> number 1) (< number -1))
370     (complex-asin number)
371     (coerce (%asin (coerce number 'double-float)) 'single-float)))
372     (((foreach single-float double-float))
373     (if (or (> number (coerce 1 '(dispatch-type number)))
374     (< number (coerce -1 '(dispatch-type number))))
375     (complex-asin number)
376     (coerce (%asin (coerce number 'double-float))
377     '(dispatch-type number))))
378 wlott 1.1 ((complex)
379 wlott 1.3 (complex-asin number))))
380 wlott 1.1
381     (defun acos (number)
382 wlott 1.3 "Return the arc cosine of NUMBER."
383 wlott 1.1 (number-dispatch ((number number))
384 wlott 1.3 ((rational)
385     (if (or (> number 1) (< number -1))
386     (complex-acos number)
387     (coerce (%acos (coerce number 'double-float)) 'single-float)))
388     (((foreach single-float double-float))
389     (if (or (> number (coerce 1 '(dispatch-type number)))
390     (< number (coerce -1 '(dispatch-type number))))
391     (complex-acos number)
392     (coerce (%acos (coerce number 'double-float))
393     '(dispatch-type number))))
394 wlott 1.1 ((complex)
395 wlott 1.3 (complex-acos number))))
396 wlott 1.1
397    
398     (defun atan (y &optional (x nil xp))
399 wlott 1.3 "Return the arc tangent of Y if X is omitted or Y/X if X is supplied."
400 wlott 1.1 (if xp
401 wlott 1.12 (flet ((atan2 (y x)
402 wlott 1.13 (declare (type double-float y x)
403     (values double-float))
404     (if (zerop x)
405     (if (zerop y)
406     (if (plusp (float-sign x))
407     y
408     (float-sign y pi))
409     (float-sign y (/ pi 2)))
410     (%atan2 y x))))
411     (number-dispatch ((y number) (x number))
412     ((double-float
413     (foreach double-float single-float fixnum bignum ratio))
414     (atan2 y (coerce x 'double-float)))
415     (((foreach single-float fixnum bignum ratio)
416     double-float)
417     (atan2 (coerce y 'double-float) x))
418     (((foreach single-float fixnum bignum ratio)
419     (foreach single-float fixnum bignum ratio))
420     (coerce (atan2 (coerce y 'double-float) (coerce x 'double-float))
421     'single-float))))
422 wlott 1.1 (number-dispatch ((y number))
423     (handle-reals %atan y)
424     ((complex)
425 ram 1.17 (complex-atan y)))))
426 wlott 1.3
427 ram 1.17 ;; It seems that everyone has a C version of sinh, cosh, and
428     ;; tanh. Let's use these for reals because the original
429     ;; implementations based on the definitions lose big in round-off
430     ;; error. These bad definitions also mean that sin and cos for
431     ;; complex numbers can also lose big.
432 wlott 1.3
433 ram 1.17 #+nil
434 wlott 1.3 (defun sinh (number)
435     "Return the hyperbolic sine of NUMBER."
436     (/ (- (exp number) (exp (- number))) 2))
437    
438 ram 1.17 (defun sinh (number)
439     "Return the hyperbolic sine of NUMBER."
440     (number-dispatch ((number number))
441     (handle-reals %sinh number)
442     ((complex)
443     (let ((x (realpart number))
444     (y (imagpart number)))
445     (complex (* (sinh x) (cos y))
446     (* (cosh x) (sin y)))))))
447    
448     #+nil
449 wlott 1.3 (defun cosh (number)
450     "Return the hyperbolic cosine of NUMBER."
451     (/ (+ (exp number) (exp (- number))) 2))
452    
453 ram 1.17 (defun cosh (number)
454     "Return the hyperbolic cosine of NUMBER."
455     (number-dispatch ((number number))
456     (handle-reals %cosh number)
457     ((complex)
458     (let ((x (realpart number))
459     (y (imagpart number)))
460     (complex (* (cosh x) (cos y))
461     (* (sinh x) (sin y)))))))
462    
463 wlott 1.3 (defun tanh (number)
464     "Return the hyperbolic tangent of NUMBER."
465 ram 1.17 (number-dispatch ((number number))
466     (handle-reals %tanh number)
467     ((complex)
468     (complex-tanh number))))
469 wlott 1.3
470     (defun asinh (number)
471     "Return the hyperbolic arc sine of NUMBER."
472 ram 1.17 (number-dispatch ((number number))
473     (handle-reals %asinh number)
474     ((complex)
475     (complex-asinh number))))
476 wlott 1.3
477     (defun acosh (number)
478     "Return the hyperbolic arc cosine of NUMBER."
479 ram 1.17 (number-dispatch ((number number))
480     ((rational)
481     ;; acosh is complex if number < 1
482     (if (< number 1)
483     (complex-acosh number)
484     (coerce (%acosh (coerce number 'double-float)) 'single-float)))
485     (((foreach single-float double-float))
486     (if (< number (coerce 1 '(dispatch-type number)))
487     (complex-acosh number)
488     (coerce (%acosh (coerce number 'double-float))
489     '(dispatch-type number))))
490     ((complex)
491     (complex-acosh number))))
492 wlott 1.3
493     (defun atanh (number)
494     "Return the hyperbolic arc tangent of NUMBER."
495 ram 1.17 (number-dispatch ((number number))
496     ((rational)
497     ;; atanh is complex if |number| > 1
498     (if (or (> number 1) (< number -1))
499     (complex-atanh number)
500     (coerce (%atanh (coerce number 'double-float)) 'single-float)))
501     (((foreach single-float double-float))
502     (if (or (> number (coerce 1 '(dispatch-type number)))
503     (< number (coerce -1 '(dispatch-type number))))
504     (complex-atanh number)
505     (coerce (%atanh (coerce number 'double-float))
506     '(dispatch-type number))))
507     ((complex)
508     (complex-atanh number))))
509 wlott 1.14
510 pw 1.18 ;;; HP-UX does not supply a C version of log1p, so
511     ;;; use the definition.
512 wlott 1.14
513     #+hpux
514 pw 1.18 (declaim (inline %log1p))
515 wlott 1.14 #+hpux
516 pw 1.18 (defun %log1p (number)
517     (declare (double-float number)
518     (optimize (speed 3) (safety 0)))
519     (the double-float (log (the (double-float 0d0) (+ number 1d0)))))
520 ram 1.17
521    
522     #+old-elfun
523     (progn
524     ;;; Here are the old definitions of the special functions, for
525     ;;; complex-valued arguments. Some of these functions suffer from
526     ;;; severe round-off error or unnecessary overflow.
527    
528     (proclaim '(inline mult-by-i))
529     (defun mult-by-i (number)
530     (complex (- (imagpart number))
531     (realpart number)))
532    
533     (defun complex-sqrt (x)
534     (exp (/ (log x) 2)))
535    
536     (defun complex-log (x)
537     (complex (log (abs x))
538     (phase x)))
539    
540     (defun complex-atanh (number)
541 ram 1.16 (/ (- (log (1+ number)) (log (- 1 number))) 2))
542 ram 1.17
543     (defun complex-tanh (number)
544     (/ (- (exp number) (exp (- number)))
545     (+ (exp number) (exp (- number)))))
546    
547     (defun complex-acos (number)
548     (* -2 (mult-by-i (log (+ (sqrt (/ (1+ number) 2))
549     (mult-by-i (sqrt (/ (- 1 number) 2))))))))
550    
551     (defun complex-acosh (number)
552     (* 2 (log (+ (sqrt (/ (1+ number) 2)) (sqrt (/ (1- number) 2))))))
553    
554     (defun complex-asin (number)
555     (- (mult-by-i (log (+ (mult-by-i number) (sqrt (- 1 (* number number))))))))
556    
557     (defun complex-asinh (number)
558     (log (+ number (sqrt (1+ (* number number))))))
559    
560     (defun complex-atan (y)
561     (let ((im (imagpart y))
562     (re (realpart y)))
563     (/ (- (log (complex (- 1 im) re))
564     (log (complex (+ 1 im) (- re))))
565     (complex 0 2))))
566    
567     (defun complex-tan (number)
568     (let* ((num (sin number))
569     (denom (cos number)))
570     (if (zerop denom) (error "~S undefined tangent." number)
571     (/ num denom))))
572     )
573    
574     #-old-specfun
575     (progn
576     ;;;;
577     ;;;; This is a set of routines that implement many elementary
578     ;;;; transcendental functions as specified by ANSI Common Lisp. The
579     ;;;; implementation is based on Kahan's paper.
580     ;;;;
581     ;;;; I believe I have accurately implemented the routines and are
582     ;;;; correct, but you may want to check for your self.
583     ;;;;
584     ;;;; These functions are written for CMU Lisp and take advantage of
585     ;;;; some of the features available there. It may be possible,
586     ;;;; however, to port this to other Lisps.
587     ;;;;
588     ;;;; Some functions are significantly more accurate than the original
589     ;;;; definitions in CMU Lisp. In fact, some functions in CMU Lisp
590     ;;;; give the wrong answer like (acos #c(-2.0 0.0)), where the true
591     ;;;; answer is pi + i*log(2-sqrt(3)).
592     ;;;;
593     ;;;; All of the implemented functions will take any number for an
594     ;;;; input, but the result will always be a either a complex
595     ;;;; single-float or a complex double-float.
596     ;;;;
597     ;;;; General functions
598     ;;;; complex-sqrt
599     ;;;; complex-log
600     ;;;; complex-atanh
601     ;;;; complex-tanh
602     ;;;; complex-acos
603     ;;;; complex-acosh
604     ;;;; complex-asin
605     ;;;; complex-asinh
606     ;;;; complex-atan
607     ;;;; complex-tan
608     ;;;;
609     ;;;; Utility functions:
610     ;;;; scalb logb
611     ;;;;
612     ;;;; Internal functions:
613     ;;;; square coerce-to-complex-type cssqs complex-log-scaled
614     ;;;; suppress-divide-by-zero
615     ;;;;
616     ;;;;
617     ;;;; Please send any bug reports, comments, or improvements to Raymond
618     ;;;; Toy at toy@rtp.ericsson.se.
619     ;;;;
620     ;;;; References
621     ;;;;
622     ;;;; Kahan, W. "Branch Cuts for Complex Elementary Functions, or Much
623     ;;;; Ado About Nothing's Sign Bit" in Iserles and Powell (eds.) "The
624     ;;;; State of the Art in Numerical Analysis", pp. 165-211, Clarendon
625     ;;;; Press, 1987
626     ;;;;
627     (declaim (inline square))
628     (declaim (ftype (function (double-float) (double-float 0d0)) square))
629     (defun square (x)
630     (declare (double-float x)
631     (values (double-float 0d0)))
632     (* x x))
633    
634     ;; If you have these functions in libm, perhaps they should be used
635     ;; instead of these Lisp versions. These versions are probably good
636     ;; enough, especially since they are portable.
637    
638     (defun scalb (x n)
639     "Compute 2^N * X without compute 2^N first (use properties of the
640     underlying floating-point format"
641     (declare (type double-float x)
642     (integer n))
643     (scale-float x n))
644    
645     (defun logb (x)
646     "Compute an integer N such that 1 <= |2^N * x| < 2.
647     For the special cases, the following values are used:
648    
649     x logb
650     NaN NaN
651     +/- infinity +infinity
652     0 -infinity
653     "
654     (declare (type double-float x))
655     (cond ((float-trapping-nan-p x)
656     x)
657     ((float-infinity-p x)
658     #.ext:double-float-positive-infinity)
659     ((zerop x)
660     ;; The answer is negative infinity, but we are supposed to
661     ;; signal divide-by-zero.
662     ;; (error 'division-by-zero :operation 'logb :operands (list x))
663     (/ -1.0d0 x)
664     )
665     (t
666     (multiple-value-bind (signif expon sign)
667     (decode-float x)
668     (declare (ignore signif sign))
669     ;; decode-float is almost right, except that the exponent
670     ;; is off by one
671     (1- expon)))))
672    
673     ;; This function is used to create a complex number of the appropriate
674     ;; type.
675    
676     (declaim (inline coerce-to-complex-type))
677     (defun coerce-to-complex-type (x y z)
678     "Create complex number with real part X and imaginary part Y such that
679     it has the same type as Z. If Z has type (complex rational), the X
680     and Y are coerced to single-float."
681     (declare (double-float x y)
682     (number z))
683     (if (subtypep (type-of (realpart z)) 'double-float)
684     (complex x y)
685     ;; Convert anything that's not a double-float to a single-float.
686     (complex (float x 1.0)
687     (float y 1.0))))
688    
689     (defun cssqs (z)
690     ;; Compute |x+i*y|^2/2^k carefully. The result is r + i*k, where k
691     ;; is an integer.
692    
693     ;; Save all FP flags
694     (let* ((fp-modes (get-floating-point-modes))
695     (x (float (realpart z) 1d0))
696     (y (float (imagpart z) 1d0))
697     (k 0)
698     (rho 0d0)
699     )
700     (declare (double-float x y)
701     (type (double-float 0d0) rho)
702     (fixnum k))
703     (unwind-protect
704     (progn
705     ;; Would this be better handled using an exception handler
706     ;; to catch the overflow or underflow signal? For now, we
707     ;; turn all traps off and look at the accrued exceptions to
708     ;; see if any signal would have been raised.
709    
710     (set-floating-point-modes :traps nil :accrued-exceptions nil)
711    
712     (setf rho (+ (square x) (square y)))
713     (cond ((and (or (float-trapping-nan-p rho)
714     (float-infinity-p rho))
715     (or (float-infinity-p (abs x))
716     (float-infinity-p (abs y))))
717     (setf rho #.ext:double-float-positive-infinity))
718     ((let* ((exceptions (getf (get-floating-point-modes)
719     :accrued-exceptions))
720     (threshold #.(/ least-positive-double-float
721     double-float-epsilon)))
722     ;; overflow raised or (underflow raised and rho <
723     ;; lambda/eps)
724     (or (member :overflow exceptions)
725     (and (member :underflow exceptions)
726     (< rho threshold))))
727     (setf k (logb (max (abs x) (abs y))))
728     (setf rho (+ (square (scalb x (- k)))
729     (square (scalb y (- k)))))))
730     (values rho k))
731     ;; restore over/underflow flags
732     (apply #'set-floating-point-modes fp-modes))))
733    
734     (defun complex-sqrt (z)
735     "Principle square root of Z
736    
737     Z may be any number, but the result is always a complex."
738     (declare (number z))
739     (multiple-value-bind (rho k)
740     (cssqs z)
741     (declare (type (double-float 0d0) rho)
742     (fixnum k))
743     (let ((x (float (realpart z) 1.0d0))
744     (y (float (imagpart z) 1.0d0))
745     (eta 0d0)
746     (nu 0d0))
747     (declare (double-float x y eta nu))
748    
749     (if (not (float-trapping-nan-p x))
750     (setf rho (+ (scalb (abs x) (- k)) (sqrt rho))))
751    
752     (cond ((oddp k)
753     (setf k (ash k -1)))
754     (t
755     (setf k (1- (ash k -1)))
756     (setf rho (+ rho rho))))
757    
758     (setf rho (scalb (sqrt rho) k))
759    
760     (setf eta rho)
761     (setf nu y)
762    
763     (when (/= rho 0d0)
764     (when (not (float-infinity-p (abs nu)))
765     (setf nu (/ (/ nu rho) 2d0)))
766     (when (< x 0d0)
767     (setf eta (abs nu))
768     (setf nu (float-sign y rho))))
769     (coerce-to-complex-type eta nu z))))
770    
771     (defun complex-log-scaled (z j)
772     "Compute log(2^j*z).
773    
774     This is for use with J /= 0 only when |z| is huge."
775     (declare (number z)
776     (fixnum j))
777     ;; The constants t0, t1, t2 should be evaluated to machine
778     ;; precision. In addition, Kahan says the accuracy of log1p
779     ;; influences the choices of these constants but doesn't say how to
780     ;; choose them. We'll just assume his choices matches our
781     ;; implementation of log1p.
782     (let ((t0 #.(/ 1 (sqrt 2.0d0)))
783     (t1 1.2d0)
784     (t2 3d0)
785     (ln2 #.(log 2d0))
786     (x (float (realpart z) 1.0d0))
787     (y (float (imagpart z) 1.0d0)))
788     (multiple-value-bind (rho k)
789     (cssqs z)
790     (declare (type (double-float 0d0) rho)
791     (fixnum k))
792     (let ((beta (max (abs x) (abs y)))
793     (theta (min (abs x) (abs y))))
794     (declare (type (double-float 0d0) beta theta))
795     (if (and (zerop k)
796     (< t0 beta)
797     (or (<= beta t1)
798     (< rho t2)))
799     (setf rho (/ (%log1p (+ (* (- beta 1.0d0)
800     (+ beta 1.0d0))
801     (* theta theta)))
802     2d0))
803     (setf rho (+ (/ (log rho) 2d0)
804     (* (+ k j) ln2))))
805     (setf theta (atan y x))
806     (coerce-to-complex-type rho theta z)))))
807    
808     (defun complex-log (z)
809     "Log of Z = log |Z| + i * arg Z
810    
811     Z may be any number, but the result is always a complex."
812     (declare (number z))
813     (complex-log-scaled z 0))
814    
815     ;; Let us note the following "strange" behavior. atanh 1.0d0 is
816     ;; +infinity, but the following code returns approx 176 + i*pi/4. The
817     ;; reason for the imaginary part is caused by the fact that arg i*y is
818     ;; never 0 since we have positive and negative zeroes.
819    
820     (defun complex-atanh (z)
821     "Compute atanh z = (log(1+z) - log(1-z))/2"
822     (declare (number z))
823     (let* (;; Constants
824     (theta #.(/ (sqrt most-positive-double-float) 4.0d0))
825     (rho #.(/ 4.0d0 (sqrt most-positive-double-float)))
826     (half-pi #.(/ pi 2.0d0))
827     (rp (float (realpart z) 1.0d0))
828     (beta (float-sign rp 1.0d0))
829     (x (* beta rp))
830     (y (* beta (- (float (imagpart z) 1.0d0))))
831     (eta 0.0d0)
832     (nu 0.0d0))
833     (declare (double-float theta rho half-pi rp beta y eta nu)
834     (type (double-float 0d0) x))
835     (cond ((or (> x theta)
836     (> (abs y) theta))
837     ;; To avoid overflow...
838     (setf eta (float-sign y half-pi))
839     ;; nu is real part of 1/(x + iy). This is x/(x^2+y^2),
840     ;; which can cause overflow. Arrange this computation so
841     ;; that it won't overflow.
842     (setf nu (let* ((x-bigger (> x (abs y)))
843     (r (if x-bigger (/ y x) (/ x y)))
844     (d (+ 1.0d0 (* r r))))
845     (declare (double-float r d))
846     (if x-bigger
847     (/ (/ x) d)
848     (/ (/ r y) d)))))
849     ((= x 1.0d0)
850     ;; Should this be changed so that if y is zero, eta is set
851     ;; to +infinity instead of approx 176? In any case
852     ;; tanh(176) is 1.0d0 within working precision.
853     (let ((t1 (+ 4d0 (square y)))
854     (t2 (+ (abs y) rho)))
855     (declare (type (double-float 0d0) t1 t2))
856     #+nil
857     (setf eta (log (/ (sqrt (sqrt t1)))
858     (sqrt t2)))
859     (setf eta (* 0.5d0 (log (the (double-float 0.0d0)
860     (/ (sqrt t1) t2)))))
861     (setf nu (* 0.5d0
862     (float-sign y
863     (+ half-pi (atan (* 0.5d0 t2))))))))
864     (t
865     (let ((t1 (+ (abs y) rho)))
866     (declare (double-float t1))
867     ;; Normal case using log1p(x) = log(1 + x)
868     (setf eta (* 0.25d0
869     (%log1p (/ (* 4.0d0 x)
870     (+ (square (- 1.0d0 x))
871     (square t1))))))
872     (setf nu (* 0.5d0
873     (atan (* 2.0d0 y)
874     (- (* (- 1.0d0 x)
875     (+ 1.0d0 x))
876     (square t1))))))))
877     (coerce-to-complex-type (* beta eta)
878     (- (* beta nu))
879     z)))
880    
881     (defun complex-tanh (z)
882     "Compute tanh z = sinh z / cosh z"
883     (declare (number z))
884     (let ((x (float (realpart z) 1.0d0))
885     (y (float (imagpart z) 1.0d0)))
886     (declare (double-float x y))
887     (cond ((> (abs x)
888 pw 1.18 #-(or linux hpux) #.(/ (%asinh most-positive-double-float) 4d0)
889 ram 1.17 ;; This is more accurate under linux.
890 pw 1.18 #+(or linux hpux) #.(/ (+ (%log 2.0d0)
891     (%log most-positive-double-float)) 4d0))
892 ram 1.17 (complex (float-sign x)
893     (float-sign y 0.0d0)))
894     (t
895     (let* ((tv (tan y))
896     (beta (+ 1.0d0 (* tv tv)))
897     (s (sinh x))
898     (rho (sqrt (+ 1.0d0 (* s s)))))
899     (declare (double-float tv s)
900     (type (double-float 0.0d0) beta rho))
901     (if (float-infinity-p (abs tv))
902     (coerce-to-complex-type (/ rho s)
903     (/ tv)
904     z)
905     (let ((den (+ 1.0d0 (* beta s s))))
906     (coerce-to-complex-type (/ (* beta rho s)
907     den)
908     (/ tv den)
909     z))))))))
910    
911     ;; Kahan says we should only compute the parts needed. Thus, the
912     ;; realpart's below should only compute the real part, not the whole
913     ;; complex expression. Doing this can be important because we may get
914     ;; spurious signals that occur in the part that we are not using.
915     ;;
916     ;; However, we take a pragmatic approach and just use the whole
917     ;; expression.
918    
919     ;; This macro suppresses any divide-by-zero errors in the body. After
920     ;; execution of the body, the floating point modes are returned as
921     ;; they were before entry.
922    
923     (defmacro suppress-divide-by-zero (&body body)
924     (let ((fp-modes (gensym)))
925     `(let ((,fp-modes (get-floating-point-modes)))
926     (unwind-protect
927     (progn
928     ;; Remove the divide-by-zero trap but leave everything
929     ;; else undisturbed.
930     (set-floating-point-modes
931     :traps (remove :divide-by-zero (getf ,fp-modes :traps)))
932     ,@body)
933     (apply #'set-floating-point-modes ,fp-modes)))))
934    
935     ;; NOTE: The formula given by Kahan is somewhat ambiguous in whether
936     ;; it's the conjugate of the square root or the square root of the
937     ;; conjugate. This needs to be checked.
938    
939     ;; I checked. It doesn't matter because (conjugate (sqrt z)) is the
940     ;; same as (sqrt (conjugate z)) for all z. This follows because
941     ;;
942     ;; (conjugate (sqrt z)) = exp(0.5*log |z|)*exp(-0.5*j*arg z).
943     ;;
944     ;; (sqrt (conjugate z)) = exp(0.5*log|z|)*exp(0.5*j*arg conj z)
945     ;;
946     ;;.and these two expressions are equal if and only if arg conj z =
947     ;;-arg z, which is clearly true for all z.
948    
949     (defun complex-acos (z)
950     "Compute acos z = pi/2 - asin z
951    
952     Z may be any number, but the result is always a complex."
953     (declare (number z))
954     (let ((sqrt-1+z (complex-sqrt (+ 1 z)))
955     (sqrt-1-z (complex-sqrt (- 1 z))))
956     (suppress-divide-by-zero
957     (complex (* 2 (atan (/ (realpart sqrt-1-z)
958     (realpart sqrt-1+z))))
959     (asinh (imagpart (* (conjugate sqrt-1+z)
960     sqrt-1-z)))))))
961    
962     (defun complex-acosh (z)
963     "Compute acosh z = 2 * log(sqrt((z+1)/2) + sqrt((z-1)/2))
964    
965     Z may be any number, but the result is always a complex."
966     (declare (number z))
967     (let ((sqrt-z-1 (complex-sqrt (- z 1)))
968     (sqrt-z+1 (complex-sqrt (+ z 1))))
969     (suppress-divide-by-zero
970     (complex (asinh (realpart (* (conjugate sqrt-z-1)
971     sqrt-z+1)))
972     (* 2 (atan (/ (imagpart sqrt-z-1)
973     (realpart sqrt-z+1))))))))
974    
975    
976     (defun complex-asin (z)
977     "Compute asin z = asinh(i*z)/i
978    
979     Z may be any number, but the result is always a complex."
980     (declare (number z))
981     (let ((sqrt-1-z (complex-sqrt (- 1 z)))
982     (sqrt-1+z (complex-sqrt (+ 1 z))))
983     (suppress-divide-by-zero
984     (complex (atan (/ (realpart z)
985     (realpart (* sqrt-1-z sqrt-1+z))))
986     (asinh (imagpart (* (conjugate sqrt-1-z)
987     sqrt-1+z)))))))
988    
989     (defun complex-asinh (z)
990     "Compute asinh z = log(z + sqrt(1 + z*z))
991    
992     Z may be any number, but the result is always a complex."
993     (declare (number z))
994     ;; asinh z = -i * asin (i*z)
995     (let* ((iz (complex (- (imagpart z)) (realpart z)))
996     (result (complex-asin iz)))
997     (complex (imagpart result)
998     (- (realpart result)))))
999    
1000     (defun complex-atan (z)
1001     "Compute atan z = atanh (i*z) / i
1002    
1003     Z may be any number, but the result is always a complex."
1004     (declare (number z))
1005     ;; atan z = -i * atanh (i*z)
1006     (let* ((iz (complex (- (imagpart z)) (realpart z)))
1007     (result (complex-atanh iz)))
1008     (complex (imagpart result)
1009     (- (realpart result)))))
1010    
1011     (defun complex-tan (z)
1012     "Compute tan z = -i * tanh(i * z)
1013    
1014     Z may be any number, but the result is always a complex."
1015     (declare (number z))
1016     ;; tan z = -i * tanh(i*z)
1017     (let* ((iz (complex (- (imagpart z)) (realpart z)))
1018     (result (complex-tanh iz)))
1019     (complex (imagpart result)
1020     (- (realpart result)))))
1021     )

  ViewVC Help
Powered by ViewVC 1.1.5