1 ;;;; This file contains all the irrational functions. (Actually, most
2 ;;;; of the work is done by calling out to C.)
4 ;;;; This software is part of the SBCL system. See the README file for
7 ;;;; This software is derived from the CMU CL system, which was
8 ;;;; written at Carnegie Mellon University and released into the
9 ;;;; public domain. The software is in the public domain and is
10 ;;;; provided with absolutely no warranty. See the COPYING and CREDITS
11 ;;;; files for more information.
13 (in-package "SB!KERNEL")
15 ;;;; miscellaneous constants, utility functions, and macros
18 #!+long-float
3.14159265358979323846264338327950288419716939937511l0
19 #!-long-float
3.14159265358979323846264338327950288419716939937511d0
)
21 ;;; Make these INLINE, since the call to C is at least as compact as a
22 ;;; Lisp call, and saves number consing to boot.
23 (eval-when (:compile-toplevel
:execute
)
25 (sb!xc
:defmacro def-math-rtn
(name num-args
)
26 (let ((function (symbolicate "%" (string-upcase name
)))
27 (args (loop for i below num-args
28 collect
(intern (format nil
"ARG~D" i
)))))
30 (declaim (inline ,function
))
31 (defun ,function
,args
34 (function double-float
35 ,@(loop repeat num-args
36 collect
'double-float
)))
39 (defun handle-reals (function var
)
40 `((((foreach fixnum single-float bignum ratio
))
41 (coerce (,function
(coerce ,var
'double-float
)) 'single-float
))
47 #!+x86
;; for constant folding
48 (macrolet ((def (name ll
)
49 `(defun ,name
,ll
(,name
,@ll
))))
62 #!+(or x86-64 arm-vfp
) ;; for constant folding
63 (macrolet ((def (name ll
)
64 `(defun ,name
,ll
(,name
,@ll
))))
67 ;;;; stubs for the Unix math library
69 ;;;; Many of these are unnecessary on the X86 because they're built
73 #!-x86
(def-math-rtn "sin" 1)
74 #!-x86
(def-math-rtn "cos" 1)
75 #!-x86
(def-math-rtn "tan" 1)
76 #!-x86
(def-math-rtn "atan" 1)
77 #!-x86
(def-math-rtn "atan2" 2)
80 (def-math-rtn "acos" 1)
81 (def-math-rtn "asin" 1)
82 (def-math-rtn "cosh" 1)
83 (def-math-rtn "sinh" 1)
84 (def-math-rtn "tanh" 1)
87 (def-math-rtn "asinh" 1)
88 (def-math-rtn "acosh" 1)
89 (def-math-rtn "atanh" 1)))
94 (declaim (inline %asin
))
96 (%atan
(/ number
(sqrt (- 1 (* number number
))))))
97 (declaim (inline %acos
))
99 (- (/ pi
2) (%asin number
)))
100 (declaim (inline %cosh
))
101 (defun %cosh
(number)
102 (/ (+ (exp number
) (exp (- number
))) 2))
103 (declaim (inline %sinh
))
104 (defun %sinh
(number)
105 (/ (- (exp number
) (exp (- number
))) 2))
106 (declaim (inline %tanh
))
107 (defun %tanh
(number)
108 (/ (%sinh number
) (%cosh number
))))
109 (declaim (inline %asinh
))
110 (defun %asinh
(number)
111 (log (+ number
(sqrt (+ (* number number
) 1.0d0
))) #.
(exp 1.0d0
)))
112 (declaim (inline %acosh
))
113 (defun %acosh
(number)
114 (log (+ number
(sqrt (- (* number number
) 1.0d0
))) #.
(exp 1.0d0
)))
115 (declaim (inline %atanh
))
116 (defun %atanh
(number)
117 (let ((ratio (/ (+ 1 number
) (- 1 number
))))
118 ;; Were we effectively zero?
121 (/ (log ratio
#.
(exp 1.0d0
)) 2.0d0
)))))
123 ;;; exponential and logarithmic
124 #!-x86
(def-math-rtn "exp" 1)
125 #!-x86
(def-math-rtn "log" 1)
126 #!-x86
(def-math-rtn "log10" 1)
127 #!-
(and win32 x86
) (def-math-rtn "pow" 2)
128 #!-
(or x86 x86-64 arm-vfp
) (def-math-rtn "sqrt" 1)
129 #!-win32
(def-math-rtn "hypot" 2)
130 #!-x86
(def-math-rtn "log1p" 1)
134 ;; This is written in a peculiar way to avoid overflow. Note that in
135 ;; sqrt(x^2 + y^2), either square or the sum can overflow.
137 ;; Factoring x^2 out of sqrt(x^2 + y^2) gives us the expression
138 ;; |x|sqrt(1 + (y/x)^2), which, assuming |x| >= |y|, can only overflow
139 ;; if |x| is sufficiently large.
141 ;; The ZEROP test suffices (y is non-negative) to guard against
142 ;; divisions by zero: x >= y > 0.
143 (declaim (inline %hypot
))
145 (declare (type double-float x y
))
153 (* x
(sqrt (1+ (* y
/x y
/x
)))))))))
159 "Return e raised to the power NUMBER."
160 (number-dispatch ((number number
))
161 (handle-reals %exp number
)
163 (* (exp (realpart number
))
164 (cis (imagpart number
))))))
166 ;;; INTEXP -- Handle the rational base, integer power case.
168 (declaim (type (or integer null
) *intexp-maximum-exponent
*))
169 (defparameter *intexp-maximum-exponent
* nil
)
171 ;;; This function precisely calculates base raised to an integral
172 ;;; power. It separates the cases by the sign of power, for efficiency
173 ;;; reasons, as powers can be calculated more efficiently if power is
174 ;;; a positive integer. Values of power are calculated as positive
175 ;;; integers, and inverted if negative.
176 (defun intexp (base power
)
177 (when (and *intexp-maximum-exponent
*
178 (> (abs power
) *intexp-maximum-exponent
*))
179 (error "The absolute value of ~S exceeds ~S."
180 power
'*intexp-maximum-exponent
*))
181 (cond ((minusp power
)
182 (/ (intexp base
(- power
))))
186 (do ((nextn (ash power -
1) (ash power -
1))
187 (total (if (oddp power
) base
1)
188 (if (oddp power
) (* base total
) total
)))
189 ((zerop nextn
) total
)
190 (setq base
(* base base
))
191 (setq power nextn
)))))
193 ;;; If an integer power of a rational, use INTEXP above. Otherwise, do
194 ;;; floating point stuff. If both args are real, we try %POW right
195 ;;; off, assuming it will return 0 if the result may be complex. If
196 ;;; so, we call COMPLEX-POW which directly computes the complex
197 ;;; result. We also separate the complex-real and real-complex cases
198 ;;; from the general complex case.
199 (defun expt (base power
)
201 "Return BASE raised to the POWER."
203 (if (and (zerop base
) (floatp power
))
204 (error 'arguments-out-of-domain-error
205 :operands
(list base power
)
207 :references
(list '(:ansi-cl
:function expt
)))
208 (let ((result (1+ (* base power
))))
209 (if (and (floatp result
) (float-nan-p result
))
212 (labels (;; determine if the double float is an integer.
213 ;; 0 - not an integer
217 (declare (type (unsigned-byte 31) ihi
)
218 (type (unsigned-byte 32) lo
)
219 (optimize (speed 3) (safety 0)))
221 (declare (type fixnum isint
))
222 (cond ((>= ihi
#x43400000
) ; exponent >= 53
225 (let ((k (- (ash ihi -
20) #x3ff
))) ; exponent
226 (declare (type (mod 53) k
))
228 (let* ((shift (- 52 k
))
229 (j (logand (ash lo
(- shift
))))
231 (declare (type (mod 32) shift
)
232 (type (unsigned-byte 32) j j2
))
234 (setq isint
(- 2 (logand j
1))))))
236 (let* ((shift (- 20 k
))
237 (j (ash ihi
(- shift
)))
239 (declare (type (mod 32) shift
)
240 (type (unsigned-byte 31) j j2
))
242 (setq isint
(- 2 (logand j
1))))))))))
244 (real-expt (x y rtype
)
245 (let ((x (coerce x
'double-float
))
246 (y (coerce y
'double-float
)))
247 (declare (double-float x y
))
248 (let* ((x-hi (double-float-high-bits x
))
249 (x-lo (double-float-low-bits x
))
250 (x-ihi (logand x-hi
#x7fffffff
))
251 (y-hi (double-float-high-bits y
))
252 (y-lo (double-float-low-bits y
))
253 (y-ihi (logand y-hi
#x7fffffff
)))
254 (declare (type (signed-byte 32) x-hi y-hi
)
255 (type (unsigned-byte 31) x-ihi y-ihi
)
256 (type (unsigned-byte 32) x-lo y-lo
))
258 (when (zerop (logior y-ihi y-lo
))
259 (return-from real-expt
(coerce 1d0 rtype
)))
261 ;; FIXME: Hardcoded qNaN/sNaN values are not portable.
262 (when (or (> x-ihi
#x7ff00000
)
263 (and (= x-ihi
#x7ff00000
) (/= x-lo
0))
265 (and (= y-ihi
#x7ff00000
) (/= y-lo
0)))
266 (return-from real-expt
(coerce (+ x y
) rtype
)))
267 (let ((yisint (if (< x-hi
0) (isint y-ihi y-lo
) 0)))
268 (declare (type fixnum yisint
))
269 ;; special value of y
270 (when (and (zerop y-lo
) (= y-ihi
#x7ff00000
))
272 (return-from real-expt
273 (cond ((and (= x-ihi
#x3ff00000
) (zerop x-lo
))
275 (coerce (- y y
) rtype
))
276 ((>= x-ihi
#x3ff00000
)
277 ;; (|x|>1)**+-inf = inf,0
282 ;; (|x|<1)**-,+inf = inf,0
285 (coerce 0 rtype
))))))
287 (let ((abs-x (abs x
)))
288 (declare (double-float abs-x
))
289 ;; special value of x
290 (when (and (zerop x-lo
)
291 (or (= x-ihi
#x7ff00000
) (zerop x-ihi
)
292 (= x-ihi
#x3ff00000
)))
293 ;; x is +-0,+-inf,+-1
294 (let ((z (if (< y-hi
0)
295 (/ 1 abs-x
) ; z = (1/|x|)
297 (declare (double-float z
))
299 (cond ((and (= x-ihi
#x3ff00000
) (zerop yisint
))
301 (let ((y*pi
(* y pi
)))
302 (declare (double-float y
*pi
))
303 (return-from real-expt
305 (coerce (%cos y
*pi
) rtype
)
306 (coerce (%sin y
*pi
) rtype
)))))
308 ;; (x<0)**odd = -(|x|**odd)
310 (return-from real-expt
(coerce z rtype
))))
314 (coerce (%pow x y
) rtype
)
316 (let ((pow (%pow abs-x y
)))
317 (declare (double-float pow
))
320 (coerce (* -
1d0 pow
) rtype
))
324 (let ((y*pi
(* y pi
)))
325 (declare (double-float y
*pi
))
327 (coerce (* pow
(%cos y
*pi
))
329 (coerce (* pow
(%sin y
*pi
))
331 (complex-expt (base power
)
332 (if (and (zerop base
) (plusp (realpart power
)))
334 (exp (* power
(log base
))))))
335 (declare (inline real-expt complex-expt
))
336 (number-dispatch ((base number
) (power number
))
337 (((foreach fixnum
(or bignum ratio
) (complex rational
)) integer
)
339 (((foreach single-float double-float
) rational
)
340 (real-expt base power
'(dispatch-type base
)))
341 (((foreach fixnum
(or bignum ratio
) single-float
)
342 (foreach ratio single-float
))
343 (real-expt base power
'single-float
))
344 (((foreach fixnum
(or bignum ratio
) single-float double-float
)
346 (real-expt base power
'double-float
))
347 ((double-float single-float
)
348 (real-expt base power
'double-float
))
349 ;; Handle (expt <complex> <rational>), except the case dealt with
350 ;; in the first clause above, (expt <(complex rational)> <integer>).
351 (((foreach (complex rational
) (complex single-float
)
352 (complex double-float
))
354 (* (expt (abs base
) power
)
355 (cis (* power
(phase base
)))))
356 ;; The next three clauses handle (expt <real> <complex>).
357 (((foreach fixnum
(or bignum ratio
) single-float
)
358 (foreach (complex single-float
) (complex rational
)))
359 (complex-expt base power
))
360 (((foreach fixnum
(or bignum ratio
) single-float
)
361 (complex double-float
))
362 (complex-expt (coerce base
'double-float
) power
))
363 ((double-float complex
)
364 (complex-expt base power
))
365 ;; The next three clauses handle (expt <complex> <float>) and
366 ;; (expt <complex> <complex>).
367 (((foreach (complex single-float
) (complex rational
))
368 (foreach (complex single-float
) (complex rational
) single-float
))
369 (complex-expt base power
))
370 (((foreach (complex single-float
) (complex rational
))
371 (foreach (complex double-float
) double-float
))
372 (complex-expt (coerce base
'(complex double-float
)) power
))
373 (((complex double-float
)
374 (foreach complex double-float single-float
))
375 (complex-expt base power
))))))
377 ;;; FIXME: Maybe rename this so that it's clearer that it only works
380 (declare (type integer x
))
383 ;; Write x = 2^n*f where 1/2 < f <= 1. Then log2(x) = n +
384 ;; log2(f). So we grab the top few bits of x and scale that
385 ;; appropriately, take the log of it and add it to n.
387 ;; Motivated by an attempt to get LOG to work better on bignums.
388 (let ((n (integer-length x
)))
389 (if (< n sb
!vm
:double-float-digits
)
390 (log (coerce x
'double-float
) 2.0d0
)
391 (let ((f (ldb (byte sb
!vm
:double-float-digits
392 (- n sb
!vm
:double-float-digits
))
394 (+ n
(log (scale-float (coerce f
'double-float
)
395 (- sb
!vm
:double-float-digits
))
398 (defun log (number &optional
(base nil base-p
))
400 "Return the logarithm of NUMBER in the base BASE, which defaults to e."
404 (if (or (typep number
'double-float
) (typep base
'double-float
))
407 ((and (typep number
'(integer (0) *))
408 (typep base
'(integer (0) *)))
409 (coerce (/ (log2 number
) (log2 base
)) 'single-float
))
410 ((and (typep number
'integer
) (typep base
'double-float
))
411 ;; No single float intermediate result
412 (/ (log2 number
) (log base
2.0d0
)))
413 ((and (typep number
'double-float
) (typep base
'integer
))
414 (/ (log number
2.0d0
) (log2 base
)))
416 (/ (log number
) (log base
))))
417 (number-dispatch ((number number
))
418 (((foreach fixnum bignum
))
420 (complex (log (- number
)) (coerce pi
'single-float
))
421 (coerce (/ (log2 number
) (log (exp 1.0d0
) 2.0d0
)) 'single-float
)))
424 (complex (log (- number
)) (coerce pi
'single-float
))
425 (let ((numerator (numerator number
))
426 (denominator (denominator number
)))
427 (if (= (integer-length numerator
)
428 (integer-length denominator
))
429 (coerce (%log1p
(coerce (- number
1) 'double-float
))
431 (coerce (/ (- (log2 numerator
) (log2 denominator
))
432 (log (exp 1.0d0
) 2.0d0
))
434 (((foreach single-float double-float
))
435 ;; Is (log -0) -infinity (libm.a) or -infinity + i*pi (Kahan)?
436 ;; Since this doesn't seem to be an implementation issue
437 ;; I (pw) take the Kahan result.
438 (if (< (float-sign number
)
439 (coerce 0 '(dispatch-type number
)))
440 (complex (log (- number
)) (coerce pi
'(dispatch-type number
)))
441 (coerce (%log
(coerce number
'double-float
))
442 '(dispatch-type number
))))
444 (complex-log number
)))))
448 "Return the square root of NUMBER."
449 (number-dispatch ((number number
))
450 (((foreach fixnum bignum ratio
))
452 (complex-sqrt number
)
453 (coerce (%sqrt
(coerce number
'double-float
)) 'single-float
)))
454 (((foreach single-float double-float
))
456 (complex-sqrt (complex number
))
457 (coerce (%sqrt
(coerce number
'double-float
))
458 '(dispatch-type number
))))
460 (complex-sqrt number
))))
462 ;;;; trigonometic and related functions
466 "Return the absolute value of the number."
467 (number-dispatch ((number number
))
468 (((foreach single-float double-float fixnum rational
))
471 (let ((rx (realpart number
))
472 (ix (imagpart number
)))
475 (sqrt (+ (* rx rx
) (* ix ix
))))
477 (coerce (%hypot
(coerce rx
'double-float
)
478 (coerce (truly-the single-float ix
) 'double-float
))
481 (%hypot rx
(truly-the double-float ix
))))))))
483 (defun phase (number)
485 "Return the angle part of the polar representation of a complex number.
486 For complex numbers, this is (atan (imagpart number) (realpart number)).
487 For non-complex positive numbers, this is 0. For non-complex negative
492 (coerce pi
'single-float
)
495 (if (minusp (float-sign number
))
496 (coerce pi
'single-float
)
499 (if (minusp (float-sign number
))
500 (coerce pi
'double-float
)
503 (atan (imagpart number
) (realpart number
)))))
507 "Return the sine of NUMBER."
508 (number-dispatch ((number number
))
509 (handle-reals %sin number
)
511 (let ((x (realpart number
))
512 (y (imagpart number
)))
513 (complex (* (sin x
) (cosh y
))
514 (* (cos x
) (sinh y
)))))))
518 "Return the cosine of NUMBER."
519 (number-dispatch ((number number
))
520 (handle-reals %cos number
)
522 (let ((x (realpart number
))
523 (y (imagpart number
)))
524 (complex (* (cos x
) (cosh y
))
525 (- (* (sin x
) (sinh y
))))))))
529 "Return the tangent of NUMBER."
530 (number-dispatch ((number number
))
531 (handle-reals %tan number
)
533 (complex-tan number
))))
537 "Return cos(Theta) + i sin(Theta), i.e. exp(i Theta)."
538 (declare (type real theta
))
539 (complex (cos theta
) (sin theta
)))
543 "Return the arc sine of NUMBER."
544 (number-dispatch ((number number
))
546 (if (or (> number
1) (< number -
1))
547 (complex-asin number
)
548 (coerce (%asin
(coerce number
'double-float
)) 'single-float
)))
549 (((foreach single-float double-float
))
550 (if (or (> number
(coerce 1 '(dispatch-type number
)))
551 (< number
(coerce -
1 '(dispatch-type number
))))
552 (complex-asin (complex number
))
553 (coerce (%asin
(coerce number
'double-float
))
554 '(dispatch-type number
))))
556 (complex-asin number
))))
560 "Return the arc cosine of NUMBER."
561 (number-dispatch ((number number
))
563 (if (or (> number
1) (< number -
1))
564 (complex-acos number
)
565 (coerce (%acos
(coerce number
'double-float
)) 'single-float
)))
566 (((foreach single-float double-float
))
567 (if (or (> number
(coerce 1 '(dispatch-type number
)))
568 (< number
(coerce -
1 '(dispatch-type number
))))
569 (complex-acos (complex number
))
570 (coerce (%acos
(coerce number
'double-float
))
571 '(dispatch-type number
))))
573 (complex-acos number
))))
575 (defun atan (y &optional
(x nil xp
))
577 "Return the arc tangent of Y if X is omitted or Y/X if X is supplied."
580 (declare (type double-float y x
)
581 (values double-float
))
584 (if (plusp (float-sign x
))
587 (float-sign y
(/ pi
2)))
589 (number-dispatch ((y real
) (x real
))
591 (foreach double-float single-float fixnum bignum ratio
))
592 (atan2 y
(coerce x
'double-float
)))
593 (((foreach single-float fixnum bignum ratio
)
595 (atan2 (coerce y
'double-float
) x
))
596 (((foreach single-float fixnum bignum ratio
)
597 (foreach single-float fixnum bignum ratio
))
598 (coerce (atan2 (coerce y
'double-float
) (coerce x
'double-float
))
600 (number-dispatch ((y number
))
601 (handle-reals %atan y
)
605 ;;; It seems that every target system has a C version of sinh, cosh,
606 ;;; and tanh. Let's use these for reals because the original
607 ;;; implementations based on the definitions lose big in round-off
608 ;;; error. These bad definitions also mean that sin and cos for
609 ;;; complex numbers can also lose big.
613 "Return the hyperbolic sine of NUMBER."
614 (number-dispatch ((number number
))
615 (handle-reals %sinh number
)
617 (let ((x (realpart number
))
618 (y (imagpart number
)))
619 (complex (* (sinh x
) (cos y
))
620 (* (cosh x
) (sin y
)))))))
624 "Return the hyperbolic cosine of NUMBER."
625 (number-dispatch ((number number
))
626 (handle-reals %cosh number
)
628 (let ((x (realpart number
))
629 (y (imagpart number
)))
630 (complex (* (cosh x
) (cos y
))
631 (* (sinh x
) (sin y
)))))))
635 "Return the hyperbolic tangent of NUMBER."
636 (number-dispatch ((number number
))
637 (handle-reals %tanh number
)
639 (complex-tanh number
))))
641 (defun asinh (number)
643 "Return the hyperbolic arc sine of NUMBER."
644 (number-dispatch ((number number
))
645 (handle-reals %asinh number
)
647 (complex-asinh number
))))
649 (defun acosh (number)
651 "Return the hyperbolic arc cosine of NUMBER."
652 (number-dispatch ((number number
))
654 ;; acosh is complex if number < 1
656 (complex-acosh number
)
657 (coerce (%acosh
(coerce number
'double-float
)) 'single-float
)))
658 (((foreach single-float double-float
))
659 (if (< number
(coerce 1 '(dispatch-type number
)))
660 (complex-acosh (complex number
))
661 (coerce (%acosh
(coerce number
'double-float
))
662 '(dispatch-type number
))))
664 (complex-acosh number
))))
666 (defun atanh (number)
668 "Return the hyperbolic arc tangent of NUMBER."
669 (number-dispatch ((number number
))
671 ;; atanh is complex if |number| > 1
672 (if (or (> number
1) (< number -
1))
673 (complex-atanh number
)
674 (coerce (%atanh
(coerce number
'double-float
)) 'single-float
)))
675 (((foreach single-float double-float
))
676 (if (or (> number
(coerce 1 '(dispatch-type number
)))
677 (< number
(coerce -
1 '(dispatch-type number
))))
678 (complex-atanh (complex number
))
679 (coerce (%atanh
(coerce number
'double-float
))
680 '(dispatch-type number
))))
682 (complex-atanh number
))))
685 ;;;; not-OLD-SPECFUN stuff
687 ;;;; (This was conditional on #-OLD-SPECFUN in the CMU CL sources,
688 ;;;; but OLD-SPECFUN was mentioned nowhere else, so it seems to be
689 ;;;; the standard special function system.)
691 ;;;; This is a set of routines that implement many elementary
692 ;;;; transcendental functions as specified by ANSI Common Lisp. The
693 ;;;; implementation is based on Kahan's paper.
695 ;;;; I believe I have accurately implemented the routines and are
696 ;;;; correct, but you may want to check for your self.
698 ;;;; These functions are written for CMU Lisp and take advantage of
699 ;;;; some of the features available there. It may be possible,
700 ;;;; however, to port this to other Lisps.
702 ;;;; Some functions are significantly more accurate than the original
703 ;;;; definitions in CMU Lisp. In fact, some functions in CMU Lisp
704 ;;;; give the wrong answer like (acos #c(-2.0 0.0)), where the true
705 ;;;; answer is pi + i*log(2-sqrt(3)).
707 ;;;; All of the implemented functions will take any number for an
708 ;;;; input, but the result will always be a either a complex
709 ;;;; single-float or a complex double-float.
711 ;;;; general functions:
723 ;;;; utility functions:
726 ;;;; internal functions:
727 ;;;; square coerce-to-complex-type cssqs complex-log-scaled
730 ;;;; Kahan, W. "Branch Cuts for Complex Elementary Functions, or Much
731 ;;;; Ado About Nothing's Sign Bit" in Iserles and Powell (eds.) "The
732 ;;;; State of the Art in Numerical Analysis", pp. 165-211, Clarendon
735 ;;;; The original CMU CL code requested:
736 ;;;; Please send any bug reports, comments, or improvements to
737 ;;;; Raymond Toy at <email address deleted during 2002 spam avalanche>.
739 ;;; FIXME: In SBCL, the floating point infinity constants like
740 ;;; SB!EXT:DOUBLE-FLOAT-POSITIVE-INFINITY aren't available as
741 ;;; constants at cross-compile time, because the cross-compilation
742 ;;; host might not have support for floating point infinities. Thus,
743 ;;; they're effectively implemented as special variable references,
744 ;;; and the code below which uses them might be unnecessarily
745 ;;; inefficient. Perhaps some sort of MAKE-LOAD-TIME-VALUE hackery
746 ;;; should be used instead? (KLUDGED 2004-03-08 CSR, by replacing the
747 ;;; special variable references with (probably equally slow)
750 ;;; FIXME: As of 2004-05, when PFD noted that IMAGPART and COMPLEX
751 ;;; differ in their interpretations of the real line, IMAGPART was
752 ;;; patch, which without a certain amount of effort would have altered
753 ;;; all the branch cut treatment. Clients of these COMPLEX- routines
754 ;;; were patched to use explicit COMPLEX, rather than implicitly
755 ;;; passing in real numbers for treatment with IMAGPART, and these
756 ;;; COMPLEX- functions altered to require arguments of type COMPLEX;
757 ;;; however, someone needs to go back to Kahan for the definitive
758 ;;; answer for treatment of negative real floating point numbers and
759 ;;; branch cuts. If adjustment is needed, it is probably the removal
760 ;;; of explicit calls to COMPLEX in the clients of irrational
761 ;;; functions. -- a slightly bitter CSR, 2004-05-16
763 (declaim (inline square
))
765 (declare (double-float x
))
768 ;;; original CMU CL comment, apparently re. SCALB and LOGB and
770 ;;; If you have these functions in libm, perhaps they should be used
771 ;;; instead of these Lisp versions. These versions are probably good
772 ;;; enough, especially since they are portable.
774 ;;; Compute 2^N * X without computing 2^N first. (Use properties of
775 ;;; the underlying floating-point format.)
776 (declaim (inline scalb
))
778 (declare (type double-float x
)
779 (type double-float-exponent n
))
782 ;;; This is like LOGB, but X is not infinity and non-zero and not a
783 ;;; NaN, so we can always return an integer.
784 (declaim (inline logb-finite
))
785 (defun logb-finite (x)
786 (declare (type double-float x
))
787 (multiple-value-bind (signif exponent sign
)
789 (declare (ignore signif sign
))
790 ;; DECODE-FLOAT is almost right, except that the exponent is off
794 ;;; Compute an integer N such that 1 <= |2^N * x| < 2.
795 ;;; For the special cases, the following values are used:
798 ;;; +/- infinity +infinity
801 (declare (type double-float x
))
802 (cond ((float-nan-p x
)
804 ((float-infinity-p x
)
805 ;; DOUBLE-FLOAT-POSITIVE-INFINITY
806 (double-from-bits 0 (1+ sb
!vm
:double-float-normal-exponent-max
) 0))
808 ;; The answer is negative infinity, but we are supposed to
809 ;; signal divide-by-zero, so do the actual division
815 ;;; This function is used to create a complex number of the
816 ;;; appropriate type:
817 ;;; Create complex number with real part X and imaginary part Y
818 ;;; such that has the same type as Z. If Z has type (complex
819 ;;; rational), the X and Y are coerced to single-float.
820 #!+long-float
(eval-when (:compile-toplevel
:load-toplevel
:execute
)
821 (error "needs work for long float support"))
822 (declaim (inline coerce-to-complex-type
))
823 (defun coerce-to-complex-type (x y z
)
824 (declare (double-float x y
)
826 (if (typep (realpart z
) 'double-float
)
828 ;; Convert anything that's not already a DOUBLE-FLOAT (because
829 ;; the initial argument was a (COMPLEX DOUBLE-FLOAT) and we
830 ;; haven't done anything to lose precision) to a SINGLE-FLOAT.
831 (complex (float x
1f0
)
834 ;;; Compute |(x+i*y)/2^k|^2 scaled to avoid over/underflow. The
835 ;;; result is r + i*k, where k is an integer.
836 #!+long-float
(eval-when (:compile-toplevel
:load-toplevel
:execute
)
837 (error "needs work for long float support"))
839 (let ((x (float (realpart z
) 1d0
))
840 (y (float (imagpart z
) 1d0
)))
841 ;; Would this be better handled using an exception handler to
842 ;; catch the overflow or underflow signal? For now, we turn all
843 ;; traps off and look at the accrued exceptions to see if any
844 ;; signal would have been raised.
845 (with-float-traps-masked (:underflow
:overflow
)
846 (let ((rho (+ (square x
) (square y
))))
847 (declare (optimize (speed 3) (space 0)))
848 (cond ((and (or (float-nan-p rho
)
849 (float-infinity-p rho
))
850 (or (float-infinity-p (abs x
))
851 (float-infinity-p (abs y
))))
852 ;; DOUBLE-FLOAT-POSITIVE-INFINITY
854 (double-from-bits 0 (1+ sb
!vm
:double-float-normal-exponent-max
) 0)
857 ;; (/ least-positive-double-float double-float-epsilon)
860 (make-double-float #x1fffff
#xfffffffe
)
862 (error "(/ least-positive-long-float long-float-epsilon)")))
863 (traps (ldb sb
!vm
::float-sticky-bits
864 (sb!vm
:floating-point-modes
))))
865 ;; Overflow raised or (underflow raised and rho <
867 (or (not (zerop (logand sb
!vm
:float-overflow-trap-bit traps
)))
868 (and (not (zerop (logand sb
!vm
:float-underflow-trap-bit
871 ;; If we're here, neither x nor y are infinity and at
872 ;; least one is non-zero.. Thus logb returns a nice
874 (let ((k (- (logb-finite (max (abs x
) (abs y
))))))
875 (values (+ (square (scalb x k
))
876 (square (scalb y k
)))
881 ;;; principal square root of Z
883 ;;; Z may be RATIONAL or COMPLEX; the result is always a COMPLEX.
884 (defun complex-sqrt (z)
885 ;; KLUDGE: Here and below, we can't just declare Z to be of type
886 ;; COMPLEX, because one-arg COMPLEX on rationals returns a rational.
887 ;; Since there isn't a rational negative zero, this is OK from the
888 ;; point of view of getting the right answer in the face of branch
889 ;; cuts, but declarations of the form (OR RATIONAL COMPLEX) are
890 ;; still ugly. -- CSR, 2004-05-16
891 (declare (type (or complex rational
) z
))
892 (multiple-value-bind (rho k
)
894 (declare (type (or (member 0d0
) (double-float 0d0
)) rho
)
896 (let ((x (float (realpart z
) 1.0d0
))
897 (y (float (imagpart z
) 1.0d0
))
900 (declare (double-float x y eta nu
))
903 ;; space 0 to get maybe-inline functions inlined.
904 (declare (optimize (speed 3) (space 0)))
906 (if (not (float-nan-p x
))
907 (setf rho
(+ (scalb (abs x
) (- k
)) (sqrt rho
))))
912 (setf k
(1- (ash k -
1)))
913 (setf rho
(+ rho rho
))))
915 (setf rho
(scalb (sqrt rho
) k
))
921 (when (not (float-infinity-p (abs nu
)))
922 (setf nu
(/ (/ nu rho
) 2d0
)))
925 (setf nu
(float-sign y rho
))))
926 (coerce-to-complex-type eta nu z
)))))
928 ;;; Compute log(2^j*z).
930 ;;; This is for use with J /= 0 only when |z| is huge.
931 (defun complex-log-scaled (z j
)
932 (declare (type (or rational complex
) z
)
934 ;; The constants t0, t1, t2 should be evaluated to machine
935 ;; precision. In addition, Kahan says the accuracy of log1p
936 ;; influences the choices of these constants but doesn't say how to
937 ;; choose them. We'll just assume his choices matches our
938 ;; implementation of log1p.
939 (let ((t0 (load-time-value
941 (make-double-float #x3fe6a09e
#x667f3bcd
)
943 (error "(/ (sqrt 2l0))")))
944 ;; KLUDGE: if repeatable fasls start failing under some weird
945 ;; xc host, this 1.2d0 might be a good place to examine: while
946 ;; it _should_ be the same in all vaguely-IEEE754 hosts, 1.2
947 ;; is not exactly representable, so something could go wrong.
950 (ln2 (load-time-value
952 (make-double-float #x3fe62e42
#xfefa39ef
)
954 (error "(log 2l0)")))
955 (x (float (realpart z
) 1.0d0
))
956 (y (float (imagpart z
) 1.0d0
)))
957 (multiple-value-bind (rho k
)
959 (declare (optimize (speed 3)))
960 (let ((beta (max (abs x
) (abs y
)))
961 (theta (min (abs x
) (abs y
))))
962 (coerce-to-complex-type (if (and (zerop k
)
966 (/ (%log1p
(+ (* (- beta
1.0d0
)
975 ;;; log of Z = log |Z| + i * arg Z
977 ;;; Z may be any number, but the result is always a complex.
978 (defun complex-log (z)
979 (declare (type (or rational complex
) z
))
980 (complex-log-scaled z
0))
982 ;;; KLUDGE: Let us note the following "strange" behavior. atanh 1.0d0
983 ;;; is +infinity, but the following code returns approx 176 + i*pi/4.
984 ;;; The reason for the imaginary part is caused by the fact that arg
985 ;;; i*y is never 0 since we have positive and negative zeroes. -- rtoy
986 ;;; Compute atanh z = (log(1+z) - log(1-z))/2.
987 (defun complex-atanh (z)
988 (declare (type (or rational complex
) z
))
990 (theta (/ (sqrt most-positive-double-float
) 4.0d0
))
991 (rho (/ 4.0d0
(sqrt most-positive-double-float
)))
992 (half-pi (/ pi
2.0d0
))
993 (rp (float (realpart z
) 1.0d0
))
994 (beta (float-sign rp
1.0d0
))
996 (y (* beta
(- (float (imagpart z
) 1.0d0
))))
999 ;; Shouldn't need this declare.
1000 (declare (double-float x y
))
1002 (declare (optimize (speed 3)))
1003 (cond ((or (> x theta
)
1005 ;; To avoid overflow...
1006 (setf nu
(float-sign y half-pi
))
1007 ;; ETA is real part of 1/(x + iy). This is x/(x^2+y^2),
1008 ;; which can cause overflow. Arrange this computation so
1009 ;; that it won't overflow.
1010 (setf eta
(let* ((x-bigger (> x
(abs y
)))
1011 (r (if x-bigger
(/ y x
) (/ x y
)))
1012 (d (+ 1.0d0
(* r r
))))
1017 ;; Should this be changed so that if y is zero, eta is set
1018 ;; to +infinity instead of approx 176? In any case
1019 ;; tanh(176) is 1.0d0 within working precision.
1020 (let ((t1 (+ 4d0
(square y
)))
1021 (t2 (+ (abs y
) rho
)))
1022 (setf eta
(log (/ (sqrt (sqrt t1
))
1026 (+ half-pi
(atan (* 0.5d0 t2
))))))))
1028 (let ((t1 (+ (abs y
) rho
)))
1029 ;; Normal case using log1p(x) = log(1 + x)
1031 (%log1p
(/ (* 4.0d0 x
)
1032 (+ (square (- 1.0d0 x
))
1039 (coerce-to-complex-type (* beta eta
)
1043 ;;; Compute tanh z = sinh z / cosh z.
1044 (defun complex-tanh (z)
1045 (declare (type (or rational complex
) z
))
1046 (let ((x (float (realpart z
) 1.0d0
))
1047 (y (float (imagpart z
) 1.0d0
)))
1049 ;; space 0 to get maybe-inline functions inlined
1050 (declare (optimize (speed 3) (space 0)))
1054 (make-double-float #x406633ce
#x8fb9f87e
)
1056 (error "(/ (+ (log 2l0) (log most-positive-long-float)) 4l0)")))
1057 (coerce-to-complex-type (float-sign x
)
1060 (let* ((tv (%tan y
))
1061 (beta (+ 1.0d0
(* tv tv
)))
1063 (rho (sqrt (+ 1.0d0
(* s s
)))))
1064 (if (float-infinity-p (abs tv
))
1065 (coerce-to-complex-type (/ rho s
)
1068 (let ((den (+ 1.0d0
(* beta s s
))))
1069 (coerce-to-complex-type (/ (* beta rho s
)
1074 ;;; Compute acos z = pi/2 - asin z.
1076 ;;; Z may be any NUMBER, but the result is always a COMPLEX.
1077 (defun complex-acos (z)
1078 ;; Kahan says we should only compute the parts needed. Thus, the
1079 ;; REALPART's below should only compute the real part, not the whole
1080 ;; complex expression. Doing this can be important because we may get
1081 ;; spurious signals that occur in the part that we are not using.
1083 ;; However, we take a pragmatic approach and just use the whole
1086 ;; NOTE: The formula given by Kahan is somewhat ambiguous in whether
1087 ;; it's the conjugate of the square root or the square root of the
1088 ;; conjugate. This needs to be checked.
1090 ;; I checked. It doesn't matter because (conjugate (sqrt z)) is the
1091 ;; same as (sqrt (conjugate z)) for all z. This follows because
1093 ;; (conjugate (sqrt z)) = exp(0.5*log |z|)*exp(-0.5*j*arg z).
1095 ;; (sqrt (conjugate z)) = exp(0.5*log|z|)*exp(0.5*j*arg conj z)
1097 ;; and these two expressions are equal if and only if arg conj z =
1098 ;; -arg z, which is clearly true for all z.
1099 (declare (type (or rational complex
) z
))
1100 (let ((sqrt-1+z
(complex-sqrt (+ 1 z
)))
1101 (sqrt-1-z (complex-sqrt (- 1 z
))))
1102 (with-float-traps-masked (:divide-by-zero
)
1103 (complex (* 2 (atan (/ (realpart sqrt-1-z
)
1104 (realpart sqrt-1
+z
))))
1105 (asinh (imagpart (* (conjugate sqrt-1
+z
)
1108 ;;; Compute acosh z = 2 * log(sqrt((z+1)/2) + sqrt((z-1)/2))
1110 ;;; Z may be any NUMBER, but the result is always a COMPLEX.
1111 (defun complex-acosh (z)
1112 (declare (type (or rational complex
) z
))
1113 (let ((sqrt-z-1 (complex-sqrt (- z
1)))
1114 (sqrt-z+1 (complex-sqrt (+ z
1))))
1115 (with-float-traps-masked (:divide-by-zero
)
1116 (complex (asinh (realpart (* (conjugate sqrt-z-1
)
1118 (* 2 (atan (/ (imagpart sqrt-z-1
)
1119 (realpart sqrt-z
+1))))))))
1121 ;;; Compute asin z = asinh(i*z)/i.
1123 ;;; Z may be any NUMBER, but the result is always a COMPLEX.
1124 (defun complex-asin (z)
1125 (declare (type (or rational complex
) z
))
1126 (let ((sqrt-1-z (complex-sqrt (- 1 z
)))
1127 (sqrt-1+z
(complex-sqrt (+ 1 z
))))
1128 (with-float-traps-masked (:divide-by-zero
)
1129 (complex (atan (/ (realpart z
)
1130 (realpart (* sqrt-1-z sqrt-1
+z
))))
1131 (asinh (imagpart (* (conjugate sqrt-1-z
)
1134 ;;; Compute asinh z = log(z + sqrt(1 + z*z)).
1136 ;;; Z may be any number, but the result is always a complex.
1137 (defun complex-asinh (z)
1138 (declare (type (or rational complex
) z
))
1139 ;; asinh z = -i * asin (i*z)
1140 (let* ((iz (complex (- (imagpart z
)) (realpart z
)))
1141 (result (complex-asin iz
)))
1142 (complex (imagpart result
)
1143 (- (realpart result
)))))
1145 ;;; Compute atan z = atanh (i*z) / i.
1147 ;;; Z may be any number, but the result is always a complex.
1148 (defun complex-atan (z)
1149 (declare (type (or rational complex
) z
))
1150 ;; atan z = -i * atanh (i*z)
1151 (let* ((iz (complex (- (imagpart z
)) (realpart z
)))
1152 (result (complex-atanh iz
)))
1153 (complex (imagpart result
)
1154 (- (realpart result
)))))
1156 ;;; Compute tan z = -i * tanh(i * z)
1158 ;;; Z may be any number, but the result is always a complex.
1159 (defun complex-tan (z)
1160 (declare (type (or rational complex
) z
))
1161 ;; tan z = -i * tanh(i*z)
1162 (let* ((iz (complex (- (imagpart z
)) (realpart z
)))
1163 (result (complex-tanh iz
)))
1164 (complex (imagpart result
)
1165 (- (realpart result
)))))