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

Contents of /src/code/irrat.lisp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.60.2.2 - (hide annotations)
Tue Feb 9 20:23:02 2010 UTC (4 years, 2 months ago) by rtoy
Branch: intl-branch
Changes since 1.60.2.1: +39 -39 lines
Mark translatable strings; update cmucl.pot and ko/cmucl.po
accordingly.
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 rtoy 1.60.2.2 "$Header: /tiger/var/lib/cvsroots/cmucl/src/code/irrat.lisp,v 1.60.2.2 2010/02/09 20:23:02 rtoy 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 rtoy 1.60.2.1 (intl:textdomain "cmucl")
20 wlott 1.1
21    
22     ;;;; Random constants, utility functions, and macros.
23    
24     (defconstant pi 3.14159265358979323846264338327950288419716939937511L0)
25 wlott 1.2 ;(defconstant e 2.71828182845904523536028747135266249775724709369996L0)
26 wlott 1.1
27 ram 1.5 ;;; Make these INLINE, since the call to C is at least as compact as a Lisp
28     ;;; call, and saves number consing to boot.
29     ;;;
30 wlott 1.1 (defmacro def-math-rtn (name num-args)
31     (let ((function (intern (concatenate 'simple-string
32     "%"
33     (string-upcase name)))))
34 ram 1.4 `(progn
35 pw 1.31 (declaim (inline ,function))
36 ram 1.4 (export ',function)
37 wlott 1.10 (alien:def-alien-routine (,name ,function) double-float
38 ram 1.4 ,@(let ((results nil))
39     (dotimes (i num-args (nreverse results))
40     (push (list (intern (format nil "ARG-~D" i))
41     'double-float)
42     results)))))))
43 wlott 1.1
44     (eval-when (compile load eval)
45    
46     (defun handle-reals (function var)
47     `((((foreach fixnum single-float bignum ratio))
48     (coerce (,function (coerce ,var 'double-float)) 'single-float))
49     ((double-float)
50 rtoy 1.46 (,function ,var))
51     #+double-double
52     ((double-double-float)
53     (,(symbolicate "DD-" function) ,var))))
54 wlott 1.1
55     ); eval-when (compile load eval)
56    
57    
58     ;;;; Stubs for the Unix math library.
59    
60     ;;; Please refer to the Unix man pages for details about these routines.
61    
62     ;;; Trigonometric.
63 rtoy 1.58 #-(and x86 (not sse2))
64     (progn
65     ;; For x86 (without sse2), we can use x87 instructions to implement
66     ;; these. With sse2, we don't currently support that, so these
67     ;; should be disabled.
68     (def-math-rtn "sin" 1)
69     (def-math-rtn "cos" 1)
70     (def-math-rtn "tan" 1)
71     (def-math-rtn "atan" 1)
72     (def-math-rtn "atan2" 2))
73 wlott 1.1 (def-math-rtn "asin" 1)
74     (def-math-rtn "acos" 1)
75     (def-math-rtn "sinh" 1)
76     (def-math-rtn "cosh" 1)
77     (def-math-rtn "tanh" 1)
78 pw 1.18 (def-math-rtn "asinh" 1)
79     (def-math-rtn "acosh" 1)
80     (def-math-rtn "atanh" 1)
81 wlott 1.1
82     ;;; Exponential and Logarithmic.
83 rtoy 1.58 #-(and x86 (not sse2))
84     (progn
85     (def-math-rtn "exp" 1)
86     (def-math-rtn "log" 1)
87     (def-math-rtn "log10" 1))
88    
89 wlott 1.1 (def-math-rtn "pow" 2)
90 rtoy 1.58 #-(or x86 sparc-v7 sparc-v8 sparc-v9)
91     (def-math-rtn "sqrt" 1)
92 wlott 1.1 (def-math-rtn "hypot" 2)
93 ram 1.17
94 rtoy 1.58 ;; Don't want log1p to use the x87 instruction.
95     #-(or hpux (and x86 (not sse2)))
96     (def-math-rtn "log1p" 1)
97    
98     ;; These are needed for use by byte-compiled files. But don't use
99     ;; these with sse2 since we don't support using the x87 instructions
100     ;; here.
101     #+(and x86 (not sse2))
102 ram 1.17 (progn
103 rtoy 1.47 #+nil
104 ram 1.17 (defun %sin (x)
105 dtc 1.19 (declare (double-float x)
106 ram 1.17 (values double-float))
107     (%sin x))
108     (defun %sin-quick (x)
109     (declare (double-float x)
110     (values double-float))
111     (%sin-quick x))
112 rtoy 1.47 #+nil
113 ram 1.17 (defun %cos (x)
114 dtc 1.19 (declare (double-float x)
115 ram 1.17 (values double-float))
116     (%cos x))
117     (defun %cos-quick (x)
118     (declare (double-float x)
119     (values double-float))
120     (%cos-quick x))
121 rtoy 1.47 #+nil
122 ram 1.17 (defun %tan (x)
123     (declare (double-float x)
124     (values double-float))
125     (%tan x))
126     (defun %tan-quick (x)
127     (declare (double-float x)
128     (values double-float))
129     (%tan-quick x))
130     (defun %atan (x)
131     (declare (double-float x)
132     (values double-float))
133     (%atan x))
134     (defun %atan2 (x y)
135     (declare (double-float x y)
136     (values double-float))
137     (%atan2 x y))
138     (defun %exp (x)
139 dtc 1.19 (declare (double-float x)
140 ram 1.17 (values double-float))
141     (%exp x))
142     (defun %log (x)
143 dtc 1.19 (declare (double-float x)
144 ram 1.17 (values double-float))
145     (%log x))
146     (defun %log10 (x)
147 dtc 1.19 (declare (double-float x)
148 ram 1.17 (values double-float))
149     (%log10 x))
150     #+nil ;; notyet
151     (defun %pow (x y)
152     (declare (type (double-float 0d0) x)
153 dtc 1.19 (double-float y)
154 ram 1.17 (values (double-float 0d0)))
155     (%pow x y))
156     (defun %sqrt (x)
157 dtc 1.19 (declare (double-float x)
158 ram 1.17 (values double-float))
159     (%sqrt x))
160     (defun %scalbn (f ex)
161 dtc 1.19 (declare (double-float f)
162 ram 1.17 (type (signed-byte 32) ex)
163     (values double-float))
164     (%scalbn f ex))
165     (defun %scalb (f ex)
166     (declare (double-float f ex)
167     (values double-float))
168     (%scalb f ex))
169     (defun %logb (x)
170 dtc 1.19 (declare (double-float x)
171     (values double-float))
172     (%logb x))
173 dtc 1.27 (defun %log1p (x)
174     (declare (double-float x)
175     (values double-float))
176     (%log1p x))
177 ram 1.17 ) ; progn
178 wlott 1.1
179 rtoy 1.44
180     ;; As above for x86. It also seems to be needed to handle
181     ;; constant-folding in the compiler.
182     #+sparc
183     (progn
184     (defun %sqrt (x)
185     (declare (double-float x)
186     (values double-float))
187     (%sqrt x))
188     )
189    
190 rtoy 1.49 ;;; The standard libm routines for sin, cos, and tan on x86 (Linux)
191     ;;; and ppc are not very accurate for large arguments when compared to
192     ;;; sparc (and maxima). This is basically caused by the fact that
193     ;;; those libraries do not do an accurate argument reduction. The
194     ;;; following functions use some routines Sun's free fdlibm library to
195     ;;; do accurate reduction. Then we call the standard C functions (or
196     ;;; vops for x86) on the reduced argument. This produces much more
197     ;;; accurate values.
198    
199 rtoy 1.48 #+(or ppc x86)
200 rtoy 1.47 (progn
201 rtoy 1.54 (declaim (inline %%ieee754-rem-pi/2))
202 rtoy 1.49 ;; Basic argument reduction routine. It returns two values: n and y
203     ;; such that (n + 8*k)*pi/2+y = x where |y|<pi/4 and n indicates in
204     ;; which octant the arg lies. Y is actually computed in two parts,
205     ;; y[0] and y[1] such that the sum is y, for accuracy.
206    
207 rtoy 1.54 (alien:def-alien-routine ("__ieee754_rem_pio2" %%ieee754-rem-pi/2) c-call:int
208 rtoy 1.47 (x double-float)
209     (y (* double-float)))
210 rtoy 1.54
211     ;; Same as above, but instead of needing to pass an array in, the
212     ;; output array is broken up into two output values instead. This is
213     ;; easier for the user, and we don't have to wrap calls with
214     ;; without-gcing.
215     (declaim (inline %ieee754-rem-pi/2))
216     (alien:def-alien-routine ("ieee754_rem_pio2" %ieee754-rem-pi/2) c-call:int
217     (x double-float)
218     (y0 double-float :out)
219     (y1 double-float :out))
220    
221 rtoy 1.48 )
222 rtoy 1.47
223 rtoy 1.58 #+(or ppc sse2)
224 rtoy 1.48 (progn
225     (declaim (inline %%sin %%cos %%tan))
226     (macrolet ((frob (alien-name lisp-name)
227     `(alien:def-alien-routine (,alien-name ,lisp-name) double-float
228     (x double-float))))
229     (frob "sin" %%sin)
230     (frob "cos" %%cos)
231     (frob "tan" %%tan))
232 rtoy 1.47 )
233    
234 rtoy 1.48 #+(or ppc x86)
235     (macrolet
236     ((frob (sin cos tan)
237 rtoy 1.54 `(progn
238 rtoy 1.55 ;; In all of the routines below, we just compute the sum of
239     ;; y0 and y1 and use that as the (reduced) argument for the
240     ;; trig functions. This is slightly less accurate than what
241     ;; fdlibm does, which calls special functions using y0 and
242     ;; y1 separately, for greater accuracy. This isn't
243     ;; implemented, and some spot checks indicate that what we
244     ;; have here is accurate.
245 rtoy 1.49 ;;
246 rtoy 1.55 ;; For x86 with an fsin/fcos/fptan instruction, the pi/4 is
247     ;; probably too restrictive.
248 rtoy 1.48 (defun %sin (x)
249     (declare (double-float x))
250     (if (< (abs x) (/ pi 4))
251     (,sin x)
252     ;; Argument reduction needed
253 rtoy 1.54 (multiple-value-bind (n y0 y1)
254     (%ieee754-rem-pi/2 x)
255     (let ((reduced (+ y0 y1)))
256     (case (logand n 3)
257     (0 (,sin reduced))
258     (1 (,cos reduced))
259     (2 (- (,sin reduced)))
260     (3 (- (,cos reduced))))))))
261 rtoy 1.48 (defun %cos (x)
262     (declare (double-float x))
263     (if (< (abs x) (/ pi 4))
264     (,cos x)
265     ;; Argument reduction needed
266 rtoy 1.54 (multiple-value-bind (n y0 y1)
267     (%ieee754-rem-pi/2 x)
268     (let ((reduced (+ y0 y1)))
269     (case (logand n 3)
270     (0 (,cos reduced))
271     (1 (- (,sin reduced)))
272     (2 (- (,cos reduced)))
273     (3 (,sin reduced)))))))
274 rtoy 1.48 (defun %tan (x)
275     (declare (double-float x))
276     (if (< (abs x) (/ pi 4))
277     (,tan x)
278     ;; Argument reduction needed
279 rtoy 1.54 (multiple-value-bind (n y0 y1)
280     (%ieee754-rem-pi/2 x)
281     (let ((reduced (+ y0 y1)))
282     (if (evenp n)
283     (,tan reduced)
284     (- (/ (,tan reduced)))))))))))
285 rtoy 1.58 ;; Don't want %sin-quick and friends with sse2.
286     #+(and x86 (not sse2))
287 rtoy 1.48 (frob %sin-quick %cos-quick %tan-quick)
288 rtoy 1.58 #+(or ppc sse2)
289 rtoy 1.48 (frob %%sin %%cos %%tan))
290    
291    
292 wlott 1.1
293     ;;;; Power functions.
294    
295     (defun exp (number)
296 rtoy 1.60.2.2 _N"Return e raised to the power NUMBER."
297 wlott 1.1 (number-dispatch ((number number))
298     (handle-reals %exp number)
299     ((complex)
300     (* (exp (realpart number))
301     (cis (imagpart number))))))
302    
303     ;;; INTEXP -- Handle the rational base, integer power case.
304    
305     (defparameter *intexp-maximum-exponent* 10000)
306    
307 rtoy 1.57 (define-condition intexp-limit-error (error)
308     ((base :initarg :base :reader intexp-base)
309     (power :initarg :power :reader intexp-power))
310     (:report (lambda (condition stream)
311 rtoy 1.60.2.2 (format stream _"The absolute value of ~S exceeds limit ~S."
312 rtoy 1.57 (intexp-power condition)
313     *intexp-maximum-exponent*))))
314    
315 ram 1.6 ;;; This function precisely calculates base raised to an integral power. It
316     ;;; separates the cases by the sign of power, for efficiency reasons, as powers
317     ;;; can be calculated more efficiently if power is a positive integer. Values
318     ;;; of power are calculated as positive integers, and inverted if negative.
319     ;;;
320 wlott 1.1 (defun intexp (base power)
321 rtoy 1.45 ;; Handle the special case of 1^power. Maxima sometimes does this,
322     ;; and there's no need to cause a continuable error in this case.
323     ;; Should we also handle (-1)^power?
324     (when (eql base 1)
325     (return-from intexp base))
326    
327 wlott 1.1 (when (> (abs power) *intexp-maximum-exponent*)
328 rtoy 1.57 ;; Allow user the option to continue with calculation, possibly
329     ;; increasing the limit to the given power.
330     (restart-case
331     (error 'intexp-limit-error
332     :base base
333     :power power)
334     (continue ()
335 rtoy 1.60.2.2 :report _"Continue with calculation")
336 rtoy 1.57 (new-limit ()
337 rtoy 1.60.2.2 :report _"Continue with calculation, update limit"
338 rtoy 1.60 (setq *intexp-maximum-exponent* (abs power)))))
339 wlott 1.1 (cond ((minusp power)
340     (/ (intexp base (- power))))
341     ((eql base 2)
342     (ash 1 power))
343     (t
344     (do ((nextn (ash power -1) (ash power -1))
345     (total (if (oddp power) base 1)
346     (if (oddp power) (* base total) total)))
347     ((zerop nextn) total)
348     (setq base (* base base))
349     (setq power nextn)))))
350    
351    
352 ram 1.6 ;;; EXPT -- Public
353     ;;;
354     ;;; If an integer power of a rational, use INTEXP above. Otherwise, do
355     ;;; floating point stuff. If both args are real, we try %POW right off,
356     ;;; assuming it will return 0 if the result may be complex. If so, we call
357     ;;; COMPLEX-POW which directly computes the complex result. We also separate
358     ;;; the complex-real and real-complex cases from the general complex case.
359     ;;;
360 wlott 1.1 (defun expt (base power)
361 rtoy 1.60.2.2 _N"Returns BASE raised to the POWER."
362 wlott 1.1 (if (zerop power)
363 rtoy 1.40 ;; CLHS says that if the power is 0, the result is 1, subject to
364     ;; numeric contagion. But what happens if base is infinity or
365     ;; NaN? Do we silently return 1? For now, I think we should
366     ;; signal an error if the FP modes say so.
367     (let ((result (1+ (* base power))))
368     ;; If we get an NaN here, that means base*power above didn't
369     ;; produce 0 and FP traps were disabled, so we handle that
370     ;; here. Should this be a continuable restart?
371     (if (and (floatp result) (float-nan-p result))
372     (float 1 result)
373     result))
374 dtc 1.22 (labels (;; determine if the double float is an integer.
375     ;; 0 - not an integer
376     ;; 1 - an odd int
377     ;; 2 - an even int
378     (isint (ihi lo)
379     (declare (type (unsigned-byte 31) ihi)
380     (type (unsigned-byte 32) lo)
381     (optimize (speed 3) (safety 0)))
382     (let ((isint 0))
383     (declare (type fixnum isint))
384     (cond ((>= ihi #x43400000) ; exponent >= 53
385     (setq isint 2))
386     ((>= ihi #x3ff00000)
387     (let ((k (- (ash ihi -20) #x3ff))) ; exponent
388     (declare (type (mod 53) k))
389     (cond ((> k 20)
390     (let* ((shift (- 52 k))
391     (j (logand (ash lo (- shift))))
392     (j2 (ash j shift)))
393     (declare (type (mod 32) shift)
394     (type (unsigned-byte 32) j j2))
395     (when (= j2 lo)
396     (setq isint (- 2 (logand j 1))))))
397     ((= lo 0)
398     (let* ((shift (- 20 k))
399     (j (ash ihi (- shift)))
400     (j2 (ash j shift)))
401     (declare (type (mod 32) shift)
402     (type (unsigned-byte 31) j j2))
403     (when (= j2 ihi)
404     (setq isint (- 2 (logand j 1))))))))))
405     isint))
406     (real-expt (x y rtype)
407     (let ((x (coerce x 'double-float))
408     (y (coerce y 'double-float)))
409     (declare (double-float x y))
410     (let* ((x-hi (kernel:double-float-high-bits x))
411     (x-lo (kernel:double-float-low-bits x))
412     (x-ihi (logand x-hi #x7fffffff))
413     (y-hi (kernel:double-float-high-bits y))
414     (y-lo (kernel:double-float-low-bits y))
415     (y-ihi (logand y-hi #x7fffffff)))
416     (declare (type (signed-byte 32) x-hi y-hi)
417     (type (unsigned-byte 31) x-ihi y-ihi)
418     (type (unsigned-byte 32) x-lo y-lo))
419     ;; y==zero: x**0 = 1
420     (when (zerop (logior y-ihi y-lo))
421     (return-from real-expt (coerce 1d0 rtype)))
422     ;; +-NaN return x+y
423     (when (or (> x-ihi #x7ff00000)
424     (and (= x-ihi #x7ff00000) (/= x-lo 0))
425     (> y-ihi #x7ff00000)
426     (and (= y-ihi #x7ff00000) (/= y-lo 0)))
427     (return-from real-expt (coerce (+ x y) rtype)))
428     (let ((yisint (if (< x-hi 0) (isint y-ihi y-lo) 0)))
429     (declare (type fixnum yisint))
430     ;; special value of y
431     (when (and (zerop y-lo) (= y-ihi #x7ff00000))
432     ;; y is +-inf
433     (return-from real-expt
434     (cond ((and (= x-ihi #x3ff00000) (zerop x-lo))
435     ;; +-1**inf is NaN
436     (coerce (- y y) rtype))
437     ((>= x-ihi #x3ff00000)
438     ;; (|x|>1)**+-inf = inf,0
439     (if (>= y-hi 0)
440     (coerce y rtype)
441     (coerce 0 rtype)))
442     (t
443     ;; (|x|<1)**-,+inf = inf,0
444     (if (< y-hi 0)
445     (coerce (- y) rtype)
446     (coerce 0 rtype))))))
447    
448     (let ((abs-x (abs x)))
449     (declare (double-float abs-x))
450     ;; special value of x
451     (when (and (zerop x-lo)
452     (or (= x-ihi #x7ff00000) (zerop x-ihi)
453     (= x-ihi #x3ff00000)))
454     ;; x is +-0,+-inf,+-1
455     (let ((z (if (< y-hi 0)
456     (/ 1 abs-x) ; z = (1/|x|)
457     abs-x)))
458     (declare (double-float z))
459     (when (< x-hi 0)
460     (cond ((and (= x-ihi #x3ff00000) (zerop yisint))
461     ;; (-1)**non-int
462     (let ((y*pi (* y pi)))
463     (declare (double-float y*pi))
464     (return-from real-expt
465 dtc 1.24 (complex
466 dtc 1.22 (coerce (%cos y*pi) rtype)
467     (coerce (%sin y*pi) rtype)))))
468     ((= yisint 1)
469     ;; (x<0)**odd = -(|x|**odd)
470     (setq z (- z)))))
471     (return-from real-expt (coerce z rtype))))
472    
473     (if (>= x-hi 0)
474     ;; x>0
475     (coerce (kernel::%pow x y) rtype)
476     ;; x<0
477     (let ((pow (kernel::%pow abs-x y)))
478     (declare (double-float pow))
479     (case yisint
480     (1 ; Odd
481     (coerce (* -1d0 pow) rtype))
482     (2 ; Even
483     (coerce pow rtype))
484     (t ; Non-integer
485     (let ((y*pi (* y pi)))
486     (declare (double-float y*pi))
487 dtc 1.24 (complex
488 dtc 1.22 (coerce (* pow (%cos y*pi)) rtype)
489     (coerce (* pow (%sin y*pi)) rtype)))))))))))))
490 rtoy 1.56 ;; This is really messy and should be cleaned up. The easiest
491     ;; way to see if we're doing what we should is the macroexpand
492     ;; the number-dispatch and check each branch.
493     ;;
494     ;; We try to apply the rule of float precision contagion (CLHS
495     ;; 12.1.4.4): the result has the same precision has the most
496     ;; precise argument.
497 dtc 1.22 (number-dispatch ((base number) (power number))
498 rtoy 1.56 (((foreach fixnum (or bignum ratio) (complex rational))
499     integer)
500 dtc 1.22 (intexp base power))
501 rtoy 1.56 (((foreach single-float double-float)
502     rational)
503 dtc 1.22 (real-expt base power '(dispatch-type base)))
504     (((foreach fixnum (or bignum ratio) single-float)
505     (foreach ratio single-float))
506     (real-expt base power 'single-float))
507     (((foreach fixnum (or bignum ratio) single-float double-float)
508     double-float)
509     (real-expt base power 'double-float))
510     ((double-float single-float)
511     (real-expt base power 'double-float))
512 rtoy 1.46 #+double-double
513 rtoy 1.56 (((foreach fixnum (or bignum ratio) single-float double-float
514     double-double-float)
515 rtoy 1.46 double-double-float)
516     (dd-%pow (coerce base 'double-double-float) power))
517     #+double-double
518     ((double-double-float
519     (foreach fixnum (or bignum ratio) single-float double-float))
520     (dd-%pow base (coerce power 'double-double-float)))
521 rtoy 1.56 (((foreach (complex rational) (complex single-float) (complex double-float)
522     #+double-double (complex double-double-float))
523     rational)
524 dtc 1.22 (* (expt (abs base) power)
525     (cis (* power (phase base)))))
526 rtoy 1.56 #+double-double
527     ((double-double-float
528 dtc 1.22 complex)
529     (if (and (zerop base) (plusp (realpart power)))
530     (* base power)
531 rtoy 1.56 (exp (* power (* (log2 base 1w0) (log 2w0))))))
532     (((foreach fixnum (or bignum ratio) single-float double-float)
533     (foreach (complex double-float)))
534     ;; Result should have double-float accuracy. Use log2 in
535     ;; case the base won't fit in a double-float.
536     (if (and (zerop base) (plusp (realpart power)))
537     (* base power)
538     (exp (* power (* (log2 base) (log 2d0))))))
539     ((double-float
540     (foreach (complex rational) (complex single-float)))
541     (if (and (zerop base) (plusp (realpart power)))
542     (* base power)
543     (exp (* power (log base)))))
544     #+double-double
545     (((foreach fixnum (or bignum ratio) single-float double-float)
546     (foreach (complex double-double-float)))
547     ;; Result should have double-double-float accuracy. Use log2
548     ;; in case the base won't fit in a double-float.
549     (if (and (zerop base) (plusp (realpart power)))
550     (* base power)
551     (exp (* power (* (log2 base 1w0) (log 2w0))))))
552     (((foreach fixnum (or bignum ratio) single-float)
553     (foreach (complex single-float)))
554     (if (and (zerop base) (plusp (realpart power)))
555     (* base power)
556     (exp (* power (log base)))))
557     (((foreach (complex rational) (complex single-float))
558     (foreach single-float (complex single-float)))
559     (if (and (zerop base) (plusp (realpart power)))
560     (* base power)
561     (exp (* power (log base)))))
562     (((foreach (complex rational) (complex single-float))
563     (foreach double-float (complex double-float)))
564     (if (and (zerop base) (plusp (realpart power)))
565     (* base power)
566     (exp (* power (log (coerce base '(complex double-float)))))))
567     #+double-double
568     (((foreach (complex rational) (complex single-float))
569     (foreach double-double-float (complex double-double-float)))
570     (if (and (zerop base) (plusp (realpart power)))
571     (* base power)
572     (exp (* power (log (coerce base '(complex double-double-float)))))))
573     (((foreach (complex double-float))
574     (foreach single-float double-float (complex single-float)
575     (complex double-float)))
576     (if (and (zerop base) (plusp (realpart power)))
577     (* base power)
578 dtc 1.22 (exp (* power (log base)))))
579 rtoy 1.56 #+double-double
580     (((foreach (complex double-float))
581     (foreach double-double-float (complex double-double-float)))
582     (if (and (zerop base) (plusp (realpart power)))
583     (* base power)
584     (exp (* power (log (coerce base '(complex double-double-float)))))))
585     #+double-double
586     (((foreach (complex double-double-float))
587     (foreach float (complex float)))
588 dtc 1.22 (if (and (zerop base) (plusp (realpart power)))
589     (* base power)
590     (exp (* power (log base)))))))))
591 wlott 1.1
592 rtoy 1.52 ;; Log base 2 of a real number. The result is a either a double-float
593     ;; or double-double-float number (real or complex, as appropriate),
594     ;; depending on the type of FLOAT-TYPE.
595     (defun log2 (x &optional (float-type 1d0))
596     (labels ((log-of-2 (f)
597     ;; log(2), with the precision specified by the type of F
598     (number-dispatch ((f real))
599     ((double-float)
600 rtoy 1.51 #.(log 2d0))
601 rtoy 1.52 #+double-double
602     ((double-double-float)
603     #.(log 2w0))))
604     (log-2-pi (f)
605     ;; log(pi), with the precision specified by the type of F
606     (number-dispatch ((f real))
607     ((double-float)
608     #.(/ pi (log 2d0)))
609     #+double-double
610     ((double-double-float)
611     #.(/ dd-pi (log 2w0)))))
612     (log1p (x)
613     ;; log(1+x), with the precision specified by the type of
614     ;; X
615     (number-dispatch ((x real))
616     (((foreach single-float double-float))
617     (%log1p (float x 1d0)))
618     #+double-double
619     ((double-double-float)
620     (dd-%log1p x))))
621     (log2-bignum (bignum)
622 rtoy 1.51 ;; Write x = 2^n*f where 1/2 < f <= 1. Then log2(x) = n
623     ;; + log2(f).
624     ;;
625     ;; So we grab the top few bits of x and scale that
626     ;; appropriately, take the log of it and add it to n.
627     ;;
628     ;; Return n and log2(f) separately.
629     (if (minusp bignum)
630     (multiple-value-bind (n frac)
631     (log2-bignum (abs bignum))
632 rtoy 1.52 (values n (complex frac (log-2-pi float-type))))
633     (let ((n (integer-length bignum))
634     (float-bits (float-digits float-type)))
635     (if (< n float-bits)
636     (values 0 (log (float bignum float-type)
637     (float 2 float-type)))
638     (let ((exp (min float-bits n))
639     (f (ldb (byte float-bits
640     (max 0 (- n float-bits)))
641 rtoy 1.51 bignum)))
642 rtoy 1.52 (values n (log (scale-float (float f float-type) (- exp))
643     (float 2 float-type)))))))))
644 rtoy 1.51 (etypecase x
645     (float
646 rtoy 1.52 (/ (log (float x float-type)) (log-of-2 float-type)))
647 rtoy 1.51 (ratio
648     (let ((top (numerator x))
649     (bot (denominator x)))
650     ;; If the number of bits in the numerator and
651     ;; denominator are different, just use the fact
652     ;; log(x/y) = log(x) - log(y). But to preserve
653     ;; accuracy, we actually do
654     ;; (log2(x)-log2(y))/log2(e)).
655     ;;
656     ;; However, if the numerator and denominator have the
657     ;; same number of bits, implying the quotient is near
658     ;; one, we use log1p(x) = log(1+x). Since the number is
659     ;; rational, we don't lose precision subtracting 1 from
660     ;; it, and converting it to double-float is accurate.
661     (if (= (integer-length top)
662     (integer-length bot))
663 rtoy 1.52 (/ (log1p (float (- x 1) float-type))
664     (log-of-2 float-type))
665 rtoy 1.51 (multiple-value-bind (top-n top-frac)
666     (log2-bignum top)
667     (multiple-value-bind (bot-n bot-frac)
668     (log2-bignum bot)
669     (+ (- top-n bot-n)
670     (- top-frac bot-frac)))))))
671     (integer
672     (multiple-value-bind (n frac)
673     (log2-bignum x)
674     (+ n frac))))))
675 rtoy 1.43
676 wlott 1.1 (defun log (number &optional (base nil base-p))
677 rtoy 1.60.2.2 _N"Return the logarithm of NUMBER in the base BASE, which defaults to e."
678 wlott 1.1 (if base-p
679 toy 1.34 (cond ((zerop base)
680     ;; ANSI spec
681     base)
682 rtoy 1.43 ((and (realp number) (realp base))
683     ;; CLHS 12.1.4.1 says
684     ;;
685     ;; When rationals and floats are combined by a
686     ;; numerical function, the rational is first converted
687     ;; to a float of the same format.
688     ;;
689     ;; So assume this applies to floats as well convert all
690     ;; numbers to the largest float format before computing
691     ;; the log.
692     ;;
693     ;; This makes (log 17 10.0) = (log 17.0 10) and so on.
694     (number-dispatch ((number real) (base real))
695     ((double-float
696 rtoy 1.51 (foreach double-float single-float))
697     (/ (log2 number) (log2 base)))
698     (((foreach fixnum bignum ratio)
699     (foreach fixnum bignum ratio single-float))
700     (let* ((result (/ (log2 number) (log2 base))))
701     ;; Figure out the right result type
702     (if (realp result)
703     (coerce result 'single-float)
704     (coerce result '(complex single-float)))))
705     (((foreach fixnum bignum ratio)
706 rtoy 1.43 double-float)
707 rtoy 1.51 (/ (log2 number) (log2 base)))
708     ((single-float
709     (foreach fixnum bignum ratio))
710     (let* ((result (/ (log2 number) (log2 base))))
711     ;; Figure out the right result type
712     (if (realp result)
713     (coerce result 'single-float)
714     (coerce result '(complex single-float)))))
715     ((double-float
716     (foreach fixnum bignum ratio))
717     (/ (log2 number) (log2 base)))
718     ((single-float double-float)
719 rtoy 1.43 (/ (log (coerce number 'double-float)) (log base)))
720 rtoy 1.46 #+double-double
721     ((double-double-float
722 rtoy 1.51 (foreach fixnum bignum ratio))
723 rtoy 1.52 (/ (log2 number 1w0) (log2 base 1w0)))
724 rtoy 1.51 #+double-double
725     ((double-double-float
726     (foreach double-double-float double-float single-float))
727 rtoy 1.46 (/ (log number) (log (coerce base 'double-double-float))))
728     #+double-double
729 rtoy 1.51 (((foreach fixnum bignum ratio)
730     double-double-float)
731 rtoy 1.52 (/ (log2 number 1w0) (log2 base 1w0)))
732 rtoy 1.51 #+double-double
733     (((foreach double-float single-float)
734 rtoy 1.46 double-double-float)
735     (/ (log (coerce number 'double-double-float)) (log base)))
736 rtoy 1.51 (((foreach single-float)
737     (foreach single-float))
738 rtoy 1.43 ;; Converting everything to double-float helps the
739     ;; cases like (log 17 10) = (/ (log 17) (log 10)).
740     ;; This is usually handled above, but if we compute (/
741     ;; (log 17) (log 10)), we get a slightly different
742     ;; answer due to roundoff. This makes it a bit more
743     ;; consistent.
744     ;;
745     ;; FIXME: This probably needs more work.
746     (let ((result (/ (log (float number 1d0))
747     (log (float base 1d0)))))
748     (if (realp result)
749     (coerce result 'single-float)
750     (coerce result '(complex single-float)))))))
751 toy 1.34 (t
752 rtoy 1.43 ;; FIXME: This probably needs some work as well.
753 toy 1.34 (/ (log number) (log base))))
754 wlott 1.1 (number-dispatch ((number number))
755 toy 1.34 (((foreach fixnum bignum))
756 wlott 1.3 (if (minusp number)
757 toy 1.34 (complex (coerce (log (- number)) 'single-float)
758     (coerce pi 'single-float))
759 toy 1.37 (coerce (/ (log2 number) #.(log (exp 1d0) 2d0)) 'single-float)))
760 toy 1.34 ((ratio)
761     (if (minusp number)
762     (complex (coerce (log (- number)) 'single-float)
763     (coerce pi 'single-float))
764 toy 1.35 ;; What happens when the ratio is close to 1? We need to
765     ;; be careful to preserve accuracy.
766     (let ((top (numerator number))
767     (bot (denominator number)))
768     ;; If the number of bits in the numerator and
769     ;; denominator are different, just use the fact
770 toy 1.37 ;; log(x/y) = log(x) - log(y). But to preserve
771     ;; accuracy, we actually do
772     ;; (log2(x)-log2(y))/log2(e)).
773     ;;
774     ;; However, if the numerator and denominator have the
775     ;; same number of bits, implying the quotient is near
776     ;; one, we use log1p(x) = log(1+x). Since the number is
777     ;; rational, we don't lose precision subtracting 1 from
778     ;; it, and converting it to double-float is accurate.
779 toy 1.35 (if (= (integer-length top)
780 toy 1.36 (integer-length bot))
781 toy 1.35 (coerce (%log1p (coerce (- number 1) 'double-float))
782     'single-float)
783 toy 1.37 (coerce (/ (- (log2 top) (log2 bot))
784     #.(log (exp 1d0) 2d0))
785     'single-float)))))
786 ram 1.17 (((foreach single-float double-float))
787     ;; Is (log -0) -infinity (libm.a) or -infinity + i*pi (Kahan)?
788     ;; Since this doesn't seem to be an implementation issue
789     ;; I (pw) take the Kahan result.
790     (if (< (float-sign number)
791     (coerce 0 '(dispatch-type number)))
792     (complex (log (- number)) (coerce pi '(dispatch-type number)))
793     (coerce (%log (coerce number 'double-float))
794     '(dispatch-type number))))
795 rtoy 1.46 #+double-double
796     ((double-double-float)
797     (let ((hi (kernel:double-double-hi number)))
798     (if (< (float-sign hi) 0d0)
799     (complex (dd-%log (- number)) dd-pi)
800     (dd-%log number))))
801 ram 1.17 ((complex)
802     (complex-log number)))))
803 wlott 1.1
804     (defun sqrt (number)
805 rtoy 1.60.2.2 _N"Return the square root of NUMBER."
806 wlott 1.1 (number-dispatch ((number number))
807 ram 1.17 (((foreach fixnum bignum ratio))
808 wlott 1.3 (if (minusp number)
809 ram 1.17 (complex-sqrt number)
810 wlott 1.3 (coerce (%sqrt (coerce number 'double-float)) 'single-float)))
811 ram 1.17 (((foreach single-float double-float))
812 dtc 1.28 (if (minusp number)
813 ram 1.17 (complex-sqrt number)
814     (coerce (%sqrt (coerce number 'double-float))
815     '(dispatch-type number))))
816 rtoy 1.46 #+double-double
817     ((double-double-float)
818     (if (minusp number)
819     (dd-complex-sqrt number)
820     (multiple-value-bind (hi lo)
821     (c::sqrt-dd (kernel:double-double-hi number) (kernel:double-double-lo number))
822 rtoy 1.50 (kernel:%make-double-double-float hi lo))))
823 rtoy 1.46 ((complex)
824     (complex-sqrt number))))
825 wlott 1.1
826    
827     ;;;; Trigonometic and Related Functions
828    
829 wlott 1.2 (defun abs (number)
830 rtoy 1.60.2.2 _N"Returns the absolute value of the number."
831 wlott 1.2 (number-dispatch ((number number))
832 rtoy 1.46 (((foreach single-float double-float fixnum rational
833     #+double-double double-double-float))
834 wlott 1.2 (abs number))
835     ((complex)
836     (let ((rx (realpart number))
837     (ix (imagpart number)))
838     (etypecase rx
839     (rational
840     (sqrt (+ (* rx rx) (* ix ix))))
841     (single-float
842     (coerce (%hypot (coerce rx 'double-float)
843     (coerce ix 'double-float))
844     'single-float))
845     (double-float
846 rtoy 1.46 (%hypot rx ix))
847     #+double-double
848     (double-double-float
849 rtoy 1.50 (multiple-value-bind (abs^2 scale)
850     (dd-cssqs number)
851     (scale-float (sqrt abs^2) scale))))))))
852 wlott 1.1
853     (defun phase (number)
854 rtoy 1.60.2.2 _N"Returns the angle part of the polar representation of a complex number.
855 wlott 1.3 For complex numbers, this is (atan (imagpart number) (realpart number)).
856 wlott 1.1 For non-complex positive numbers, this is 0. For non-complex negative
857     numbers this is PI."
858 wlott 1.3 (etypecase number
859 ram 1.17 (rational
860 wlott 1.3 (if (minusp number)
861     (coerce pi 'single-float)
862     0.0f0))
863 ram 1.17 (single-float
864     (if (minusp (float-sign number))
865     (coerce pi 'single-float)
866     0.0f0))
867 wlott 1.3 (double-float
868 ram 1.17 (if (minusp (float-sign number))
869 wlott 1.3 (coerce pi 'double-float)
870     0.0d0))
871 rtoy 1.46 #+double-double
872     (double-double-float
873     (if (minusp (float-sign number))
874     dd-pi
875     0w0))
876 wlott 1.3 (complex
877     (atan (imagpart number) (realpart number)))))
878 wlott 1.1
879    
880     (defun sin (number)
881 rtoy 1.60.2.2 _N"Return the sine of NUMBER."
882 wlott 1.1 (number-dispatch ((number number))
883     (handle-reals %sin number)
884     ((complex)
885     (let ((x (realpart number))
886     (y (imagpart number)))
887 ram 1.17 (complex (* (sin x) (cosh y))
888     (* (cos x) (sinh y)))))))
889 wlott 1.1
890     (defun cos (number)
891 rtoy 1.60.2.2 _N"Return the cosine of NUMBER."
892 wlott 1.1 (number-dispatch ((number number))
893     (handle-reals %cos number)
894     ((complex)
895     (let ((x (realpart number))
896     (y (imagpart number)))
897 ram 1.17 (complex (* (cos x) (cosh y))
898     (- (* (sin x) (sinh y))))))))
899 wlott 1.1
900     (defun tan (number)
901 rtoy 1.60.2.2 _N"Return the tangent of NUMBER."
902 wlott 1.1 (number-dispatch ((number number))
903     (handle-reals %tan number)
904     ((complex)
905 ram 1.17 (complex-tan number))))
906 wlott 1.1
907 wlott 1.2 (defun cis (theta)
908 rtoy 1.60.2.2 _N"Return cos(Theta) + i sin(Theta), AKA exp(i Theta)."
909 wlott 1.2 (if (complexp theta)
910 rtoy 1.60.2.2 (error _"Argument to CIS is complex: ~S" theta)
911 wlott 1.2 (complex (cos theta) (sin theta))))
912 wlott 1.1
913     (defun asin (number)
914 rtoy 1.60.2.2 _N"Return the arc sine of NUMBER."
915 wlott 1.1 (number-dispatch ((number number))
916 wlott 1.3 ((rational)
917     (if (or (> number 1) (< number -1))
918     (complex-asin number)
919     (coerce (%asin (coerce number 'double-float)) 'single-float)))
920     (((foreach single-float double-float))
921 rtoy 1.42 (if (or (float-nan-p number)
922     (and (<= number (coerce 1 '(dispatch-type number)))
923     (>= number (coerce -1 '(dispatch-type number)))))
924 wlott 1.3 (coerce (%asin (coerce number 'double-float))
925 rtoy 1.42 '(dispatch-type number))
926 rtoy 1.46 (complex-asin number)))
927     #+double-double
928     ((double-double-float)
929     (if (or (float-nan-p number)
930     (and (<= number 1w0)
931     (>= number -1w0)))
932     (dd-%asin number)
933     (dd-complex-asin number)))
934 wlott 1.1 ((complex)
935 wlott 1.3 (complex-asin number))))
936 wlott 1.1
937     (defun acos (number)
938 rtoy 1.60.2.2 _N"Return the arc cosine of NUMBER."
939 wlott 1.1 (number-dispatch ((number number))
940 wlott 1.3 ((rational)
941     (if (or (> number 1) (< number -1))
942     (complex-acos number)
943     (coerce (%acos (coerce number 'double-float)) 'single-float)))
944     (((foreach single-float double-float))
945 rtoy 1.42 (if (or (float-nan-p number)
946     (and (<= number (coerce 1 '(dispatch-type number)))
947     (>= number (coerce -1 '(dispatch-type number)))))
948 wlott 1.3 (coerce (%acos (coerce number 'double-float))
949 rtoy 1.42 '(dispatch-type number))
950 rtoy 1.46 (complex-acos number)))
951     #+double-double
952     ((double-double-float)
953     (if (or (float-nan-p number)
954     (and (<= number 1w0)
955     (>= number -1w0)))
956     (dd-%acos number)
957     (complex-acos number)))
958 wlott 1.1 ((complex)
959 wlott 1.3 (complex-acos number))))
960 wlott 1.1
961    
962     (defun atan (y &optional (x nil xp))
963 rtoy 1.60.2.2 _N"Return the arc tangent of Y if X is omitted or Y/X if X is supplied."
964 wlott 1.1 (if xp
965 wlott 1.12 (flet ((atan2 (y x)
966 wlott 1.13 (declare (type double-float y x)
967     (values double-float))
968     (if (zerop x)
969     (if (zerop y)
970     (if (plusp (float-sign x))
971     y
972     (float-sign y pi))
973     (float-sign y (/ pi 2)))
974     (%atan2 y x))))
975 toy 1.33 ;; If X is given, both X and Y must be real numbers.
976     (number-dispatch ((y real) (x real))
977 wlott 1.13 ((double-float
978     (foreach double-float single-float fixnum bignum ratio))
979     (atan2 y (coerce x 'double-float)))
980     (((foreach single-float fixnum bignum ratio)
981     double-float)
982     (atan2 (coerce y 'double-float) x))
983     (((foreach single-float fixnum bignum ratio)
984     (foreach single-float fixnum bignum ratio))
985     (coerce (atan2 (coerce y 'double-float) (coerce x 'double-float))
986 rtoy 1.46 'single-float))
987     #+double-double
988     ((double-double-float
989     (foreach double-double-float double-float single-float fixnum bignum ratio))
990     (dd-%atan2 y (coerce x 'double-double-float)))
991     #+double-double
992     (((foreach double-float single-float fixnum bignum ratio)
993     double-double-float)
994     (dd-%atan2 (coerce y 'double-double-float) x))))
995 wlott 1.1 (number-dispatch ((y number))
996     (handle-reals %atan y)
997     ((complex)
998 ram 1.17 (complex-atan y)))))
999 wlott 1.3
1000 ram 1.17 (defun sinh (number)
1001 rtoy 1.60.2.2 _N"Return the hyperbolic sine of NUMBER."
1002 ram 1.17 (number-dispatch ((number number))
1003     (handle-reals %sinh number)
1004     ((complex)
1005     (let ((x (realpart number))
1006     (y (imagpart number)))
1007     (complex (* (sinh x) (cos y))
1008     (* (cosh x) (sin y)))))))
1009    
1010     (defun cosh (number)
1011 rtoy 1.60.2.2 _N"Return the hyperbolic cosine of NUMBER."
1012 ram 1.17 (number-dispatch ((number number))
1013     (handle-reals %cosh number)
1014     ((complex)
1015     (let ((x (realpart number))
1016     (y (imagpart number)))
1017     (complex (* (cosh x) (cos y))
1018     (* (sinh x) (sin y)))))))
1019    
1020 wlott 1.3 (defun tanh (number)
1021 rtoy 1.60.2.2 _N"Return the hyperbolic tangent of NUMBER."
1022 ram 1.17 (number-dispatch ((number number))
1023     (handle-reals %tanh number)
1024     ((complex)
1025     (complex-tanh number))))
1026 wlott 1.3
1027     (defun asinh (number)
1028 rtoy 1.60.2.2 _N"Return the hyperbolic arc sine of NUMBER."
1029 ram 1.17 (number-dispatch ((number number))
1030     (handle-reals %asinh number)
1031     ((complex)
1032     (complex-asinh number))))
1033 wlott 1.3
1034     (defun acosh (number)
1035 rtoy 1.60.2.2 _N"Return the hyperbolic arc cosine of NUMBER."
1036 ram 1.17 (number-dispatch ((number number))
1037     ((rational)
1038     ;; acosh is complex if number < 1
1039     (if (< number 1)
1040     (complex-acosh number)
1041     (coerce (%acosh (coerce number 'double-float)) 'single-float)))
1042     (((foreach single-float double-float))
1043     (if (< number (coerce 1 '(dispatch-type number)))
1044     (complex-acosh number)
1045     (coerce (%acosh (coerce number 'double-float))
1046     '(dispatch-type number))))
1047 rtoy 1.46 #+double-double
1048     ((double-double-float)
1049     (if (< number 1w0)
1050     (complex-acosh number)
1051     (dd-%acosh number)))
1052 ram 1.17 ((complex)
1053     (complex-acosh number))))
1054 wlott 1.3
1055     (defun atanh (number)
1056 rtoy 1.60.2.2 _N"Return the hyperbolic arc tangent of NUMBER."
1057 ram 1.17 (number-dispatch ((number number))
1058     ((rational)
1059     ;; atanh is complex if |number| > 1
1060     (if (or (> number 1) (< number -1))
1061     (complex-atanh number)
1062     (coerce (%atanh (coerce number 'double-float)) 'single-float)))
1063     (((foreach single-float double-float))
1064     (if (or (> number (coerce 1 '(dispatch-type number)))
1065     (< number (coerce -1 '(dispatch-type number))))
1066     (complex-atanh number)
1067     (coerce (%atanh (coerce number 'double-float))
1068     '(dispatch-type number))))
1069 rtoy 1.46 #+double-double
1070     ((double-double-float)
1071     (if (or (> number 1w0)
1072     (< number -1w0))
1073     (complex-atanh number)
1074     (dd-%atanh (coerce number 'double-double-float))))
1075 ram 1.17 ((complex)
1076     (complex-atanh number))))
1077 wlott 1.14
1078 toy 1.32 ;;; HP-UX does not supply a C version of log1p, so use the definition.
1079     ;;; We really need to fix this. The definition really loses big-time
1080     ;;; in roundoff as x gets small.
1081 wlott 1.14
1082     #+hpux
1083 pw 1.18 (declaim (inline %log1p))
1084 wlott 1.14 #+hpux
1085 pw 1.18 (defun %log1p (number)
1086     (declare (double-float number)
1087     (optimize (speed 3) (safety 0)))
1088     (the double-float (log (the (double-float 0d0) (+ number 1d0)))))
1089 ram 1.17
1090    
1091     ;;;;
1092     ;;;; This is a set of routines that implement many elementary
1093     ;;;; transcendental functions as specified by ANSI Common Lisp. The
1094     ;;;; implementation is based on Kahan's paper.
1095     ;;;;
1096     ;;;; I believe I have accurately implemented the routines and are
1097     ;;;; correct, but you may want to check for your self.
1098     ;;;;
1099     ;;;; These functions are written for CMU Lisp and take advantage of
1100     ;;;; some of the features available there. It may be possible,
1101     ;;;; however, to port this to other Lisps.
1102     ;;;;
1103     ;;;; Some functions are significantly more accurate than the original
1104     ;;;; definitions in CMU Lisp. In fact, some functions in CMU Lisp
1105     ;;;; give the wrong answer like (acos #c(-2.0 0.0)), where the true
1106     ;;;; answer is pi + i*log(2-sqrt(3)).
1107     ;;;;
1108     ;;;; All of the implemented functions will take any number for an
1109     ;;;; input, but the result will always be a either a complex
1110     ;;;; single-float or a complex double-float.
1111     ;;;;
1112     ;;;; General functions
1113     ;;;; complex-sqrt
1114     ;;;; complex-log
1115     ;;;; complex-atanh
1116     ;;;; complex-tanh
1117     ;;;; complex-acos
1118     ;;;; complex-acosh
1119     ;;;; complex-asin
1120     ;;;; complex-asinh
1121     ;;;; complex-atan
1122     ;;;; complex-tan
1123     ;;;;
1124     ;;;; Utility functions:
1125     ;;;; scalb logb
1126     ;;;;
1127     ;;;; Internal functions:
1128     ;;;; square coerce-to-complex-type cssqs complex-log-scaled
1129     ;;;;
1130     ;;;;
1131     ;;;; Please send any bug reports, comments, or improvements to Raymond
1132     ;;;; Toy at toy@rtp.ericsson.se.
1133     ;;;;
1134     ;;;; References
1135     ;;;;
1136     ;;;; Kahan, W. "Branch Cuts for Complex Elementary Functions, or Much
1137     ;;;; Ado About Nothing's Sign Bit" in Iserles and Powell (eds.) "The
1138     ;;;; State of the Art in Numerical Analysis", pp. 165-211, Clarendon
1139     ;;;; Press, 1987
1140     ;;;;
1141 toy 1.32
1142 ram 1.17 (declaim (inline square))
1143     (defun square (x)
1144 rtoy 1.46 (declare (float x))
1145 ram 1.17 (* x x))
1146    
1147     ;; If you have these functions in libm, perhaps they should be used
1148     ;; instead of these Lisp versions. These versions are probably good
1149     ;; enough, especially since they are portable.
1150    
1151 dtc 1.19 (declaim (inline scalb))
1152 ram 1.17 (defun scalb (x n)
1153 rtoy 1.60.2.2 _N"Compute 2^N * X without compute 2^N first (use properties of the
1154 ram 1.17 underlying floating-point format"
1155 rtoy 1.46 (declare (type float x)
1156 dtc 1.19 (type double-float-exponent n))
1157 ram 1.17 (scale-float x n))
1158    
1159 toy 1.32 (declaim (inline logb-finite))
1160     (defun logb-finite (x)
1161 rtoy 1.60.2.2 _N"Same as logb but X is not infinity and non-zero and not a NaN, so
1162 toy 1.32 that we can always return an integer"
1163 rtoy 1.46 (declare (type float x))
1164 toy 1.32 (multiple-value-bind (signif expon sign)
1165     (decode-float x)
1166     (declare (ignore signif sign))
1167     ;; decode-float is almost right, except that the exponent
1168     ;; is off by one
1169     (1- expon)))
1170    
1171 ram 1.17 (defun logb (x)
1172 rtoy 1.60.2.2 _N"Compute an integer N such that 1 <= |2^(-N) * x| < 2.
1173 ram 1.17 For the special cases, the following values are used:
1174    
1175     x logb
1176     NaN NaN
1177     +/- infinity +infinity
1178     0 -infinity
1179     "
1180 rtoy 1.46 (declare (type float x))
1181 dtc 1.26 (cond ((float-nan-p x)
1182 ram 1.17 x)
1183     ((float-infinity-p x)
1184     #.ext:double-float-positive-infinity)
1185     ((zerop x)
1186     ;; The answer is negative infinity, but we are supposed to
1187 toy 1.32 ;; signal divide-by-zero, so do the actual division
1188 rtoy 1.46 (/ -1 x)
1189 ram 1.17 )
1190     (t
1191 toy 1.32 (logb-finite x))))
1192    
1193    
1194 ram 1.17
1195     ;; This function is used to create a complex number of the appropriate
1196     ;; type.
1197    
1198     (declaim (inline coerce-to-complex-type))
1199     (defun coerce-to-complex-type (x y z)
1200 rtoy 1.60.2.2 _N"Create complex number with real part X and imaginary part Y such that
1201 ram 1.17 it has the same type as Z. If Z has type (complex rational), the X
1202     and Y are coerced to single-float."
1203     (declare (double-float x y)
1204 toy 1.32 (number z)
1205     (optimize (extensions:inhibit-warnings 3)))
1206     (if (typep (realpart z) 'double-float)
1207 ram 1.17 (complex x y)
1208     ;; Convert anything that's not a double-float to a single-float.
1209 toy 1.32 (complex (float x 1f0)
1210     (float y 1f0))))
1211 ram 1.17
1212     (defun cssqs (z)
1213 dtc 1.21 ;; Compute |(x+i*y)/2^k|^2 scaled to avoid over/underflow. The
1214     ;; result is r + i*k, where k is an integer.
1215 ram 1.17
1216     ;; Save all FP flags
1217 dtc 1.21 (let ((x (float (realpart z) 1d0))
1218 toy 1.32 (y (float (imagpart z) 1d0)))
1219 dtc 1.21 ;; Would this be better handled using an exception handler to
1220     ;; catch the overflow or underflow signal? For now, we turn all
1221     ;; traps off and look at the accrued exceptions to see if any
1222     ;; signal would have been raised.
1223     (with-float-traps-masked (:underflow :overflow)
1224 toy 1.32 (let ((rho (+ (square x) (square y))))
1225     (declare (optimize (speed 3) (space 0)))
1226     (cond ((and (or (float-nan-p rho)
1227     (float-infinity-p rho))
1228     (or (float-infinity-p (abs x))
1229     (float-infinity-p (abs y))))
1230     (values ext:double-float-positive-infinity 0))
1231     ((let ((threshold #.(/ least-positive-double-float
1232     double-float-epsilon))
1233     (traps (ldb vm::float-sticky-bits
1234     (vm:floating-point-modes))))
1235     ;; Overflow raised or (underflow raised and rho <
1236     ;; lambda/eps)
1237     (or (not (zerop (logand vm:float-overflow-trap-bit traps)))
1238     (and (not (zerop (logand vm:float-underflow-trap-bit traps)))
1239     (< rho threshold))))
1240     ;; If we're here, neither x nor y are infinity and at
1241     ;; least one is non-zero.. Thus logb returns a nice
1242     ;; integer.
1243     (let ((k (- (logb-finite (max (abs x) (abs y))))))
1244     (values (+ (square (scalb x k))
1245     (square (scalb y k)))
1246     (- k))))
1247     (t
1248     (values rho 0)))))))
1249 ram 1.17
1250     (defun complex-sqrt (z)
1251 rtoy 1.60.2.2 _N"Principle square root of Z
1252 ram 1.17
1253     Z may be any number, but the result is always a complex."
1254     (declare (number z))
1255 rtoy 1.46 #+double-double
1256     (when (typep z '(or double-double-float (complex double-double-float)))
1257     (return-from complex-sqrt (dd-complex-sqrt z)))
1258 ram 1.17 (multiple-value-bind (rho k)
1259     (cssqs z)
1260 toy 1.32 (declare (type (or (member 0d0) (double-float 0d0)) rho)
1261     (type fixnum k))
1262 ram 1.17 (let ((x (float (realpart z) 1.0d0))
1263     (y (float (imagpart z) 1.0d0))
1264     (eta 0d0)
1265     (nu 0d0))
1266     (declare (double-float x y eta nu))
1267    
1268 toy 1.32 (locally
1269     ;; space 0 to get maybe-inline functions inlined.
1270     (declare (optimize (speed 3) (space 0)))
1271    
1272 rtoy 1.53 (if (not (locally (declare (optimize (inhibit-warnings 3)))
1273     (float-nan-p x)))
1274 toy 1.32 (setf rho (+ (scalb (abs x) (- k)) (sqrt rho))))
1275    
1276     (cond ((oddp k)
1277     (setf k (ash k -1)))
1278     (t
1279     (setf k (1- (ash k -1)))
1280     (setf rho (+ rho rho))))
1281    
1282     (setf rho (scalb (sqrt rho) k))
1283    
1284     (setf eta rho)
1285     (setf nu y)
1286    
1287     (when (/= rho 0d0)
1288     (when (not (float-infinity-p (abs nu)))
1289     (setf nu (/ (/ nu rho) 2d0)))
1290     (when (< x 0d0)
1291     (setf eta (abs nu))
1292     (setf nu (float-sign y rho))))
1293     (coerce-to-complex-type eta nu z)))))
1294 ram 1.17
1295     (defun complex-log-scaled (z j)
1296 rtoy 1.60.2.2 _N"Compute log(2^j*z).
1297 ram 1.17
1298     This is for use with J /= 0 only when |z| is huge."
1299     (declare (number z)
1300     (fixnum j))
1301     ;; The constants t0, t1, t2 should be evaluated to machine
1302     ;; precision. In addition, Kahan says the accuracy of log1p
1303     ;; influences the choices of these constants but doesn't say how to
1304     ;; choose them. We'll just assume his choices matches our
1305     ;; implementation of log1p.
1306     (let ((t0 #.(/ 1 (sqrt 2.0d0)))
1307     (t1 1.2d0)
1308     (t2 3d0)
1309     (ln2 #.(log 2d0))
1310     (x (float (realpart z) 1.0d0))
1311     (y (float (imagpart z) 1.0d0)))
1312     (multiple-value-bind (rho k)
1313     (cssqs z)
1314 toy 1.32 (declare (optimize (speed 3)))
1315 ram 1.17 (let ((beta (max (abs x) (abs y)))
1316     (theta (min (abs x) (abs y))))
1317 toy 1.32 (coerce-to-complex-type (if (and (zerop k)
1318     (< t0 beta)
1319     (or (<= beta t1)
1320     (< rho t2)))
1321     (/ (%log1p (+ (* (- beta 1.0d0)
1322     (+ beta 1.0d0))
1323     (* theta theta)))
1324     2d0)
1325     (+ (/ (log rho) 2d0)
1326     (* (+ k j) ln2)))
1327     (atan y x)
1328     z)))))
1329 ram 1.17
1330     (defun complex-log (z)
1331 rtoy 1.60.2.2 _N"Log of Z = log |Z| + i * arg Z
1332 ram 1.17
1333     Z may be any number, but the result is always a complex."
1334     (declare (number z))
1335 rtoy 1.46 #+double-double
1336     (when (typep z '(or double-double-float (complex double-double-float)))
1337     (return-from complex-log (dd-complex-log-scaled z 0)))
1338 ram 1.17 (complex-log-scaled z 0))
1339    
1340     ;; Let us note the following "strange" behavior. atanh 1.0d0 is
1341     ;; +infinity, but the following code returns approx 176 + i*pi/4. The
1342     ;; reason for the imaginary part is caused by the fact that arg i*y is
1343     ;; never 0 since we have positive and negative zeroes.
1344    
1345     (defun complex-atanh (z)
1346 rtoy 1.60.2.2 _N"Compute atanh z = (log(1+z) - log(1-z))/2"
1347 ram 1.17 (declare (number z))
1348 rtoy 1.46 #+double-double
1349     (when (typep z '(or double-double-float (complex double-double-float)))
1350     (return-from complex-atanh (dd-complex-atanh z)))
1351    
1352 rtoy 1.41 (if (and (realp z) (< z -1))
1353     ;; atanh is continuous in quadrant III in this case.
1354     (complex-atanh (complex z -0f0))
1355     (let* ( ;; Constants
1356     (theta (/ (sqrt most-positive-double-float) 4.0d0))
1357     (rho (/ 4.0d0 (sqrt most-positive-double-float)))
1358     (half-pi (/ pi 2.0d0))
1359     (rp (float (realpart z) 1.0d0))
1360     (beta (float-sign rp 1.0d0))
1361     (x (* beta rp))
1362     (y (* beta (- (float (imagpart z) 1.0d0))))
1363     (eta 0.0d0)
1364     (nu 0.0d0))
1365     ;; Shouldn't need this declare.
1366     (declare (double-float x y))
1367     (locally
1368     (declare (optimize (speed 3)))
1369     (cond ((or (> x theta)
1370     (> (abs y) theta))
1371     ;; To avoid overflow...
1372     (setf nu (float-sign y half-pi))
1373     ;; eta is real part of 1/(x + iy). This is x/(x^2+y^2),
1374     ;; which can cause overflow. Arrange this computation so
1375     ;; that it won't overflow.
1376     (setf eta (let* ((x-bigger (> x (abs y)))
1377     (r (if x-bigger (/ y x) (/ x y)))
1378     (d (+ 1.0d0 (* r r))))
1379     (if x-bigger
1380     (/ (/ x) d)
1381     (/ (/ r y) d)))))
1382     ((= x 1.0d0)
1383     ;; Should this be changed so that if y is zero, eta is set
1384     ;; to +infinity instead of approx 176? In any case
1385     ;; tanh(176) is 1.0d0 within working precision.
1386     (let ((t1 (+ 4d0 (square y)))
1387     (t2 (+ (abs y) rho)))
1388     (setf eta (log (/ (sqrt (sqrt t1))
1389     (sqrt t2))))
1390     (setf nu (* 0.5d0
1391     (float-sign y
1392     (+ half-pi (atan (* 0.5d0 t2))))))))
1393     (t
1394     (let ((t1 (+ (abs y) rho)))
1395     ;; Normal case using log1p(x) = log(1 + x)
1396     (setf eta (* 0.25d0
1397     (%log1p (/ (* 4.0d0 x)
1398     (+ (square (- 1.0d0 x))
1399     (square t1))))))
1400     (setf nu (* 0.5d0
1401     (atan (* 2.0d0 y)
1402     (- (* (- 1.0d0 x)
1403     (+ 1.0d0 x))
1404     (square t1))))))))
1405     (coerce-to-complex-type (* beta eta)
1406     (- (* beta nu))
1407     z)))))
1408 ram 1.17
1409     (defun complex-tanh (z)
1410 rtoy 1.60.2.2 _N"Compute tanh z = sinh z / cosh z"
1411 ram 1.17 (declare (number z))
1412 rtoy 1.46 #+double-double
1413     (when (typep z '(or double-double-float (complex double-double-float)))
1414     (return-from complex-tanh (dd-complex-tanh z)))
1415    
1416 ram 1.17 (let ((x (float (realpart z) 1.0d0))
1417     (y (float (imagpart z) 1.0d0)))
1418 toy 1.32 (locally
1419     ;; space 0 to get maybe-inline functions inlined
1420     (declare (optimize (speed 3) (space 0)))
1421     (cond ((> (abs x)
1422     #-(or linux hpux) #.(/ (%asinh most-positive-double-float) 4d0)
1423     ;; This is more accurate under linux.
1424     #+(or linux hpux) #.(/ (+ (%log 2.0d0)
1425     (%log most-positive-double-float)) 4d0))
1426     (coerce-to-complex-type (float-sign x)
1427     (float-sign y) z))
1428     (t
1429     (let* ((tv (%tan y))
1430     (beta (+ 1.0d0 (* tv tv)))
1431     (s (sinh x))
1432     (rho (sqrt (+ 1.0d0 (* s s)))))
1433     (if (float-infinity-p (abs tv))
1434     (coerce-to-complex-type (/ rho s)
1435     (/ tv)
1436     z)
1437     (let ((den (+ 1.0d0 (* beta s s))))
1438     (coerce-to-complex-type (/ (* beta rho s)
1439     den)
1440     (/ tv den)
1441     z)))))))))
1442 ram 1.17
1443     ;; Kahan says we should only compute the parts needed. Thus, the
1444     ;; realpart's below should only compute the real part, not the whole
1445     ;; complex expression. Doing this can be important because we may get
1446     ;; spurious signals that occur in the part that we are not using.
1447     ;;
1448     ;; However, we take a pragmatic approach and just use the whole
1449     ;; expression.
1450    
1451     ;; NOTE: The formula given by Kahan is somewhat ambiguous in whether
1452     ;; it's the conjugate of the square root or the square root of the
1453     ;; conjugate. This needs to be checked.
1454    
1455     ;; I checked. It doesn't matter because (conjugate (sqrt z)) is the
1456     ;; same as (sqrt (conjugate z)) for all z. This follows because
1457     ;;
1458     ;; (conjugate (sqrt z)) = exp(0.5*log |z|)*exp(-0.5*j*arg z).
1459     ;;
1460     ;; (sqrt (conjugate z)) = exp(0.5*log|z|)*exp(0.5*j*arg conj z)
1461     ;;
1462 toy 1.35 ;; and these two expressions are equal if and only if arg conj z =
1463     ;; -arg z, which is clearly true for all z.
1464 ram 1.17
1465 rtoy 1.38 ;; NOTE: The rules of Common Lisp says that if you mix a real with a
1466     ;; complex, the real is converted to a complex before performing the
1467     ;; operation. However, Kahan says in this paper (pg 176):
1468     ;;
1469     ;; (iii) Careless handling can turn infinity or the sign of zero into
1470     ;; misinformation that subsequently disappears leaving behind
1471     ;; only a plausible but incorrect result. That is why compilers
1472     ;; must not transform z-1 into z-(1+i*0), as we have seen above,
1473     ;; nor -(-x-x^2) into (x+x^2), as we shall see below, lest a
1474     ;; subsequent logarithm or square root produce a non-zero
1475     ;; imaginary part whose sign is opposite to what was intended.
1476     ;;
1477     ;; The interesting examples are too long and complicated to reproduce
1478     ;; here. We refer the reader to his paper.
1479     ;;
1480     ;; The functions below are intended to handle the cases where a real
1481     ;; is mixed with a complex and we don't want CL complex contagion to
1482     ;; occur..
1483    
1484     (declaim (inline 1+z 1-z z-1 z+1))
1485     (defun 1+z (z)
1486     (complex (+ 1 (realpart z)) (imagpart z)))
1487     (defun 1-z (z)
1488     (complex (- 1 (realpart z)) (- (imagpart z))))
1489     (defun z-1 (z)
1490     (complex (- (realpart z) 1) (imagpart z)))
1491     (defun z+1 (z)
1492     (complex (+ (realpart z) 1) (imagpart z)))
1493    
1494 ram 1.17 (defun complex-acos (z)
1495 rtoy 1.60.2.2 _N"Compute acos z = pi/2 - asin z
1496 ram 1.17
1497     Z may be any number, but the result is always a complex."
1498     (declare (number z))
1499 rtoy 1.46 #+double-double
1500     (when (typep z '(or double-double-float (complex double-double-float)))
1501     (return-from complex-acos (dd-complex-acos z)))
1502 rtoy 1.41 (if (and (realp z) (> z 1))
1503     ;; acos is continuous in quadrant IV in this case.
1504     (complex-acos (complex z -0f0))
1505     (let ((sqrt-1+z (complex-sqrt (1+z z)))
1506     (sqrt-1-z (complex-sqrt (1-z z))))
1507     (with-float-traps-masked (:divide-by-zero)
1508     (complex (* 2 (atan (/ (realpart sqrt-1-z)
1509     (realpart sqrt-1+z))))
1510     (asinh (imagpart (* (conjugate sqrt-1+z)
1511     sqrt-1-z))))))))
1512 ram 1.17
1513     (defun complex-acosh (z)
1514 rtoy 1.60.2.2 _N"Compute acosh z = 2 * log(sqrt((z+1)/2) + sqrt((z-1)/2))
1515 ram 1.17
1516     Z may be any number, but the result is always a complex."
1517     (declare (number z))
1518 rtoy 1.38 (let ((sqrt-z-1 (complex-sqrt (z-1 z)))
1519     (sqrt-z+1 (complex-sqrt (z+1 z))))
1520 dtc 1.21 (with-float-traps-masked (:divide-by-zero)
1521     (complex (asinh (realpart (* (conjugate sqrt-z-1)
1522     sqrt-z+1)))
1523     (* 2 (atan (/ (imagpart sqrt-z-1)
1524     (realpart sqrt-z+1))))))))
1525 ram 1.17
1526    
1527     (defun complex-asin (z)
1528 rtoy 1.60.2.2 _N"Compute asin z = asinh(i*z)/i
1529 ram 1.17
1530     Z may be any number, but the result is always a complex."
1531     (declare (number z))
1532 rtoy 1.46 #+double-double
1533     (when (typep z '(or double-double-float (complex double-double-float)))
1534     (return-from complex-asin (dd-complex-asin z)))
1535 rtoy 1.41 (if (and (realp z) (> z 1))
1536     ;; asin is continuous in quadrant IV in this case.
1537     (complex-asin (complex z -0f0))
1538     (let ((sqrt-1-z (complex-sqrt (1-z z)))
1539     (sqrt-1+z (complex-sqrt (1+z z))))
1540     (with-float-traps-masked (:divide-by-zero)
1541     (complex (atan (/ (realpart z)
1542     (realpart (* sqrt-1-z sqrt-1+z))))
1543     (asinh (imagpart (* (conjugate sqrt-1-z)
1544     sqrt-1+z))))))))
1545 ram 1.17
1546     (defun complex-asinh (z)
1547 rtoy 1.60.2.2 _N"Compute asinh z = log(z + sqrt(1 + z*z))
1548 ram 1.17
1549     Z may be any number, but the result is always a complex."
1550     (declare (number z))
1551     ;; asinh z = -i * asin (i*z)
1552 rtoy 1.46 #+double-double
1553     (when (typep z '(or double-double-float (complex double-double-float)))
1554     (return-from complex-asinh (dd-complex-asinh z)))
1555 ram 1.17 (let* ((iz (complex (- (imagpart z)) (realpart z)))
1556     (result (complex-asin iz)))
1557     (complex (imagpart result)
1558     (- (realpart result)))))
1559    
1560     (defun complex-atan (z)
1561 rtoy 1.60.2.2 _N"Compute atan z = atanh (i*z) / i
1562 ram 1.17
1563     Z may be any number, but the result is always a complex."
1564     (declare (number z))
1565     ;; atan z = -i * atanh (i*z)
1566 rtoy 1.46 #+double-double
1567     (when (typep z '(or double-double-float (complex double-double-float)))
1568     (return-from complex-atan (dd-complex-atan z)))
1569 ram 1.17 (let* ((iz (complex (- (imagpart z)) (realpart z)))
1570     (result (complex-atanh iz)))
1571     (complex (imagpart result)
1572     (- (realpart result)))))
1573    
1574     (defun complex-tan (z)
1575 rtoy 1.60.2.2 _N"Compute tan z = -i * tanh(i * z)
1576 ram 1.17
1577     Z may be any number, but the result is always a complex."
1578     (declare (number z))
1579     ;; tan z = -i * tanh(i*z)
1580 rtoy 1.46 #+double-double
1581     (when (typep z '(or double-double-float (complex double-double-float)))
1582     (return-from complex-tan (dd-complex-tan z)))
1583 ram 1.17 (let* ((iz (complex (- (imagpart z)) (realpart z)))
1584     (result (complex-tanh iz)))
1585     (complex (imagpart result)
1586     (- (realpart result)))))
1587 toy 1.32

  ViewVC Help
Powered by ViewVC 1.1.5