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

Contents of /src/code/irrat.lisp

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.5