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

Contents of /src/code/irrat.lisp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.19 - (hide annotations)
Mon Feb 10 23:24:21 1997 UTC (17 years, 2 months ago) by dtc
Branch: MAIN
Branch point for: RELENG_18
Changes since 1.18: +17 -17 lines
Inline scalb and better declare its 'n' argument so the compiler can
generate inline code for this operation. Few other little cleanups.
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.19 "$Header: /tiger/var/lib/cvsroots/cmucl/src/code/irrat.lisp,v 1.19 1997/02/10 23:24:21 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 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 dtc 1.19 (declare (double-float x)
87 ram 1.17 (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 dtc 1.19 (declare (double-float x)
95 ram 1.17 (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 dtc 1.19 (declare (double-float x)
119 ram 1.17 (values double-float))
120     (%exp x))
121     (defun %log (x)
122 dtc 1.19 (declare (double-float x)
123 ram 1.17 (values double-float))
124     (%log x))
125     (defun %log10 (x)
126 dtc 1.19 (declare (double-float x)
127 ram 1.17 (values double-float))
128     (%log10 x))
129     #+nil ;; notyet
130     (defun %pow (x y)
131     (declare (type (double-float 0d0) x)
132 dtc 1.19 (double-float y)
133 ram 1.17 (values (double-float 0d0)))
134     (%pow x y))
135     (defun %sqrt (x)
136 dtc 1.19 (declare (double-float x)
137 ram 1.17 (values double-float))
138     (%sqrt x))
139     (defun %scalbn (f ex)
140 dtc 1.19 (declare (double-float f)
141 ram 1.17 (type (signed-byte 32) ex)
142     (values double-float))
143     (%scalbn f ex))
144     (defun %scalb (f ex)
145     (declare (double-float f ex)
146     (values double-float))
147     (%scalb f ex))
148     (defun %logb (x)
149 dtc 1.19 (declare (double-float x)
150     (values double-float))
151     (%logb x))
152 ram 1.17 ) ; progn
153 wlott 1.1
154    
155     ;;;; Power functions.
156    
157     (defun exp (number)
158     "Return e raised to the power NUMBER."
159     (number-dispatch ((number number))
160     (handle-reals %exp number)
161     ((complex)
162     (* (exp (realpart number))
163     (cis (imagpart number))))))
164    
165     ;;; INTEXP -- Handle the rational base, integer power case.
166    
167     (defparameter *intexp-maximum-exponent* 10000)
168    
169 ram 1.6 ;;; This function precisely calculates base raised to an integral power. It
170     ;;; separates the cases by the sign of power, for efficiency reasons, as powers
171     ;;; can be calculated more efficiently if power is a positive integer. Values
172     ;;; of power are calculated as positive integers, and inverted if negative.
173     ;;;
174 wlott 1.1 (defun intexp (base power)
175     (when (> (abs power) *intexp-maximum-exponent*)
176     (cerror "Continue with calculation."
177     "The absolute value of ~S exceeds ~S."
178     power '*intexp-maximum-exponent* base power))
179     (cond ((minusp power)
180     (/ (intexp base (- power))))
181     ((eql base 2)
182     (ash 1 power))
183     (t
184     (do ((nextn (ash power -1) (ash power -1))
185     (total (if (oddp power) base 1)
186     (if (oddp power) (* base total) total)))
187     ((zerop nextn) total)
188     (setq base (* base base))
189     (setq power nextn)))))
190    
191    
192 ram 1.6 ;;; EXPT -- Public
193     ;;;
194     ;;; If an integer power of a rational, use INTEXP above. Otherwise, do
195     ;;; floating point stuff. If both args are real, we try %POW right off,
196     ;;; assuming it will return 0 if the result may be complex. If so, we call
197     ;;; COMPLEX-POW which directly computes the complex result. We also separate
198     ;;; the complex-real and real-complex cases from the general complex case.
199     ;;;
200 wlott 1.1 (defun expt (base power)
201 wlott 1.3 "Returns BASE raised to the POWER."
202 wlott 1.1 (if (zerop power)
203 ram 1.6 (1+ (* base power))
204     (labels ((real-expt (base power rtype)
205     (let* ((fbase (coerce base 'double-float))
206     (fpower (coerce power 'double-float))
207     (res (coerce (%pow fbase fpower) rtype)))
208     (if (and (zerop res) (minusp fbase))
209     (multiple-value-bind (re im)
210     (complex-pow fbase fpower)
211     (%make-complex (coerce re rtype) (coerce im rtype)))
212     res)))
213     (complex-pow (fbase fpower)
214     (let ((pow (%pow (- fbase) fpower))
215     (fpower*pi (* fpower pi)))
216     (values (* pow (%cos fpower*pi))
217     (* pow (%sin fpower*pi))))))
218     (declare (inline real-expt))
219     (number-dispatch ((base number) (power number))
220     (((foreach fixnum (or bignum ratio) (complex rational)) integer)
221     (intexp base power))
222 ram 1.8 (((foreach single-float double-float) rational)
223 ram 1.6 (real-expt base power '(dispatch-type base)))
224 ram 1.9 (((foreach fixnum (or bignum ratio) single-float)
225 ram 1.6 (foreach ratio single-float))
226     (real-expt base power 'single-float))
227     (((foreach fixnum (or bignum ratio) single-float double-float)
228     double-float)
229     (real-expt base power 'double-float))
230     ((double-float single-float)
231     (real-expt base power 'double-float))
232     (((foreach (complex rational) (complex float)) rational)
233     (* (expt (abs base) power)
234     (cis (* power (phase base)))))
235     (((foreach fixnum (or bignum ratio) single-float double-float)
236     complex)
237 ram 1.17 (exp (* power (log base))))
238     (((foreach (complex float) (complex rational)) number)
239 ram 1.6 (exp (* power (log base))))))))
240 wlott 1.1
241     (defun log (number &optional (base nil base-p))
242 wlott 1.3 "Return the logarithm of NUMBER in the base BASE, which defaults to e."
243 wlott 1.1 (if base-p
244 pw 1.18 (if (zerop base)
245     base ; ANSI spec
246     (/ (log number) (log base)))
247 wlott 1.1 (number-dispatch ((number number))
248 ram 1.17 (((foreach fixnum bignum ratio))
249 wlott 1.3 (if (minusp number)
250     (complex (log (- number)) (coerce pi 'single-float))
251     (coerce (%log (coerce number 'double-float)) 'single-float)))
252 ram 1.17 (((foreach single-float double-float))
253     ;; Is (log -0) -infinity (libm.a) or -infinity + i*pi (Kahan)?
254     ;; Since this doesn't seem to be an implementation issue
255     ;; I (pw) take the Kahan result.
256     (if (< (float-sign number)
257     (coerce 0 '(dispatch-type number)))
258     (complex (log (- number)) (coerce pi '(dispatch-type number)))
259     (coerce (%log (coerce number 'double-float))
260     '(dispatch-type number))))
261     ((complex)
262     (complex-log number)))))
263 wlott 1.1
264     (defun sqrt (number)
265     "Return the square root of NUMBER."
266     (number-dispatch ((number number))
267 ram 1.17 (((foreach fixnum bignum ratio))
268 wlott 1.3 (if (minusp number)
269 ram 1.17 (complex-sqrt number)
270 wlott 1.3 (coerce (%sqrt (coerce number 'double-float)) 'single-float)))
271 ram 1.17 (((foreach single-float double-float))
272     ;; NOTE there is a problem with (at least x86 NPX) of what result
273     ;; should be returned for (sqrt -0.0). The x86 hardware FSQRT
274     ;; instruction returns -0d0. The result is that Python will perhaps
275     ;; note the following test in generic sqrt, non-negatively constrained
276     ;; float types will be passed to FSQRT (or libm on other boxes).
277     ;; So, in the interest of consistency of compiled and interpreted
278     ;; codes, the following test is disabled for now. Maybe the float-sign
279     ;; test could be moved to the optimization codes.
280     (if (< (#+nil float-sign #-nil identity number)
281     (coerce 0 '(dispatch-type number)))
282     (complex-sqrt number)
283     (coerce (%sqrt (coerce number 'double-float))
284     '(dispatch-type number))))
285     ((complex)
286     (complex-sqrt number))))
287 wlott 1.1
288    
289     ;;;; Trigonometic and Related Functions
290    
291 wlott 1.2 (defun abs (number)
292     "Returns the absolute value of the number."
293     (number-dispatch ((number number))
294     (((foreach single-float double-float fixnum rational))
295     (abs number))
296     ((complex)
297     (let ((rx (realpart number))
298     (ix (imagpart number)))
299     (etypecase rx
300     (rational
301     (sqrt (+ (* rx rx) (* ix ix))))
302     (single-float
303     (coerce (%hypot (coerce rx 'double-float)
304     (coerce ix 'double-float))
305     'single-float))
306     (double-float
307     (%hypot rx ix)))))))
308 wlott 1.1
309     (defun phase (number)
310     "Returns the angle part of the polar representation of a complex number.
311 wlott 1.3 For complex numbers, this is (atan (imagpart number) (realpart number)).
312 wlott 1.1 For non-complex positive numbers, this is 0. For non-complex negative
313     numbers this is PI."
314 wlott 1.3 (etypecase number
315 ram 1.17 (rational
316 wlott 1.3 (if (minusp number)
317     (coerce pi 'single-float)
318     0.0f0))
319 ram 1.17 (single-float
320     (if (minusp (float-sign number))
321     (coerce pi 'single-float)
322     0.0f0))
323 wlott 1.3 (double-float
324 ram 1.17 (if (minusp (float-sign number))
325 wlott 1.3 (coerce pi 'double-float)
326     0.0d0))
327     (complex
328     (atan (imagpart number) (realpart number)))))
329 wlott 1.1
330    
331     (defun sin (number)
332 wlott 1.3 "Return the sine of NUMBER."
333 wlott 1.1 (number-dispatch ((number number))
334     (handle-reals %sin number)
335     ((complex)
336     (let ((x (realpart number))
337     (y (imagpart number)))
338 ram 1.17 (complex (* (sin x) (cosh y))
339     (* (cos x) (sinh y)))))))
340 wlott 1.1
341     (defun cos (number)
342 wlott 1.3 "Return the cosine of NUMBER."
343 wlott 1.1 (number-dispatch ((number number))
344     (handle-reals %cos number)
345     ((complex)
346     (let ((x (realpart number))
347     (y (imagpart number)))
348 ram 1.17 (complex (* (cos x) (cosh y))
349     (- (* (sin x) (sinh y))))))))
350 wlott 1.1
351     (defun tan (number)
352 wlott 1.3 "Return the tangent of NUMBER."
353 wlott 1.1 (number-dispatch ((number number))
354     (handle-reals %tan number)
355     ((complex)
356 ram 1.17 (complex-tan number))))
357 wlott 1.1
358 wlott 1.2 (defun cis (theta)
359 wlott 1.3 "Return cos(Theta) + i sin(Theta), AKA exp(i Theta)."
360 wlott 1.2 (if (complexp theta)
361     (error "Argument to CIS is complex: ~S" theta)
362     (complex (cos theta) (sin theta))))
363 wlott 1.1
364     (defun asin (number)
365 wlott 1.3 "Return the arc sine of NUMBER."
366 wlott 1.1 (number-dispatch ((number number))
367 wlott 1.3 ((rational)
368     (if (or (> number 1) (< number -1))
369     (complex-asin number)
370     (coerce (%asin (coerce number 'double-float)) 'single-float)))
371     (((foreach single-float double-float))
372     (if (or (> number (coerce 1 '(dispatch-type number)))
373     (< number (coerce -1 '(dispatch-type number))))
374     (complex-asin number)
375     (coerce (%asin (coerce number 'double-float))
376     '(dispatch-type number))))
377 wlott 1.1 ((complex)
378 wlott 1.3 (complex-asin number))))
379 wlott 1.1
380     (defun acos (number)
381 wlott 1.3 "Return the arc cosine of NUMBER."
382 wlott 1.1 (number-dispatch ((number number))
383 wlott 1.3 ((rational)
384     (if (or (> number 1) (< number -1))
385     (complex-acos number)
386     (coerce (%acos (coerce number 'double-float)) 'single-float)))
387     (((foreach single-float double-float))
388     (if (or (> number (coerce 1 '(dispatch-type number)))
389     (< number (coerce -1 '(dispatch-type number))))
390     (complex-acos number)
391     (coerce (%acos (coerce number 'double-float))
392     '(dispatch-type number))))
393 wlott 1.1 ((complex)
394 wlott 1.3 (complex-acos number))))
395 wlott 1.1
396    
397     (defun atan (y &optional (x nil xp))
398 wlott 1.3 "Return the arc tangent of Y if X is omitted or Y/X if X is supplied."
399 wlott 1.1 (if xp
400 wlott 1.12 (flet ((atan2 (y x)
401 wlott 1.13 (declare (type double-float y x)
402     (values double-float))
403     (if (zerop x)
404     (if (zerop y)
405     (if (plusp (float-sign x))
406     y
407     (float-sign y pi))
408     (float-sign y (/ pi 2)))
409     (%atan2 y x))))
410     (number-dispatch ((y number) (x number))
411     ((double-float
412     (foreach double-float single-float fixnum bignum ratio))
413     (atan2 y (coerce x 'double-float)))
414     (((foreach single-float fixnum bignum ratio)
415     double-float)
416     (atan2 (coerce y 'double-float) x))
417     (((foreach single-float fixnum bignum ratio)
418     (foreach single-float fixnum bignum ratio))
419     (coerce (atan2 (coerce y 'double-float) (coerce x 'double-float))
420     'single-float))))
421 wlott 1.1 (number-dispatch ((y number))
422     (handle-reals %atan y)
423     ((complex)
424 ram 1.17 (complex-atan y)))))
425 wlott 1.3
426 ram 1.17 ;; It seems that everyone has a C version of sinh, cosh, and
427     ;; tanh. Let's use these for reals because the original
428     ;; implementations based on the definitions lose big in round-off
429     ;; error. These bad definitions also mean that sin and cos for
430     ;; complex numbers can also lose big.
431 wlott 1.3
432 ram 1.17 #+nil
433 wlott 1.3 (defun sinh (number)
434     "Return the hyperbolic sine of NUMBER."
435     (/ (- (exp number) (exp (- number))) 2))
436    
437 ram 1.17 (defun sinh (number)
438     "Return the hyperbolic sine of NUMBER."
439     (number-dispatch ((number number))
440     (handle-reals %sinh number)
441     ((complex)
442     (let ((x (realpart number))
443     (y (imagpart number)))
444     (complex (* (sinh x) (cos y))
445     (* (cosh x) (sin y)))))))
446    
447     #+nil
448 wlott 1.3 (defun cosh (number)
449     "Return the hyperbolic cosine of NUMBER."
450     (/ (+ (exp number) (exp (- number))) 2))
451    
452 ram 1.17 (defun cosh (number)
453     "Return the hyperbolic cosine of NUMBER."
454     (number-dispatch ((number number))
455     (handle-reals %cosh number)
456     ((complex)
457     (let ((x (realpart number))
458     (y (imagpart number)))
459     (complex (* (cosh x) (cos y))
460     (* (sinh x) (sin y)))))))
461    
462 wlott 1.3 (defun tanh (number)
463     "Return the hyperbolic tangent of NUMBER."
464 ram 1.17 (number-dispatch ((number number))
465     (handle-reals %tanh number)
466     ((complex)
467     (complex-tanh number))))
468 wlott 1.3
469     (defun asinh (number)
470     "Return the hyperbolic arc sine of NUMBER."
471 ram 1.17 (number-dispatch ((number number))
472     (handle-reals %asinh number)
473     ((complex)
474     (complex-asinh number))))
475 wlott 1.3
476     (defun acosh (number)
477     "Return the hyperbolic arc cosine of NUMBER."
478 ram 1.17 (number-dispatch ((number number))
479     ((rational)
480     ;; acosh is complex if number < 1
481     (if (< number 1)
482     (complex-acosh number)
483     (coerce (%acosh (coerce number 'double-float)) 'single-float)))
484     (((foreach single-float double-float))
485     (if (< number (coerce 1 '(dispatch-type number)))
486     (complex-acosh number)
487     (coerce (%acosh (coerce number 'double-float))
488     '(dispatch-type number))))
489     ((complex)
490     (complex-acosh number))))
491 wlott 1.3
492     (defun atanh (number)
493     "Return the hyperbolic arc tangent of NUMBER."
494 ram 1.17 (number-dispatch ((number number))
495     ((rational)
496     ;; atanh is complex if |number| > 1
497     (if (or (> number 1) (< number -1))
498     (complex-atanh number)
499     (coerce (%atanh (coerce number 'double-float)) 'single-float)))
500     (((foreach single-float double-float))
501     (if (or (> number (coerce 1 '(dispatch-type number)))
502     (< number (coerce -1 '(dispatch-type number))))
503     (complex-atanh number)
504     (coerce (%atanh (coerce number 'double-float))
505     '(dispatch-type number))))
506     ((complex)
507     (complex-atanh number))))
508 wlott 1.14
509 pw 1.18 ;;; HP-UX does not supply a C version of log1p, so
510     ;;; use the definition.
511 wlott 1.14
512     #+hpux
513 pw 1.18 (declaim (inline %log1p))
514 wlott 1.14 #+hpux
515 pw 1.18 (defun %log1p (number)
516     (declare (double-float number)
517     (optimize (speed 3) (safety 0)))
518     (the double-float (log (the (double-float 0d0) (+ number 1d0)))))
519 ram 1.17
520    
521     #+old-elfun
522     (progn
523     ;;; Here are the old definitions of the special functions, for
524     ;;; complex-valued arguments. Some of these functions suffer from
525     ;;; severe round-off error or unnecessary overflow.
526    
527     (proclaim '(inline mult-by-i))
528     (defun mult-by-i (number)
529     (complex (- (imagpart number))
530     (realpart number)))
531    
532     (defun complex-sqrt (x)
533     (exp (/ (log x) 2)))
534    
535     (defun complex-log (x)
536     (complex (log (abs x))
537     (phase x)))
538    
539     (defun complex-atanh (number)
540 ram 1.16 (/ (- (log (1+ number)) (log (- 1 number))) 2))
541 ram 1.17
542     (defun complex-tanh (number)
543     (/ (- (exp number) (exp (- number)))
544     (+ (exp number) (exp (- number)))))
545    
546     (defun complex-acos (number)
547     (* -2 (mult-by-i (log (+ (sqrt (/ (1+ number) 2))
548     (mult-by-i (sqrt (/ (- 1 number) 2))))))))
549    
550     (defun complex-acosh (number)
551     (* 2 (log (+ (sqrt (/ (1+ number) 2)) (sqrt (/ (1- number) 2))))))
552    
553     (defun complex-asin (number)
554     (- (mult-by-i (log (+ (mult-by-i number) (sqrt (- 1 (* number number))))))))
555    
556     (defun complex-asinh (number)
557     (log (+ number (sqrt (1+ (* number number))))))
558    
559     (defun complex-atan (y)
560     (let ((im (imagpart y))
561     (re (realpart y)))
562     (/ (- (log (complex (- 1 im) re))
563     (log (complex (+ 1 im) (- re))))
564     (complex 0 2))))
565    
566     (defun complex-tan (number)
567     (let* ((num (sin number))
568     (denom (cos number)))
569     (if (zerop denom) (error "~S undefined tangent." number)
570     (/ num denom))))
571     )
572    
573     #-old-specfun
574     (progn
575     ;;;;
576     ;;;; This is a set of routines that implement many elementary
577     ;;;; transcendental functions as specified by ANSI Common Lisp. The
578     ;;;; implementation is based on Kahan's paper.
579     ;;;;
580     ;;;; I believe I have accurately implemented the routines and are
581     ;;;; correct, but you may want to check for your self.
582     ;;;;
583     ;;;; These functions are written for CMU Lisp and take advantage of
584     ;;;; some of the features available there. It may be possible,
585     ;;;; however, to port this to other Lisps.
586     ;;;;
587     ;;;; Some functions are significantly more accurate than the original
588     ;;;; definitions in CMU Lisp. In fact, some functions in CMU Lisp
589     ;;;; give the wrong answer like (acos #c(-2.0 0.0)), where the true
590     ;;;; answer is pi + i*log(2-sqrt(3)).
591     ;;;;
592     ;;;; All of the implemented functions will take any number for an
593     ;;;; input, but the result will always be a either a complex
594     ;;;; single-float or a complex double-float.
595     ;;;;
596     ;;;; General functions
597     ;;;; complex-sqrt
598     ;;;; complex-log
599     ;;;; complex-atanh
600     ;;;; complex-tanh
601     ;;;; complex-acos
602     ;;;; complex-acosh
603     ;;;; complex-asin
604     ;;;; complex-asinh
605     ;;;; complex-atan
606     ;;;; complex-tan
607     ;;;;
608     ;;;; Utility functions:
609     ;;;; scalb logb
610     ;;;;
611     ;;;; Internal functions:
612     ;;;; square coerce-to-complex-type cssqs complex-log-scaled
613     ;;;; suppress-divide-by-zero
614     ;;;;
615     ;;;;
616     ;;;; Please send any bug reports, comments, or improvements to Raymond
617     ;;;; Toy at toy@rtp.ericsson.se.
618     ;;;;
619     ;;;; References
620     ;;;;
621     ;;;; Kahan, W. "Branch Cuts for Complex Elementary Functions, or Much
622     ;;;; Ado About Nothing's Sign Bit" in Iserles and Powell (eds.) "The
623     ;;;; State of the Art in Numerical Analysis", pp. 165-211, Clarendon
624     ;;;; Press, 1987
625     ;;;;
626     (declaim (inline square))
627     (declaim (ftype (function (double-float) (double-float 0d0)) square))
628     (defun square (x)
629     (declare (double-float x)
630     (values (double-float 0d0)))
631     (* x x))
632    
633     ;; If you have these functions in libm, perhaps they should be used
634     ;; instead of these Lisp versions. These versions are probably good
635     ;; enough, especially since they are portable.
636    
637 dtc 1.19 (declaim (inline scalb))
638 ram 1.17 (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 dtc 1.19 (type double-float-exponent n))
643 ram 1.17 (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