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

Contents of /src/code/irrat.lisp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.56 - (hide annotations)
Wed Oct 8 17:03:07 2008 UTC (5 years, 6 months ago) by rtoy
Branch: MAIN
Changes since 1.55: +75 -9 lines
Ticket #24.

Take the rule of float precision contagion (CLHS 12.1.4.4) to also
mean that the result should be as accurate as the most accurate
argument.  Effectively, all args are coerced to the highest precision
first before computing expt.

There's a simple test program to check that every case is covered with
the expected precision.  (I think).

(defun test-expt ()
  (dolist (base '(2 2f0 2d0 2w0) t)
    (dolist (power '(1/2 .5f0 .5d0 .5w0))
      (flet ((test-it (b p a eps expected)
	       (let* ((res (expt b p))
		      (absdiff (abs (- res a))))
		 (unless
		     (or (typep (realpart res) (type-of expected))
			 (<= absdiff (* 10 eps)))
		   (format t "FAILED: ~A^~A = ~A (~A)~%" b p res a)))))

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

  ViewVC Help
Powered by ViewVC 1.1.5