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
))))
28 (declaim (inline ,function
))
29 (sb!alien
:define-alien-routine
(,name
,function
) double-float
30 ,@(let ((results nil
))
31 (dotimes (i num-args
(nreverse results
))
32 (push (list (intern (format nil
"ARG-~D" i
))
36 (defun handle-reals (function var
)
37 `((((foreach fixnum single-float bignum ratio
))
38 (coerce (,function
(coerce ,var
'double-float
)) 'single-float
))
44 #!+x86
;; for constant folding
45 (macrolet ((def (name ll
)
46 `(defun ,name
,ll
(,name
,@ll
))))
59 #!+x86-64
;; for constant folding
60 (macrolet ((def (name ll
)
61 `(defun ,name
,ll
(,name
,@ll
))))
64 ;;;; stubs for the Unix math library
66 ;;;; Many of these are unnecessary on the X86 because they're built
70 #!-x86
(def-math-rtn "sin" 1)
71 #!-x86
(def-math-rtn "cos" 1)
72 #!-x86
(def-math-rtn "tan" 1)
73 (def-math-rtn "asin" 1)
74 (def-math-rtn "acos" 1)
75 #!-x86
(def-math-rtn "atan" 1)
76 #!-x86
(def-math-rtn "atan2" 2)
77 (def-math-rtn "sinh" 1)
78 (def-math-rtn "cosh" 1)
81 (def-math-rtn "tanh" 1)
82 (def-math-rtn "asinh" 1)
83 (def-math-rtn "acosh" 1)
84 (def-math-rtn "atanh" 1))
87 (declaim (inline %tanh
))
89 (/ (%sinh number
) (%cosh number
)))
90 (declaim (inline %asinh
))
91 (defun %asinh
(number)
92 (log (+ number
(sqrt (+ (* number number
) 1.0d0
))) #.
(exp 1.0d0
)))
93 (declaim (inline %acosh
))
94 (defun %acosh
(number)
95 (log (+ number
(sqrt (- (* number number
) 1.0d0
))) #.
(exp 1.0d0
)))
96 (declaim (inline %atanh
))
97 (defun %atanh
(number)
98 (let ((ratio (/ (+ 1 number
) (- 1 number
))))
99 ;; Were we effectively zero?
102 (/ (log ratio
#.
(exp 1.0d0
)) 2.0d0
)))))
104 ;;; exponential and logarithmic
105 #!-x86
(def-math-rtn "exp" 1)
106 #!-x86
(def-math-rtn "log" 1)
107 #!-x86
(def-math-rtn "log10" 1)
108 #!-win32
(def-math-rtn "pow" 2)
109 #!-
(or x86 x86-64
) (def-math-rtn "sqrt" 1)
110 (def-math-rtn "hypot" 2)
111 #!-
(or hpux x86
) (def-math-rtn "log1p" 1)
117 "Return e raised to the power NUMBER."
118 (number-dispatch ((number number
))
119 (handle-reals %exp number
)
121 (* (exp (realpart number
))
122 (cis (imagpart number
))))))
124 ;;; INTEXP -- Handle the rational base, integer power case.
126 (declaim (type (or integer null
) *intexp-maximum-exponent
*))
127 (defparameter *intexp-maximum-exponent
* nil
)
129 ;;; This function precisely calculates base raised to an integral
130 ;;; power. It separates the cases by the sign of power, for efficiency
131 ;;; reasons, as powers can be calculated more efficiently if power is
132 ;;; a positive integer. Values of power are calculated as positive
133 ;;; integers, and inverted if negative.
134 (defun intexp (base power
)
135 (when (and *intexp-maximum-exponent
*
136 (> (abs power
) *intexp-maximum-exponent
*))
137 (error "The absolute value of ~S exceeds ~S."
138 power
'*intexp-maximum-exponent
*))
139 (cond ((minusp power
)
140 (/ (intexp base
(- power
))))
144 (do ((nextn (ash power -
1) (ash power -
1))
145 (total (if (oddp power
) base
1)
146 (if (oddp power
) (* base total
) total
)))
147 ((zerop nextn
) total
)
148 (setq base
(* base base
))
149 (setq power nextn
)))))
151 ;;; If an integer power of a rational, use INTEXP above. Otherwise, do
152 ;;; floating point stuff. If both args are real, we try %POW right
153 ;;; off, assuming it will return 0 if the result may be complex. If
154 ;;; so, we call COMPLEX-POW which directly computes the complex
155 ;;; result. We also separate the complex-real and real-complex cases
156 ;;; from the general complex case.
157 (defun expt (base power
)
159 "Return BASE raised to the POWER."
161 (let ((result (1+ (* base power
))))
162 (if (and (floatp result
) (float-nan-p result
))
165 (labels (;; determine if the double float is an integer.
166 ;; 0 - not an integer
170 (declare (type (unsigned-byte 31) ihi
)
171 (type (unsigned-byte 32) lo
)
172 (optimize (speed 3) (safety 0)))
174 (declare (type fixnum isint
))
175 (cond ((>= ihi
#x43400000
) ; exponent >= 53
178 (let ((k (- (ash ihi -
20) #x3ff
))) ; exponent
179 (declare (type (mod 53) k
))
181 (let* ((shift (- 52 k
))
182 (j (logand (ash lo
(- shift
))))
184 (declare (type (mod 32) shift
)
185 (type (unsigned-byte 32) j j2
))
187 (setq isint
(- 2 (logand j
1))))))
189 (let* ((shift (- 20 k
))
190 (j (ash ihi
(- shift
)))
192 (declare (type (mod 32) shift
)
193 (type (unsigned-byte 31) j j2
))
195 (setq isint
(- 2 (logand j
1))))))))))
197 (real-expt (x y rtype
)
198 (let ((x (coerce x
'double-float
))
199 (y (coerce y
'double-float
)))
200 (declare (double-float x y
))
201 (let* ((x-hi (sb!kernel
:double-float-high-bits x
))
202 (x-lo (sb!kernel
:double-float-low-bits x
))
203 (x-ihi (logand x-hi
#x7fffffff
))
204 (y-hi (sb!kernel
:double-float-high-bits y
))
205 (y-lo (sb!kernel
:double-float-low-bits y
))
206 (y-ihi (logand y-hi
#x7fffffff
)))
207 (declare (type (signed-byte 32) x-hi y-hi
)
208 (type (unsigned-byte 31) x-ihi y-ihi
)
209 (type (unsigned-byte 32) x-lo y-lo
))
211 (when (zerop (logior y-ihi y-lo
))
212 (return-from real-expt
(coerce 1d0 rtype
)))
214 (when (or (> x-ihi
#x7ff00000
)
215 (and (= x-ihi
#x7ff00000
) (/= x-lo
0))
217 (and (= y-ihi
#x7ff00000
) (/= y-lo
0)))
218 (return-from real-expt
(coerce (+ x y
) rtype
)))
219 (let ((yisint (if (< x-hi
0) (isint y-ihi y-lo
) 0)))
220 (declare (type fixnum yisint
))
221 ;; special value of y
222 (when (and (zerop y-lo
) (= y-ihi
#x7ff00000
))
224 (return-from real-expt
225 (cond ((and (= x-ihi
#x3ff00000
) (zerop x-lo
))
227 (coerce (- y y
) rtype
))
228 ((>= x-ihi
#x3ff00000
)
229 ;; (|x|>1)**+-inf = inf,0
234 ;; (|x|<1)**-,+inf = inf,0
237 (coerce 0 rtype
))))))
239 (let ((abs-x (abs x
)))
240 (declare (double-float abs-x
))
241 ;; special value of x
242 (when (and (zerop x-lo
)
243 (or (= x-ihi
#x7ff00000
) (zerop x-ihi
)
244 (= x-ihi
#x3ff00000
)))
245 ;; x is +-0,+-inf,+-1
246 (let ((z (if (< y-hi
0)
247 (/ 1 abs-x
) ; z = (1/|x|)
249 (declare (double-float z
))
251 (cond ((and (= x-ihi
#x3ff00000
) (zerop yisint
))
253 (let ((y*pi
(* y pi
)))
254 (declare (double-float y
*pi
))
255 (return-from real-expt
257 (coerce (%cos y
*pi
) rtype
)
258 (coerce (%sin y
*pi
) rtype
)))))
260 ;; (x<0)**odd = -(|x|**odd)
262 (return-from real-expt
(coerce z rtype
))))
266 (coerce (sb!kernel
::%pow x y
) rtype
)
268 (let ((pow (sb!kernel
::%pow abs-x y
)))
269 (declare (double-float pow
))
272 (coerce (* -
1d0 pow
) rtype
))
276 (let ((y*pi
(* y pi
)))
277 (declare (double-float y
*pi
))
279 (coerce (* pow
(%cos y
*pi
))
281 (coerce (* pow
(%sin y
*pi
))
283 (declare (inline real-expt
))
284 (number-dispatch ((base number
) (power number
))
285 (((foreach fixnum
(or bignum ratio
) (complex rational
)) integer
)
287 (((foreach single-float double-float
) rational
)
288 (real-expt base power
'(dispatch-type base
)))
289 (((foreach fixnum
(or bignum ratio
) single-float
)
290 (foreach ratio single-float
))
291 (real-expt base power
'single-float
))
292 (((foreach fixnum
(or bignum ratio
) single-float double-float
)
294 (real-expt base power
'double-float
))
295 ((double-float single-float
)
296 (real-expt base power
'double-float
))
297 (((foreach (complex rational
) (complex float
)) rational
)
298 (* (expt (abs base
) power
)
299 (cis (* power
(phase base
)))))
300 (((foreach fixnum
(or bignum ratio
) single-float double-float
)
302 (if (and (zerop base
) (plusp (realpart power
)))
304 (exp (* power
(log base
)))))
305 (((foreach (complex float
) (complex rational
))
306 (foreach complex double-float single-float
))
307 (if (and (zerop base
) (plusp (realpart power
)))
309 (exp (* power
(log base
)))))))))
311 ;;; FIXME: Maybe rename this so that it's clearer that it only works
314 (declare (type integer x
))
317 ;; Write x = 2^n*f where 1/2 < f <= 1. Then log2(x) = n +
318 ;; log2(f). So we grab the top few bits of x and scale that
319 ;; appropriately, take the log of it and add it to n.
321 ;; Motivated by an attempt to get LOG to work better on bignums.
322 (let ((n (integer-length x
)))
323 (if (< n sb
!vm
:double-float-digits
)
324 (log (coerce x
'double-float
) 2.0d0
)
325 (let ((f (ldb (byte sb
!vm
:double-float-digits
326 (- n sb
!vm
:double-float-digits
))
328 (+ n
(log (scale-float (coerce f
'double-float
)
329 (- sb
!vm
:double-float-digits
))
332 (defun log (number &optional
(base nil base-p
))
334 "Return the logarithm of NUMBER in the base BASE, which defaults to e."
337 ((zerop base
) 0f0
) ; FIXME: type
338 ((and (typep number
'(integer (0) *))
339 (typep base
'(integer (0) *)))
340 (coerce (/ (log2 number
) (log2 base
)) 'single-float
))
341 (t (/ (log number
) (log base
))))
342 (number-dispatch ((number number
))
343 (((foreach fixnum bignum
))
345 (complex (log (- number
)) (coerce pi
'single-float
))
346 (coerce (/ (log2 number
) (log (exp 1.0d0
) 2.0d0
)) 'single-float
)))
349 (complex (log (- number
)) (coerce pi
'single-float
))
350 (let ((numerator (numerator number
))
351 (denominator (denominator number
)))
352 (if (= (integer-length numerator
)
353 (integer-length denominator
))
354 (coerce (%log1p
(coerce (- number
1) 'double-float
))
356 (coerce (/ (- (log2 numerator
) (log2 denominator
))
357 (log (exp 1.0d0
) 2.0d0
))
359 (((foreach single-float double-float
))
360 ;; Is (log -0) -infinity (libm.a) or -infinity + i*pi (Kahan)?
361 ;; Since this doesn't seem to be an implementation issue
362 ;; I (pw) take the Kahan result.
363 (if (< (float-sign number
)
364 (coerce 0 '(dispatch-type number
)))
365 (complex (log (- number
)) (coerce pi
'(dispatch-type number
)))
366 (coerce (%log
(coerce number
'double-float
))
367 '(dispatch-type number
))))
369 (complex-log number
)))))
373 "Return the square root of NUMBER."
374 (number-dispatch ((number number
))
375 (((foreach fixnum bignum ratio
))
377 (complex-sqrt number
)
378 (coerce (%sqrt
(coerce number
'double-float
)) 'single-float
)))
379 (((foreach single-float double-float
))
381 (complex-sqrt (complex number
))
382 (coerce (%sqrt
(coerce number
'double-float
))
383 '(dispatch-type number
))))
385 (complex-sqrt number
))))
387 ;;;; trigonometic and related functions
391 "Return the absolute value of the number."
392 (number-dispatch ((number number
))
393 (((foreach single-float double-float fixnum rational
))
396 (let ((rx (realpart number
))
397 (ix (imagpart number
)))
400 (sqrt (+ (* rx rx
) (* ix ix
))))
402 (coerce (%hypot
(coerce rx
'double-float
)
403 (coerce ix
'double-float
))
408 (defun phase (number)
410 "Return the angle part of the polar representation of a complex number.
411 For complex numbers, this is (atan (imagpart number) (realpart number)).
412 For non-complex positive numbers, this is 0. For non-complex negative
417 (coerce pi
'single-float
)
420 (if (minusp (float-sign number
))
421 (coerce pi
'single-float
)
424 (if (minusp (float-sign number
))
425 (coerce pi
'double-float
)
428 (atan (imagpart number
) (realpart number
)))))
432 "Return the sine of NUMBER."
433 (number-dispatch ((number number
))
434 (handle-reals %sin number
)
436 (let ((x (realpart number
))
437 (y (imagpart number
)))
438 (complex (* (sin x
) (cosh y
))
439 (* (cos x
) (sinh y
)))))))
443 "Return the cosine of NUMBER."
444 (number-dispatch ((number number
))
445 (handle-reals %cos number
)
447 (let ((x (realpart number
))
448 (y (imagpart number
)))
449 (complex (* (cos x
) (cosh y
))
450 (- (* (sin x
) (sinh y
))))))))
454 "Return the tangent of NUMBER."
455 (number-dispatch ((number number
))
456 (handle-reals %tan number
)
458 (complex-tan number
))))
462 "Return cos(Theta) + i sin(Theta), i.e. exp(i Theta)."
463 (declare (type real theta
))
464 (complex (cos theta
) (sin theta
)))
468 "Return the arc sine of NUMBER."
469 (number-dispatch ((number number
))
471 (if (or (> number
1) (< number -
1))
472 (complex-asin number
)
473 (coerce (%asin
(coerce number
'double-float
)) 'single-float
)))
474 (((foreach single-float double-float
))
475 (if (or (> number
(coerce 1 '(dispatch-type number
)))
476 (< number
(coerce -
1 '(dispatch-type number
))))
477 (complex-asin (complex number
))
478 (coerce (%asin
(coerce number
'double-float
))
479 '(dispatch-type number
))))
481 (complex-asin number
))))
485 "Return the arc cosine of NUMBER."
486 (number-dispatch ((number number
))
488 (if (or (> number
1) (< number -
1))
489 (complex-acos number
)
490 (coerce (%acos
(coerce number
'double-float
)) 'single-float
)))
491 (((foreach single-float double-float
))
492 (if (or (> number
(coerce 1 '(dispatch-type number
)))
493 (< number
(coerce -
1 '(dispatch-type number
))))
494 (complex-acos (complex number
))
495 (coerce (%acos
(coerce number
'double-float
))
496 '(dispatch-type number
))))
498 (complex-acos number
))))
500 (defun atan (y &optional
(x nil xp
))
502 "Return the arc tangent of Y if X is omitted or Y/X if X is supplied."
505 (declare (type double-float y x
)
506 (values double-float
))
509 (if (plusp (float-sign x
))
512 (float-sign y
(/ pi
2)))
514 (number-dispatch ((y real
) (x real
))
516 (foreach double-float single-float fixnum bignum ratio
))
517 (atan2 y
(coerce x
'double-float
)))
518 (((foreach single-float fixnum bignum ratio
)
520 (atan2 (coerce y
'double-float
) x
))
521 (((foreach single-float fixnum bignum ratio
)
522 (foreach single-float fixnum bignum ratio
))
523 (coerce (atan2 (coerce y
'double-float
) (coerce x
'double-float
))
525 (number-dispatch ((y number
))
526 (handle-reals %atan y
)
530 ;;; It seems that every target system has a C version of sinh, cosh,
531 ;;; and tanh. Let's use these for reals because the original
532 ;;; implementations based on the definitions lose big in round-off
533 ;;; error. These bad definitions also mean that sin and cos for
534 ;;; complex numbers can also lose big.
538 "Return the hyperbolic sine of NUMBER."
539 (number-dispatch ((number number
))
540 (handle-reals %sinh number
)
542 (let ((x (realpart number
))
543 (y (imagpart number
)))
544 (complex (* (sinh x
) (cos y
))
545 (* (cosh x
) (sin y
)))))))
549 "Return the hyperbolic cosine of NUMBER."
550 (number-dispatch ((number number
))
551 (handle-reals %cosh number
)
553 (let ((x (realpart number
))
554 (y (imagpart number
)))
555 (complex (* (cosh x
) (cos y
))
556 (* (sinh x
) (sin y
)))))))
560 "Return the hyperbolic tangent of NUMBER."
561 (number-dispatch ((number number
))
562 (handle-reals %tanh number
)
564 (complex-tanh number
))))
566 (defun asinh (number)
568 "Return the hyperbolic arc sine of NUMBER."
569 (number-dispatch ((number number
))
570 (handle-reals %asinh number
)
572 (complex-asinh number
))))
574 (defun acosh (number)
576 "Return the hyperbolic arc cosine of NUMBER."
577 (number-dispatch ((number number
))
579 ;; acosh is complex if number < 1
581 (complex-acosh number
)
582 (coerce (%acosh
(coerce number
'double-float
)) 'single-float
)))
583 (((foreach single-float double-float
))
584 (if (< number
(coerce 1 '(dispatch-type number
)))
585 (complex-acosh (complex number
))
586 (coerce (%acosh
(coerce number
'double-float
))
587 '(dispatch-type number
))))
589 (complex-acosh number
))))
591 (defun atanh (number)
593 "Return the hyperbolic arc tangent of NUMBER."
594 (number-dispatch ((number number
))
596 ;; atanh is complex if |number| > 1
597 (if (or (> number
1) (< number -
1))
598 (complex-atanh number
)
599 (coerce (%atanh
(coerce number
'double-float
)) 'single-float
)))
600 (((foreach single-float double-float
))
601 (if (or (> number
(coerce 1 '(dispatch-type number
)))
602 (< number
(coerce -
1 '(dispatch-type number
))))
603 (complex-atanh (complex number
))
604 (coerce (%atanh
(coerce number
'double-float
))
605 '(dispatch-type number
))))
607 (complex-atanh number
))))
609 ;;; HP-UX does not supply a C version of log1p, so use the definition.
611 ;;; FIXME: This is really not a good definition. As per Raymond Toy
612 ;;; working on CMU CL, "The definition really loses big-time in
613 ;;; roundoff as x gets small."
615 #!-sb-fluid
(declaim (inline %log1p
))
617 (defun %log1p
(number)
618 (declare (double-float number
)
619 (optimize (speed 3) (safety 0)))
620 (the double-float
(log (the (double-float 0d0
) (+ number
1d0
)))))
622 ;;;; not-OLD-SPECFUN stuff
624 ;;;; (This was conditional on #-OLD-SPECFUN in the CMU CL sources,
625 ;;;; but OLD-SPECFUN was mentioned nowhere else, so it seems to be
626 ;;;; the standard special function system.)
628 ;;;; This is a set of routines that implement many elementary
629 ;;;; transcendental functions as specified by ANSI Common Lisp. The
630 ;;;; implementation is based on Kahan's paper.
632 ;;;; I believe I have accurately implemented the routines and are
633 ;;;; correct, but you may want to check for your self.
635 ;;;; These functions are written for CMU Lisp and take advantage of
636 ;;;; some of the features available there. It may be possible,
637 ;;;; however, to port this to other Lisps.
639 ;;;; Some functions are significantly more accurate than the original
640 ;;;; definitions in CMU Lisp. In fact, some functions in CMU Lisp
641 ;;;; give the wrong answer like (acos #c(-2.0 0.0)), where the true
642 ;;;; answer is pi + i*log(2-sqrt(3)).
644 ;;;; All of the implemented functions will take any number for an
645 ;;;; input, but the result will always be a either a complex
646 ;;;; single-float or a complex double-float.
648 ;;;; general functions:
660 ;;;; utility functions:
663 ;;;; internal functions:
664 ;;;; square coerce-to-complex-type cssqs complex-log-scaled
667 ;;;; Kahan, W. "Branch Cuts for Complex Elementary Functions, or Much
668 ;;;; Ado About Nothing's Sign Bit" in Iserles and Powell (eds.) "The
669 ;;;; State of the Art in Numerical Analysis", pp. 165-211, Clarendon
672 ;;;; The original CMU CL code requested:
673 ;;;; Please send any bug reports, comments, or improvements to
674 ;;;; Raymond Toy at <email address deleted during 2002 spam avalanche>.
676 ;;; FIXME: In SBCL, the floating point infinity constants like
677 ;;; SB!EXT:DOUBLE-FLOAT-POSITIVE-INFINITY aren't available as
678 ;;; constants at cross-compile time, because the cross-compilation
679 ;;; host might not have support for floating point infinities. Thus,
680 ;;; they're effectively implemented as special variable references,
681 ;;; and the code below which uses them might be unnecessarily
682 ;;; inefficient. Perhaps some sort of MAKE-LOAD-TIME-VALUE hackery
683 ;;; should be used instead? (KLUDGED 2004-03-08 CSR, by replacing the
684 ;;; special variable references with (probably equally slow)
687 ;;; FIXME: As of 2004-05, when PFD noted that IMAGPART and COMPLEX
688 ;;; differ in their interpretations of the real line, IMAGPART was
689 ;;; patch, which without a certain amount of effort would have altered
690 ;;; all the branch cut treatment. Clients of these COMPLEX- routines
691 ;;; were patched to use explicit COMPLEX, rather than implicitly
692 ;;; passing in real numbers for treatment with IMAGPART, and these
693 ;;; COMPLEX- functions altered to require arguments of type COMPLEX;
694 ;;; however, someone needs to go back to Kahan for the definitive
695 ;;; answer for treatment of negative real floating point numbers and
696 ;;; branch cuts. If adjustment is needed, it is probably the removal
697 ;;; of explicit calls to COMPLEX in the clients of irrational
698 ;;; functions. -- a slightly bitter CSR, 2004-05-16
700 (declaim (inline square
))
702 (declare (double-float x
))
705 ;;; original CMU CL comment, apparently re. SCALB and LOGB and
707 ;;; If you have these functions in libm, perhaps they should be used
708 ;;; instead of these Lisp versions. These versions are probably good
709 ;;; enough, especially since they are portable.
711 ;;; Compute 2^N * X without computing 2^N first. (Use properties of
712 ;;; the underlying floating-point format.)
713 (declaim (inline scalb
))
715 (declare (type double-float x
)
716 (type double-float-exponent n
))
719 ;;; This is like LOGB, but X is not infinity and non-zero and not a
720 ;;; NaN, so we can always return an integer.
721 (declaim (inline logb-finite
))
722 (defun logb-finite (x)
723 (declare (type double-float x
))
724 (multiple-value-bind (signif exponent sign
)
726 (declare (ignore signif sign
))
727 ;; DECODE-FLOAT is almost right, except that the exponent is off
731 ;;; Compute an integer N such that 1 <= |2^N * x| < 2.
732 ;;; For the special cases, the following values are used:
735 ;;; +/- infinity +infinity
738 (declare (type double-float x
))
739 (cond ((float-nan-p x
)
741 ((float-infinity-p x
)
742 ;; DOUBLE-FLOAT-POSITIVE-INFINITY
743 (double-from-bits 0 (1+ sb
!vm
:double-float-normal-exponent-max
) 0))
745 ;; The answer is negative infinity, but we are supposed to
746 ;; signal divide-by-zero, so do the actual division
752 ;;; This function is used to create a complex number of the
753 ;;; appropriate type:
754 ;;; Create complex number with real part X and imaginary part Y
755 ;;; such that has the same type as Z. If Z has type (complex
756 ;;; rational), the X and Y are coerced to single-float.
757 #!+long-float
(eval-when (:compile-toplevel
:load-toplevel
:execute
)
758 (error "needs work for long float support"))
759 (declaim (inline coerce-to-complex-type
))
760 (defun coerce-to-complex-type (x y z
)
761 (declare (double-float x y
)
763 (if (typep (realpart z
) 'double-float
)
765 ;; Convert anything that's not already a DOUBLE-FLOAT (because
766 ;; the initial argument was a (COMPLEX DOUBLE-FLOAT) and we
767 ;; haven't done anything to lose precision) to a SINGLE-FLOAT.
768 (complex (float x
1f0
)
771 ;;; Compute |(x+i*y)/2^k|^2 scaled to avoid over/underflow. The
772 ;;; result is r + i*k, where k is an integer.
773 #!+long-float
(eval-when (:compile-toplevel
:load-toplevel
:execute
)
774 (error "needs work for long float support"))
776 (let ((x (float (realpart z
) 1d0
))
777 (y (float (imagpart z
) 1d0
)))
778 ;; Would this be better handled using an exception handler to
779 ;; catch the overflow or underflow signal? For now, we turn all
780 ;; traps off and look at the accrued exceptions to see if any
781 ;; signal would have been raised.
782 (with-float-traps-masked (:underflow
:overflow
)
783 (let ((rho (+ (square x
) (square y
))))
784 (declare (optimize (speed 3) (space 0)))
785 (cond ((and (or (float-nan-p rho
)
786 (float-infinity-p rho
))
787 (or (float-infinity-p (abs x
))
788 (float-infinity-p (abs y
))))
789 ;; DOUBLE-FLOAT-POSITIVE-INFINITY
791 (double-from-bits 0 (1+ sb
!vm
:double-float-normal-exponent-max
) 0)
793 ((let ((threshold #.
(/ least-positive-double-float
794 double-float-epsilon
))
795 (traps (ldb sb
!vm
::float-sticky-bits
796 (sb!vm
:floating-point-modes
))))
797 ;; Overflow raised or (underflow raised and rho <
799 (or (not (zerop (logand sb
!vm
:float-overflow-trap-bit traps
)))
800 (and (not (zerop (logand sb
!vm
:float-underflow-trap-bit
803 ;; If we're here, neither x nor y are infinity and at
804 ;; least one is non-zero.. Thus logb returns a nice
806 (let ((k (- (logb-finite (max (abs x
) (abs y
))))))
807 (values (+ (square (scalb x k
))
808 (square (scalb y k
)))
813 ;;; principal square root of Z
815 ;;; Z may be RATIONAL or COMPLEX; the result is always a COMPLEX.
816 (defun complex-sqrt (z)
817 ;; KLUDGE: Here and below, we can't just declare Z to be of type
818 ;; COMPLEX, because one-arg COMPLEX on rationals returns a rational.
819 ;; Since there isn't a rational negative zero, this is OK from the
820 ;; point of view of getting the right answer in the face of branch
821 ;; cuts, but declarations of the form (OR RATIONAL COMPLEX) are
822 ;; still ugly. -- CSR, 2004-05-16
823 (declare (type (or complex rational
) z
))
824 (multiple-value-bind (rho k
)
826 (declare (type (or (member 0d0
) (double-float 0d0
)) rho
)
828 (let ((x (float (realpart z
) 1.0d0
))
829 (y (float (imagpart z
) 1.0d0
))
832 (declare (double-float x y eta nu
))
835 ;; space 0 to get maybe-inline functions inlined.
836 (declare (optimize (speed 3) (space 0)))
838 (if (not (float-nan-p x
))
839 (setf rho
(+ (scalb (abs x
) (- k
)) (sqrt rho
))))
844 (setf k
(1- (ash k -
1)))
845 (setf rho
(+ rho rho
))))
847 (setf rho
(scalb (sqrt rho
) k
))
853 (when (not (float-infinity-p (abs nu
)))
854 (setf nu
(/ (/ nu rho
) 2d0
)))
857 (setf nu
(float-sign y rho
))))
858 (coerce-to-complex-type eta nu z
)))))
860 ;;; Compute log(2^j*z).
862 ;;; This is for use with J /= 0 only when |z| is huge.
863 (defun complex-log-scaled (z j
)
864 (declare (type (or rational complex
) z
)
866 ;; The constants t0, t1, t2 should be evaluated to machine
867 ;; precision. In addition, Kahan says the accuracy of log1p
868 ;; influences the choices of these constants but doesn't say how to
869 ;; choose them. We'll just assume his choices matches our
870 ;; implementation of log1p.
871 (let ((t0 #.
(/ 1 (sqrt 2.0d0
)))
875 (x (float (realpart z
) 1.0d0
))
876 (y (float (imagpart z
) 1.0d0
)))
877 (multiple-value-bind (rho k
)
879 (declare (optimize (speed 3)))
880 (let ((beta (max (abs x
) (abs y
)))
881 (theta (min (abs x
) (abs y
))))
882 (coerce-to-complex-type (if (and (zerop k
)
886 (/ (%log1p
(+ (* (- beta
1.0d0
)
895 ;;; log of Z = log |Z| + i * arg Z
897 ;;; Z may be any number, but the result is always a complex.
898 (defun complex-log (z)
899 (declare (type (or rational complex
) z
))
900 (complex-log-scaled z
0))
902 ;;; KLUDGE: Let us note the following "strange" behavior. atanh 1.0d0
903 ;;; is +infinity, but the following code returns approx 176 + i*pi/4.
904 ;;; The reason for the imaginary part is caused by the fact that arg
905 ;;; i*y is never 0 since we have positive and negative zeroes. -- rtoy
906 ;;; Compute atanh z = (log(1+z) - log(1-z))/2.
907 (defun complex-atanh (z)
908 (declare (type (or rational complex
) z
))
910 (theta (/ (sqrt most-positive-double-float
) 4.0d0
))
911 (rho (/ 4.0d0
(sqrt most-positive-double-float
)))
912 (half-pi (/ pi
2.0d0
))
913 (rp (float (realpart z
) 1.0d0
))
914 (beta (float-sign rp
1.0d0
))
916 (y (* beta
(- (float (imagpart z
) 1.0d0
))))
919 ;; Shouldn't need this declare.
920 (declare (double-float x y
))
922 (declare (optimize (speed 3)))
923 (cond ((or (> x theta
)
925 ;; To avoid overflow...
926 (setf nu
(float-sign y half-pi
))
927 ;; ETA is real part of 1/(x + iy). This is x/(x^2+y^2),
928 ;; which can cause overflow. Arrange this computation so
929 ;; that it won't overflow.
930 (setf eta
(let* ((x-bigger (> x
(abs y
)))
931 (r (if x-bigger
(/ y x
) (/ x y
)))
932 (d (+ 1.0d0
(* r r
))))
937 ;; Should this be changed so that if y is zero, eta is set
938 ;; to +infinity instead of approx 176? In any case
939 ;; tanh(176) is 1.0d0 within working precision.
940 (let ((t1 (+ 4d0
(square y
)))
941 (t2 (+ (abs y
) rho
)))
942 (setf eta
(log (/ (sqrt (sqrt t1
))
946 (+ half-pi
(atan (* 0.5d0 t2
))))))))
948 (let ((t1 (+ (abs y
) rho
)))
949 ;; Normal case using log1p(x) = log(1 + x)
951 (%log1p
(/ (* 4.0d0 x
)
952 (+ (square (- 1.0d0 x
))
959 (coerce-to-complex-type (* beta eta
)
963 ;;; Compute tanh z = sinh z / cosh z.
964 (defun complex-tanh (z)
965 (declare (type (or rational complex
) z
))
966 (let ((x (float (realpart z
) 1.0d0
))
967 (y (float (imagpart z
) 1.0d0
)))
969 ;; space 0 to get maybe-inline functions inlined
970 (declare (optimize (speed 3) (space 0)))
972 ;; FIXME: this form is hideously broken wrt
973 ;; cross-compilation portability. Much else in this
974 ;; file is too, of course, sometimes hidden by
975 ;; constant-folding, but this one in particular clearly
976 ;; depends on host and target
977 ;; MOST-POSITIVE-DOUBLE-FLOATs being equal. -- CSR,
980 (log most-positive-double-float
))
982 (coerce-to-complex-type (float-sign x
)
986 (beta (+ 1.0d0
(* tv tv
)))
988 (rho (sqrt (+ 1.0d0
(* s s
)))))
989 (if (float-infinity-p (abs tv
))
990 (coerce-to-complex-type (/ rho s
)
993 (let ((den (+ 1.0d0
(* beta s s
))))
994 (coerce-to-complex-type (/ (* beta rho s
)
999 ;;; Compute acos z = pi/2 - asin z.
1001 ;;; Z may be any NUMBER, but the result is always a COMPLEX.
1002 (defun complex-acos (z)
1003 ;; Kahan says we should only compute the parts needed. Thus, the
1004 ;; REALPART's below should only compute the real part, not the whole
1005 ;; complex expression. Doing this can be important because we may get
1006 ;; spurious signals that occur in the part that we are not using.
1008 ;; However, we take a pragmatic approach and just use the whole
1011 ;; NOTE: The formula given by Kahan is somewhat ambiguous in whether
1012 ;; it's the conjugate of the square root or the square root of the
1013 ;; conjugate. This needs to be checked.
1015 ;; I checked. It doesn't matter because (conjugate (sqrt z)) is the
1016 ;; same as (sqrt (conjugate z)) for all z. This follows because
1018 ;; (conjugate (sqrt z)) = exp(0.5*log |z|)*exp(-0.5*j*arg z).
1020 ;; (sqrt (conjugate z)) = exp(0.5*log|z|)*exp(0.5*j*arg conj z)
1022 ;; and these two expressions are equal if and only if arg conj z =
1023 ;; -arg z, which is clearly true for all z.
1024 (declare (type (or rational complex
) z
))
1025 (let ((sqrt-1+z
(complex-sqrt (+ 1 z
)))
1026 (sqrt-1-z (complex-sqrt (- 1 z
))))
1027 (with-float-traps-masked (:divide-by-zero
)
1028 (complex (* 2 (atan (/ (realpart sqrt-1-z
)
1029 (realpart sqrt-1
+z
))))
1030 (asinh (imagpart (* (conjugate sqrt-1
+z
)
1033 ;;; Compute acosh z = 2 * log(sqrt((z+1)/2) + sqrt((z-1)/2))
1035 ;;; Z may be any NUMBER, but the result is always a COMPLEX.
1036 (defun complex-acosh (z)
1037 (declare (type (or rational complex
) z
))
1038 (let ((sqrt-z-1 (complex-sqrt (- z
1)))
1039 (sqrt-z+1 (complex-sqrt (+ z
1))))
1040 (with-float-traps-masked (:divide-by-zero
)
1041 (complex (asinh (realpart (* (conjugate sqrt-z-1
)
1043 (* 2 (atan (/ (imagpart sqrt-z-1
)
1044 (realpart sqrt-z
+1))))))))
1046 ;;; Compute asin z = asinh(i*z)/i.
1048 ;;; Z may be any NUMBER, but the result is always a COMPLEX.
1049 (defun complex-asin (z)
1050 (declare (type (or rational complex
) z
))
1051 (let ((sqrt-1-z (complex-sqrt (- 1 z
)))
1052 (sqrt-1+z
(complex-sqrt (+ 1 z
))))
1053 (with-float-traps-masked (:divide-by-zero
)
1054 (complex (atan (/ (realpart z
)
1055 (realpart (* sqrt-1-z sqrt-1
+z
))))
1056 (asinh (imagpart (* (conjugate sqrt-1-z
)
1059 ;;; Compute asinh z = log(z + sqrt(1 + z*z)).
1061 ;;; Z may be any number, but the result is always a complex.
1062 (defun complex-asinh (z)
1063 (declare (type (or rational complex
) z
))
1064 ;; asinh z = -i * asin (i*z)
1065 (let* ((iz (complex (- (imagpart z
)) (realpart z
)))
1066 (result (complex-asin iz
)))
1067 (complex (imagpart result
)
1068 (- (realpart result
)))))
1070 ;;; Compute atan z = atanh (i*z) / i.
1072 ;;; Z may be any number, but the result is always a complex.
1073 (defun complex-atan (z)
1074 (declare (type (or rational complex
) z
))
1075 ;; atan z = -i * atanh (i*z)
1076 (let* ((iz (complex (- (imagpart z
)) (realpart z
)))
1077 (result (complex-atanh iz
)))
1078 (complex (imagpart result
)
1079 (- (realpart result
)))))
1081 ;;; Compute tan z = -i * tanh(i * z)
1083 ;;; Z may be any number, but the result is always a complex.
1084 (defun complex-tan (z)
1085 (declare (type (or rational complex
) z
))
1086 ;; tan z = -i * tanh(i*z)
1087 (let* ((iz (complex (- (imagpart z
)) (realpart z
)))
1088 (result (complex-tanh iz
)))
1089 (complex (imagpart result
)
1090 (- (realpart result
)))))