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
19 (eval-when (:compile-toplevel
:load-toplevel
:execute
)
20 (defun handle-reals (function var
)
21 `((((foreach fixnum single-float bignum ratio
))
22 (coerce (,function
(coerce ,var
'double-float
)) 'single-float
))
26 (defun handle-complex (form)
27 `((((foreach (complex double-float
) (complex single-float
) (complex rational
)))
30 ;;; Make these INLINE, since the call to C is at least as compact as a
31 ;;; Lisp call, and saves number consing to boot.
32 (defmacro def-math-rtn
(name num-args
&optional wrapper
)
33 (let ((function (symbolicate "%" (string-upcase name
)))
34 (args (loop for i below num-args
35 collect
(intern (format nil
"ARG~D" i
)))))
37 (declaim (inline ,function
))
38 (defun ,function
,args
39 (declare (sb-c:flushable sb-c
:%alien-funcall
))
40 (truly-the ;; avoid checking the result
41 ,(type-specifier (fun-type-returns (info :function
:type function
)))
43 (extern-alien ,(format nil
"~:[~;sb_~]~a" wrapper name
)
44 (function double-float
45 ,@(loop repeat num-args
46 collect
'double-float
)))
50 #+x86
;; for constant folding
51 (macrolet ((def (name ll
)
52 `(defun ,name
,ll
(,name
,@ll
))))
67 #+(or x86-64 arm-vfp arm64 riscv
) ;; for constant folding
68 (macrolet ((def (name ll
)
69 `(defun ,name
,ll
(,name
,@ll
))))
72 ;;;; stubs for the Unix math library
74 ;;;; Many of these are unnecessary on the X86 because they're built
78 #-x86
(def-math-rtn "sin" 1)
79 #-x86
(def-math-rtn "cos" 1)
80 #-x86
(def-math-rtn "tan" 1)
81 #-x86
(def-math-rtn "atan" 1)
82 #-x86
(def-math-rtn "atan2" 2)
84 ;;; See src/runtime/wrap.c for the definitions of the "sb_"-prefixed things.
85 (def-math-rtn "acos" 1 #+win32 t
)
86 (def-math-rtn "asin" 1 #+win32 t
)
87 (def-math-rtn "cosh" 1 #+win32 t
)
88 (def-math-rtn "sinh" 1 #+win32 t
)
89 (def-math-rtn "tanh" 1 #+win32 t
)
90 (def-math-rtn "asinh" 1 #+win32 t
)
91 (def-math-rtn "acosh" 1 #+win32 t
)
92 (def-math-rtn "atanh" 1 #+win32 t
)
94 ;;; exponential and logarithmic
95 (def-math-rtn "hypot" 2 #+win32 t
)
96 #-x86
(def-math-rtn "exp" 1)
97 #-x86
(def-math-rtn "log" 1)
98 #-x86
(def-math-rtn "log10" 1)
99 (def-math-rtn "pow" 2)
100 #-
(or x86 x86-64 arm-vfp arm64 riscv
) (def-math-rtn "sqrt" 1)
101 #-x86
(def-math-rtn "log1p" 1)
102 #-x86
(def-math-rtn "log2" 1)
108 "Return e raised to the power NUMBER."
109 (declare (explicit-check))
110 (number-dispatch ((number number
))
111 (handle-reals %exp number
)
113 (* (exp (realpart number
))
114 (cis (imagpart number
))))))
116 ;;; INTEXP -- Handle the rational base, integer power case.
117 ;;; This function precisely calculates base raised to an integral
118 ;;; power. It separates the cases by the sign of power, for efficiency
119 ;;; reasons, as powers can be calculated more efficiently if power is
120 ;;; a positive integer. Values of power are calculated as positive
121 ;;; integers, and inverted if negative.
122 (defun intexp (base power
)
123 (declare (explicit-check))
131 (cond ((= power
0) 1)
133 (t (error 'division-by-zero
134 :operands
(list 0 power
)
137 (let ((den (denominator base
))
138 (num (numerator base
)))
140 (let ((negated (- power
)))
142 (intexp den negated
))
144 (intexp (- den
) negated
))
146 (build-ratio (intexp den negated
)
147 (intexp num negated
)))))
148 (build-ratio (intexp num power
)
149 (intexp den power
)))))
151 (/ (intexp base
(- power
))))
157 (nextn (ash power -
1) (ash power -
1))
158 (total (if (oddp power
) base
1)
159 (if (oddp power
) (* base total
) total
)))
160 ((zerop nextn
) total
)
161 (setq base
(* base base
))
162 (setq power nextn
)))))
164 ;;; If an integer power of a rational, use INTEXP above. Otherwise, do
165 ;;; floating point stuff. If both args are real, we try %POW right
166 ;;; off, assuming it will return 0 if the result may be complex. If
167 ;;; so, we call COMPLEX-POW which directly computes the complex
168 ;;; result. We also separate the complex-real and real-complex cases
169 ;;; from the general complex case.
170 (defun expt (base power
)
171 "Return BASE raised to the POWER."
172 (declare (explicit-check))
174 (if (and (zerop base
) (floatp power
))
175 (error 'arguments-out-of-domain-error
176 :operands
(list base power
)
178 :references
'((:ansi-cl
:function expt
)))
179 (let ((result (1+ (* base power
))))
180 (if (and (floatp result
) (float-nan-p result
))
183 (labels (;; determine if the double float is an integer.
184 ;; 0 - not an integer
188 (declare (type (unsigned-byte 31) ihi
)
189 (type (unsigned-byte 32) lo
)
190 (optimize (speed 3) (safety 0)))
192 (declare (type fixnum isint
))
193 (cond ((>= ihi
#x43400000
) ; exponent >= 53
196 (let ((k (- (ash ihi -
20) #x3ff
))) ; exponent
197 (declare (type (mod 53) k
))
199 (let* ((shift (- 52 k
))
200 (j (logand (ash lo
(- shift
))))
202 (declare (type (mod 32) shift
)
203 (type (unsigned-byte 32) j j2
))
205 (setq isint
(- 2 (logand j
1))))))
207 (let* ((shift (- 20 k
))
208 (j (ash ihi
(- shift
)))
210 (declare (type (mod 32) shift
)
211 (type (unsigned-byte 31) j j2
))
213 (setq isint
(- 2 (logand j
1))))))))))
215 (real-expt (x y rtype
)
216 (let ((x (coerce x
'double-float
))
217 (y (coerce y
'double-float
)))
218 (declare (double-float x y
))
219 (let* ((x-hi (double-float-high-bits x
))
220 (x-lo (double-float-low-bits x
))
221 (x-ihi (logand x-hi
#x7fffffff
))
222 (y-hi (double-float-high-bits y
))
223 (y-lo (double-float-low-bits y
))
224 (y-ihi (logand y-hi
#x7fffffff
)))
225 (declare (type (signed-byte 32) x-hi y-hi
)
226 (type (unsigned-byte 31) x-ihi y-ihi
)
227 (type (unsigned-byte 32) x-lo y-lo
))
229 (when (zerop (logior y-ihi y-lo
))
230 (return-from real-expt
(coerce 1d0 rtype
)))
232 ;; FIXME: Hardcoded qNaN/sNaN values are not portable.
233 (when (or (> x-ihi
#x7ff00000
)
234 (and (= x-ihi
#x7ff00000
) (/= x-lo
0))
236 (and (= y-ihi
#x7ff00000
) (/= y-lo
0)))
237 (return-from real-expt
(coerce (+ x y
) rtype
)))
238 (let ((yisint (if (< x-hi
0) (isint y-ihi y-lo
) 0)))
239 (declare (type fixnum yisint
))
240 ;; special value of y
241 (when (and (zerop y-lo
) (= y-ihi
#x7ff00000
))
243 (return-from real-expt
244 (cond ((and (= x-ihi
#x3ff00000
) (zerop x-lo
))
246 (coerce (- y y
) rtype
))
247 ((>= x-ihi
#x3ff00000
)
248 ;; (|x|>1)**+-inf = inf,0
253 ;; (|x|<1)**-,+inf = inf,0
256 (coerce 0 rtype
))))))
258 (let ((abs-x (abs x
)))
259 (declare (double-float abs-x
))
260 ;; special value of x
261 (when (and (zerop x-lo
)
262 (or (= x-ihi
#x7ff00000
) (zerop x-ihi
)
263 (= x-ihi
#x3ff00000
)))
264 ;; x is +-0,+-inf,+-1
265 (let ((z (if (< y-hi
0)
266 (/ 1 abs-x
) ; z = (1/|x|)
268 (declare (double-float z
))
270 (cond ((and (= x-ihi
#x3ff00000
) (zerop yisint
))
272 (let ((y*pi
(* y pi
)))
273 (declare (double-float y
*pi
))
274 (return-from real-expt
276 (coerce (%cos y
*pi
) rtype
)
277 (coerce (%sin y
*pi
) rtype
)))))
279 ;; (x<0)**odd = -(|x|**odd)
281 (return-from real-expt
(coerce z rtype
))))
285 (coerce (%pow x y
) rtype
)
287 (let ((pow (%pow abs-x y
)))
288 (declare (double-float pow
))
291 (coerce (* -
1d0 pow
) rtype
))
295 (let ((y*pi
(* y pi
)))
296 (declare (double-float y
*pi
))
298 (coerce (* pow
(%cos y
*pi
))
300 (coerce (* pow
(%sin y
*pi
))
302 (complex-expt (base power
)
303 (if (and (zerop base
) (plusp (realpart power
)))
305 (exp (* power
(log base
))))))
306 (declare (inline real-expt complex-expt
))
307 (number-dispatch ((base number
) (power number
))
308 (((foreach fixnum
(or bignum ratio
) (complex rational
)) integer
)
310 (((foreach single-float double-float
) rational
)
311 (real-expt base power
'(dispatch-type base
)))
312 (((foreach fixnum
(or bignum ratio
) single-float
)
313 (foreach ratio single-float
))
314 (real-expt base power
'single-float
))
315 (((foreach fixnum
(or bignum ratio
) single-float double-float
)
317 (real-expt base power
'double-float
))
318 ((double-float single-float
)
319 (real-expt base power
'double-float
))
320 ;; Handle (expt <complex> <rational>), except the case dealt with
321 ;; in the first clause above, (expt <(complex rational)> <integer>).
322 (((foreach (complex rational
))
324 (* (expt (abs base
) power
)
325 (cis (* power
(phase base
)))))
326 (((foreach (complex single-float
) (complex double-float
))
327 (foreach fixnum
(or bignum ratio
)))
328 (* (expt (abs base
) power
)
329 (cis (* power
(phase base
)))))
330 ;; The next three clauses handle (expt <real> <complex>).
331 (((foreach fixnum
(or bignum ratio
) single-float
)
332 (foreach (complex single-float
) (complex rational
)))
333 (complex-expt base power
))
334 (((foreach fixnum
(or bignum ratio
) single-float
)
335 (complex double-float
))
336 (complex-expt (coerce base
'double-float
) power
))
337 ((double-float complex
)
338 (complex-expt base power
))
339 ;; The next three clauses handle (expt <complex> <float>) and
340 ;; (expt <complex> <complex>).
341 (((foreach (complex single-float
) (complex rational
))
342 (foreach (complex single-float
) (complex rational
) single-float
))
343 (complex-expt base power
))
344 (((foreach (complex single-float
) (complex rational
))
345 (foreach (complex double-float
) double-float
))
346 (complex-expt (coerce base
'(complex double-float
)) power
))
347 (((complex double-float
)
348 (foreach complex double-float single-float
))
349 (complex-expt base power
))))))
351 (declaim (start-block log
))
352 (defun log2/nonnegative-integer
(x)
353 (declare (type (integer 0) x
))
356 ;; Write x = 2^n*f where 1/2 <= f < 1. Then log2(x) = n +
357 ;; log2(f). So we grab the top few bits of x and scale that
358 ;; appropriately, take the log of it and add it to n.
360 ;; Motivated by an attempt to get LOG to work better on bignums.
361 (cond ((typep x
'sb-vm
:signed-word
)
362 (%log2
(coerce x
'double-float
)))
364 (let ((n (integer-length x
)))
365 (cond ((< n sb-vm
:double-float-digits
)
366 (%log2
(coerce x
'double-float
)))
368 (let ((f (ldb (byte sb-vm
:double-float-digits
369 (- n sb-vm
:double-float-digits
))
371 (+ n
(%log2
(scale-float (coerce f
'double-float
)
372 (- sb-vm
:double-float-digits
)))))))))))
374 (defun log2/double-float
(x)
375 (declare (type double-float x
))
376 (if (float-sign-bit-set-p x
)
377 (complex (%log2
(- x
)) (sb-xc:/ pi
(log 2d0
)))
380 (defun log2/nonnegative-ratio
(x)
381 (declare (type (and (rational 0) (not integer
)) x
))
382 (truly-the double-float
383 (- (truly-the double-float
(log2/nonnegative-integer
(numerator x
)))
384 (truly-the double-float
(log2/nonnegative-integer
(denominator x
))))))
386 (defun log2/nonnegative-rational
(x)
387 (declare (type (rational 0) x
))
388 (if (typep x
'integer
)
389 (log2/nonnegative-integer x
)
390 (log2/nonnegative-ratio x
)))
392 (defun log2/rational
(x)
393 (declare (type rational x
))
395 (complex (log2/nonnegative-rational
(- x
)) (sb-xc:/ pi
(log 2d0
)))
396 (log2/nonnegative-rational x
)))
398 (defun log (number &optional
(base nil base-p
))
399 "Return the logarithm of NUMBER in the base BASE, which defaults to e."
400 (declare (explicit-check))
404 (if (or (typep number
'(or double-float
(complex double-float
)))
405 (typep base
'(or double-float
(complex double-float
))))
409 (error 'division-by-zero
:operands
(list number base
) :operation
'log
))
410 ((and (typep number
'(rational 0))
411 (typep base
'(rational 0)))
412 (coerce (/ (truly-the double-float
(log2/nonnegative-rational number
))
413 (truly-the double-float
(log2/nonnegative-rational base
)))
416 (number-dispatch ((number number
) (base number
))
417 ;; (log <real> <real>)
419 ;; we excluded real results from rational arguments above
420 (coerce (/ (log2/rational number
) (log2/rational base
))
421 '(complex single-float
)))
422 ((rational (foreach single-float double-float
))
423 (let ((result (/ (log2/rational number
) (log2/double-float
(coerce base
'double-float
)))))
425 (coerce result
'(dispatch-type base
))
426 (coerce result
'(complex (dispatch-type base
))))))
427 (((foreach single-float double-float
) rational
)
428 (let ((result (/ (log2/double-float
(coerce number
'double-float
)) (log2/rational base
))))
430 (coerce result
'(dispatch-type number
))
431 (coerce result
'(complex (dispatch-type number
))))))
432 ((single-float single-float
)
433 (/ (log number
) (log base
)))
434 ((single-float double-float
)
435 (/ (log (coerce number
'double-float
)) (log base
)))
436 ((double-float single-float
)
437 (/ (log number
) (log (coerce base
'double-float
))))
438 ((double-float double-float
)
439 (/ (log number
) (log base
)))
440 ;; complex single-float result
441 (((foreach rational single-float
) (foreach (complex rational
) (complex single-float
)))
442 (/ (log number
) (log base
)))
443 (((foreach (complex rational
) (complex single-float
))
444 (foreach rational single-float
(complex rational
) (complex single-float
)))
445 (/ (log number
) (log base
)))
446 ;; complex double-float result, from contagion
447 (((foreach double-float
(complex double-float
))
448 (foreach (complex rational
) (complex single-float
)))
449 (/ (log number
) (log (coerce base
'(complex double-float
)))))
450 (((foreach (complex rational
) (complex single-float
))
451 (foreach double-float
(complex double-float
)))
452 (/ (log (coerce number
'(complex double-float
))) (log base
)))
453 (((foreach rational single-float
) (complex double-float
))
454 (/ (log (coerce number
'double-float
)) (log base
)))
455 (((complex double-float
) (foreach rational single-float
))
456 (/ (log number
) (log (coerce base
'double-float
))))
457 ;; complex double-float result, no contagion
458 ((double-float (complex double-float
))
459 (/ (log number
) (log base
)))
460 (((complex double-float
) (foreach double-float
(complex double-float
)))
461 (/ (log number
) (log base
))))))
462 (let ((log2e 1.4426950408889634d0
))
464 (error 'division-by-zero
:operands
(list number
) :operation
'log
))
465 (number-dispatch ((number number
))
466 (((foreach fixnum bignum
))
468 (complex (log (- number
)) (coerce pi
'single-float
))
469 (coerce (/ (truly-the double-float
(log2/nonnegative-integer number
))
474 (complex (log (- number
)) (coerce pi
'single-float
))
475 (let ((numerator (numerator number
))
476 (denominator (denominator number
)))
477 (if (and (<= -
1 (- (integer-length numerator
) (integer-length denominator
)) 1)
479 (coerce (%log1p
(coerce (- number
1) 'double-float
))
481 (coerce (/ (- (truly-the double-float
(log2/nonnegative-integer numerator
))
482 (truly-the double-float
(log2/nonnegative-integer denominator
)))
485 (((foreach single-float double-float
))
486 ;; Is (log -0) -infinity (libm.a) or -infinity + i*pi (Kahan)?
487 ;; Since this doesn't seem to be an implementation issue
488 ;; I (pw) take the Kahan result.
489 (if (float-sign-bit-set-p number
)
490 (complex (log (- number
)) (coerce pi
'(dispatch-type number
)))
491 (coerce (%log
(coerce number
'double-float
))
492 '(dispatch-type number
))))
494 (complex-log number
))))))
495 (declaim (end-block))
498 "Return the square root of NUMBER."
499 (declare (explicit-check))
500 (number-dispatch ((number number
))
501 (((foreach fixnum bignum ratio
))
504 (coerce (%sqrt
(- (coerce number
'double-float
))) 'single-float
))
505 (coerce (%sqrt
(coerce number
'double-float
)) 'single-float
)))
506 (((foreach single-float double-float
))
508 (complex (coerce 0.0 '(dispatch-type number
))
509 (coerce (%sqrt
(- (coerce number
'double-float
)))
510 '(dispatch-type number
)))
511 (coerce (%sqrt
(coerce number
'double-float
))
512 '(dispatch-type number
))))
514 (complex-sqrt number
))))
516 ;;;; trigonometic and related functions
519 "Return the absolute value of the number."
520 (declare (explicit-check))
521 (number-dispatch ((number number
))
522 (((foreach single-float double-float fixnum rational
))
525 (let ((rx (realpart number
))
526 (ix (imagpart number
)))
529 (sqrt (+ (* rx rx
) (* ix ix
))))
531 (coerce (%hypot
(coerce rx
'double-float
)
532 (coerce (truly-the single-float ix
) 'double-float
))
535 (%hypot rx
(truly-the double-float ix
))))))))
537 (defun phase (number)
538 "Return the angle part of the polar representation of a complex number.
539 For complex numbers, this is (atan (imagpart number) (realpart number)).
540 For non-complex positive numbers, this is 0. For non-complex negative
542 (declare (explicit-check))
543 (number-dispatch ((number number
))
546 (coerce pi
'single-float
)
549 (if (minusp (float-sign number
))
550 (coerce pi
'single-float
)
553 (if (minusp (float-sign number
))
554 (coerce pi
'double-float
)
557 (atan (imagpart number
) (realpart number
)))))
560 "Return the sine of NUMBER."
561 (declare (explicit-check))
562 (number-dispatch ((number number
))
563 (handle-reals %sin number
)
565 (let ((x (realpart number
))
566 (y (imagpart number
)))
567 (complex (* (sin x
) (cosh y
))
568 (* (cos x
) (sinh y
)))))))
571 "Return the cosine of NUMBER."
572 (declare (explicit-check))
573 (number-dispatch ((number number
))
574 (handle-reals %cos number
)
576 (let ((x (realpart number
))
577 (y (imagpart number
)))
578 (complex (* (cos x
) (cosh y
))
579 (- (* (sin x
) (sinh y
))))))))
582 "Return the tangent of NUMBER."
583 (declare (explicit-check))
584 (number-dispatch ((number number
))
585 (handle-reals %tan number
)
587 ;; tan z = -i * tanh(i*z)
588 (let* ((result (complex-tanh (complex (- (imagpart number
))
589 (realpart number
)))))
590 (complex (imagpart result
)
591 (- (realpart result
)))))))
594 "Return cos(Theta) + i sin(Theta), i.e. exp(i Theta)."
595 (declare (explicit-check ))
596 (number-dispatch ((theta real
))
597 (((foreach single-float double-float rational
))
598 (complex (cos theta
) (sin theta
)))))
601 "Return the arc sine of NUMBER."
602 (declare (explicit-check))
603 (number-dispatch ((number number
))
605 (if (or (> number
1) (< number -
1))
606 (complex-asin number
)
607 (coerce (%asin
(coerce number
'double-float
)) 'single-float
)))
608 (((foreach single-float double-float
))
609 (if (or (> number
(coerce 1 '(dispatch-type number
)))
610 (< number
(coerce -
1 '(dispatch-type number
))))
611 (complex-asin (complex number
))
612 (coerce (%asin
(coerce number
'double-float
))
613 '(dispatch-type number
))))
615 (complex-asin number
))))
618 "Return the arc cosine of NUMBER."
619 (declare (explicit-check))
620 (number-dispatch ((number number
))
622 (if (or (> number
1) (< number -
1))
623 (complex-acos number
)
624 (coerce (%acos
(coerce number
'double-float
)) 'single-float
)))
625 (((foreach single-float double-float
))
626 (if (or (> number
(coerce 1 '(dispatch-type number
)))
627 (< number
(coerce -
1 '(dispatch-type number
))))
628 (complex-acos (complex number
))
629 (coerce (%acos
(coerce number
'double-float
))
630 '(dispatch-type number
))))
632 (complex-acos number
))))
634 (defun atan (y &optional
(x nil xp
))
635 "Return the arc tangent of Y if X is omitted or Y/X if X is supplied."
636 (declare (explicit-check))
639 (declare (type double-float y x
)
640 (values double-float
))
643 (if (not (float-sign-bit-set-p x
))
646 (float-sign y
(sb-xc:/ pi
2)))
648 (number-dispatch ((y real
) (x real
))
650 (foreach double-float single-float fixnum bignum ratio
))
651 (atan2 y
(coerce x
'double-float
)))
652 (((foreach single-float fixnum bignum ratio
)
654 (atan2 (coerce y
'double-float
) x
))
655 (((foreach single-float fixnum bignum ratio
)
656 (foreach single-float fixnum bignum ratio
))
657 (coerce (atan2 (coerce y
'double-float
) (coerce x
'double-float
))
659 (number-dispatch ((y number
))
660 (handle-reals %atan y
)
664 ;;; It seems that every target system has a C version of sinh, cosh,
665 ;;; and tanh. Let's use these for reals because the original
666 ;;; implementations based on the definitions lose big in round-off
667 ;;; error. These bad definitions also mean that sin and cos for
668 ;;; complex numbers can also lose big.
671 "Return the hyperbolic sine of NUMBER."
672 (declare (explicit-check))
673 (number-dispatch ((number number
))
674 (handle-reals %sinh number
)
676 (let ((x (realpart number
))
677 (y (imagpart number
)))
678 (complex (* (sinh x
) (cos y
))
679 (* (cosh x
) (sin y
)))))))
682 "Return the hyperbolic cosine of NUMBER."
683 (declare (explicit-check))
684 (number-dispatch ((number number
))
685 (handle-reals %cosh number
)
687 (let ((x (realpart number
))
688 (y (imagpart number
)))
689 (complex (* (cosh x
) (cos y
))
690 (* (sinh x
) (sin y
)))))))
693 "Return the hyperbolic tangent of NUMBER."
694 (declare (explicit-check))
695 (number-dispatch ((number number
))
696 (handle-reals %tanh number
)
698 (complex-tanh number
))))
700 (defun asinh (number)
701 "Return the hyperbolic arc sine of NUMBER."
702 (declare (explicit-check))
703 (number-dispatch ((number number
))
704 (handle-reals %asinh number
)
706 (complex-asinh number
))))
708 (defun acosh (number)
709 "Return the hyperbolic arc cosine of NUMBER."
710 (declare (explicit-check))
711 (number-dispatch ((number number
))
713 ;; acosh is complex if number < 1
715 (complex-acosh number
)
716 (coerce (%acosh
(coerce number
'double-float
)) 'single-float
)))
717 (((foreach single-float double-float
))
718 (if (< number
(coerce 1 '(dispatch-type number
)))
719 (complex-acosh (complex number
))
720 (coerce (%acosh
(coerce number
'double-float
))
721 '(dispatch-type number
))))
723 (complex-acosh number
))))
725 (defun atanh (number)
726 "Return the hyperbolic arc tangent of NUMBER."
727 (declare (explicit-check))
728 (number-dispatch ((number number
))
730 ;; atanh is complex if |number| > 1
731 (if (or (> number
1) (< number -
1))
732 (complex-atanh number
)
733 (coerce (%atanh
(coerce number
'double-float
)) 'single-float
)))
734 (((foreach single-float double-float
))
735 (if (or (> number
(coerce 1 '(dispatch-type number
)))
736 (< number
(coerce -
1 '(dispatch-type number
))))
737 (complex-atanh (complex number
))
738 (coerce (%atanh
(coerce number
'double-float
))
739 '(dispatch-type number
))))
741 (complex-atanh number
))))
744 ;;;; not-OLD-SPECFUN stuff
746 ;;;; (This was conditional on #-OLD-SPECFUN in the CMU CL sources,
747 ;;;; but OLD-SPECFUN was mentioned nowhere else, so it seems to be
748 ;;;; the standard special function system.)
750 ;;;; This is a set of routines that implement many elementary
751 ;;;; transcendental functions as specified by ANSI Common Lisp. The
752 ;;;; implementation is based on Kahan's paper.
754 ;;;; I believe I have accurately implemented the routines and are
755 ;;;; correct, but you may want to check for your self.
757 ;;;; These functions are written for CMU Lisp and take advantage of
758 ;;;; some of the features available there. It may be possible,
759 ;;;; however, to port this to other Lisps.
761 ;;;; Some functions are significantly more accurate than the original
762 ;;;; definitions in CMU Lisp. In fact, some functions in CMU Lisp
763 ;;;; give the wrong answer like (acos #c(-2.0 0.0)), where the true
764 ;;;; answer is pi + i*log(2-sqrt(3)).
766 ;;;; All of the implemented functions will take any number for an
767 ;;;; input, but the result will always be a either a complex
768 ;;;; single-float or a complex double-float.
770 ;;;; general functions:
781 ;;;; utility functions:
784 ;;;; internal functions:
785 ;;;; square coerce-to-complex-type cssqs
788 ;;;; Kahan, W. "Branch Cuts for Complex Elementary Functions, or Much
789 ;;;; Ado About Nothing's Sign Bit" in Iserles and Powell (eds.) "The
790 ;;;; State of the Art in Numerical Analysis", pp. 165-211, Clarendon
793 ;;;; The original CMU CL code requested:
794 ;;;; Please send any bug reports, comments, or improvements to
795 ;;;; Raymond Toy at <email address deleted during 2002 spam avalanche>.
797 ;;; FIXME: As of 2004-05, when PFD noted that IMAGPART and COMPLEX
798 ;;; differ in their interpretations of the real line, IMAGPART was
799 ;;; patch, which without a certain amount of effort would have altered
800 ;;; all the branch cut treatment. Clients of these COMPLEX- routines
801 ;;; were patched to use explicit COMPLEX, rather than implicitly
802 ;;; passing in real numbers for treatment with IMAGPART, and these
803 ;;; COMPLEX- functions altered to require arguments of type COMPLEX;
804 ;;; however, someone needs to go back to Kahan for the definitive
805 ;;; answer for treatment of negative real floating point numbers and
806 ;;; branch cuts. If adjustment is needed, it is probably the removal
807 ;;; of explicit calls to COMPLEX in the clients of irrational
808 ;;; functions. -- a slightly bitter CSR, 2004-05-16
810 (declaim (inline square
))
812 (declare (double-float x
))
815 ;;; original CMU CL comment, apparently re. LOGB and
817 ;;; If you have these functions in libm, perhaps they should be used
818 ;;; instead of these Lisp versions. These versions are probably good
819 ;;; enough, especially since they are portable.
821 ;;; This is like LOGB, but X is not infinity and non-zero and not a
822 ;;; NaN, so we can always return an integer.
823 (declaim (inline logb-finite
))
824 (defun logb-finite (x)
825 (declare (type double-float x
))
826 (multiple-value-bind (signif exponent sign
)
828 (declare (ignore signif sign
))
829 ;; DECODE-FLOAT is almost right, except that the exponent is off
833 ;;; Compute an integer N such that 1 <= |2^N * x| < 2.
834 ;;; For the special cases, the following values are used:
837 ;;; +/- infinity +infinity
840 (declare (type double-float x
))
841 (cond ((float-nan-p x
)
843 ((float-infinity-p x
) double-float-positive-infinity
)
845 ;; The answer is negative infinity, but we are supposed to
846 ;; signal divide-by-zero, so do the actual division
851 ;;; This function is used to create a complex number of the
852 ;;; appropriate type:
853 ;;; Create complex number with real part X and imaginary part Y
854 ;;; such that has the same type as Z. If Z has type (complex
855 ;;; rational), the X and Y are coerced to single-float.
856 #+long-float
(eval-when (:compile-toplevel
:load-toplevel
:execute
)
857 (error "needs work for long float support"))
858 (declaim (inline coerce-to-complex-type
))
859 (defun coerce-to-complex-type (x y z
)
860 (declare (double-float x y
)
862 (if (typep (realpart z
) 'double-float
)
864 ;; Convert anything that's not already a DOUBLE-FLOAT (because
865 ;; the initial argument was a (COMPLEX DOUBLE-FLOAT) and we
866 ;; haven't done anything to lose precision) to a SINGLE-FLOAT.
867 (complex (float x
1f0
)
870 ;;; Compute |(x+i*y)/2^k|^2 scaled to avoid over/underflow. The
871 ;;; result is r + i*k, where k is an integer.
872 #+long-float
(eval-when (:compile-toplevel
:load-toplevel
:execute
)
873 (error "needs work for long float support"))
875 (declare (muffle-conditions compiler-note
))
876 (let ((x (float (realpart z
) 1d0
))
877 (y (float (imagpart z
) 1d0
)))
878 ;; Would this be better handled using an exception handler to
879 ;; catch the overflow or underflow signal? For now, we turn all
880 ;; traps off and look at the accrued exceptions to see if any
881 ;; signal would have been raised.
882 (with-float-traps-masked (:underflow
:overflow
)
883 (let ((rho (+ (square x
) (square y
))))
884 (declare (optimize (speed 3) (space 0)))
885 (cond ((and (float-infinity-or-nan-p rho
)
886 (or (float-infinity-p x
) (float-infinity-p y
)))
887 (values double-float-positive-infinity
0))
889 ;; (/ least-positive-double-float double-float-epsilon)
891 (make-double-float #x1fffff
#xfffffffe
)
893 (error "(/ least-positive-long-float long-float-epsilon)"))
894 (traps (ldb sb-vm
:float-sticky-bits
895 (sb-vm:floating-point-modes
))))
896 ;; Overflow raised or (underflow raised and rho <
898 (or (not (zerop (logand sb-vm
:float-overflow-trap-bit traps
)))
899 (and (not (zerop (logand sb-vm
:float-underflow-trap-bit
902 ;; If we're here, neither x nor y are infinity and at
903 ;; least one is non-zero.. Thus logb returns a nice
905 (let ((k (- (logb-finite (max (abs x
) (abs y
))))))
906 (values (+ (square (scale-float x k
))
907 (square (scale-float y k
)))
912 ;;; principal square root of Z
914 ;;; Z may be RATIONAL or COMPLEX; the result is always a COMPLEX.
915 (defun complex-sqrt (z)
916 ;; KLUDGE: Here and below, we can't just declare Z to be of type
917 ;; COMPLEX, because one-arg COMPLEX on rationals returns a rational.
918 ;; Since there isn't a rational negative zero, this is OK from the
919 ;; point of view of getting the right answer in the face of branch
920 ;; cuts, but declarations of the form (OR RATIONAL COMPLEX) are
921 ;; still ugly. -- CSR, 2004-05-16
922 (declare (type (or complex rational
) z
))
923 (multiple-value-bind (rho k
)
925 (declare (type (or (member 0d0
) (double-float 0d0
)) rho
)
927 (let ((x (float (realpart z
) 1.0d0
))
928 (y (float (imagpart z
) 1.0d0
))
931 (declare (double-float x y eta nu
)
932 ;; get maybe-inline functions inlined.
933 (optimize (space 0)))
934 (if (not (float-nan-p x
))
935 (setf rho
(+ (scale-float (abs x
) (- k
)) (sqrt rho
))))
940 (setf k
(1- (ash k -
1)))
941 (setf rho
(+ rho rho
))))
943 (setf rho
(scale-float (sqrt rho
) k
))
949 (when (not (float-infinity-p nu
))
950 (setf nu
(/ (/ nu rho
) 2d0
)))
953 (setf nu
(float-sign y rho
))))
954 (coerce-to-complex-type eta nu z
))))
956 ;;; log of Z = log |Z| + i * arg Z
958 ;;; Z may be any number, but the result is always a complex.
959 (defun complex-log (z)
960 (declare (muffle-conditions compiler-note
))
961 (declare (type (or rational complex
) z
))
962 ;; The constants t0, t1, t2 should be evaluated to machine
963 ;; precision. In addition, Kahan says the accuracy of log1p
964 ;; influences the choices of these constants but doesn't say how to
965 ;; choose them. We'll just assume his choices matches our
966 ;; implementation of log1p.
967 (let ((t0 #-long-float
(make-double-float #x3fe6a09e
#x667f3bcd
)
968 #+long-float
(error "(/ (sqrt 2l0) 2)"))
971 (ln2 #-long-float
(make-double-float #x3fe62e42
#xfefa39ef
)
972 #+long-float
(error "(log 2l0)"))
973 (x (float (realpart z
) 1.0d0
))
974 (y (float (imagpart z
) 1.0d0
)))
975 (multiple-value-bind (rho k
)
977 (declare (optimize (speed 3)))
978 (let ((beta (max (abs x
) (abs y
)))
979 (theta (min (abs x
) (abs y
))))
980 (coerce-to-complex-type (if (and (zerop k
)
984 (/ (%log1p
(+ (* (- beta
1.0d0
)
993 ;;; KLUDGE: Let us note the following "strange" behavior. atanh 1.0d0
994 ;;; is +infinity, but the following code returns approx 176 + i*pi/4.
995 ;;; The reason for the imaginary part is caused by the fact that arg
996 ;;; i*y is never 0 since we have positive and negative zeroes. -- rtoy
997 ;;; Compute atanh z = (log(1+z) - log(1-z))/2.
998 (defun complex-atanh (z)
999 (declare (muffle-conditions compiler-note
))
1000 (declare (type (or rational complex
) z
))
1002 (theta (sb-xc:/ (sb-xc:sqrt most-positive-double-float
) 4.0d0
))
1003 (rho (sb-xc:/ 4.0d0
(sb-xc:sqrt most-positive-double-float
)))
1004 (half-pi (sb-xc:/ pi
2.0d0
))
1005 (rp (float (realpart z
) 1.0d0
))
1006 (beta (float-sign rp
1.0d0
))
1008 (y (* beta
(- (float (imagpart z
) 1.0d0
))))
1011 ;; Shouldn't need this declare.
1012 (declare (double-float x y
))
1014 (declare (optimize (speed 3)))
1015 (cond ((or (> x theta
)
1017 ;; To avoid overflow...
1018 (setf nu
(float-sign y half-pi
))
1019 ;; ETA is real part of 1/(x + iy). This is x/(x^2+y^2),
1020 ;; which can cause overflow. Arrange this computation so
1021 ;; that it won't overflow.
1022 (setf eta
(let* ((x-bigger (> x
(abs y
)))
1023 (r (if x-bigger
(/ y x
) (/ x y
)))
1024 (d (+ 1.0d0
(* r r
))))
1029 ;; Should this be changed so that if y is zero, eta is set
1030 ;; to +infinity instead of approx 176? In any case
1031 ;; tanh(176) is 1.0d0 within working precision.
1032 (let ((t1 (+ 4d0
(square y
)))
1033 (t2 (+ (abs y
) rho
)))
1034 (setf eta
(log (/ (sqrt (sqrt t1
))
1038 (+ half-pi
(atan (* 0.5d0 t2
))))))))
1040 (let ((t1 (+ (abs y
) rho
)))
1041 ;; Normal case using log1p(x) = log(1 + x)
1043 (%log1p
(/ (* 4.0d0 x
)
1044 (+ (square (- 1.0d0 x
))
1051 (coerce-to-complex-type (* beta eta
)
1056 (format t
"~x~%" (double-float-bits (/ (+ (log 2l0) (log most-positive-long-float
)) 4l0)))
1064 double most_pos_dbl
= 1.7976931348623157e308
;
1065 double thing
= (log(most_pos_dbl) + log
(2.0e0
)) / 4.0e0
;
1067 memcpy
(&word
, &thing
, 8);
1068 printf
("%lX = %20.15lf\n", word
, thing
);
1070 prints
: 406633CE8FB9F87D
= 177.618965018485966
1073 ;;; Compute tanh z = sinh z / cosh z.
1074 (defun complex-tanh (z)
1075 (declare (muffle-conditions compiler-note
))
1076 (declare (type (or rational complex
) z
))
1077 (let ((x (float (realpart z
) 1.0d0
))
1078 (y (float (imagpart z
) 1.0d0
)))
1080 ;; space 0 to get maybe-inline functions inlined
1081 (declare (optimize (speed 3) (space 0)))
1083 #-long-float
(make-double-float #x406633ce
#x8fb9f87d
)
1084 #+long-float
(error "(/ (+ (log 2l0) (log most-positive-long-float)) 4l0)"))
1085 (coerce-to-complex-type (float-sign x
)
1088 (let* ((tv (%tan y
))
1089 (beta (+ 1.0d0
(* tv tv
)))
1091 (rho (sqrt (+ 1.0d0
(* s s
)))))
1092 (if (float-infinity-p tv
)
1093 (coerce-to-complex-type (/ rho s
)
1096 (let ((den (+ 1.0d0
(* beta s s
))))
1097 (coerce-to-complex-type (/ (* beta rho s
)
1102 ;;; Compute acos z = pi/2 - asin z.
1104 ;;; Z may be any NUMBER, but the result is always a COMPLEX.
1105 (defun complex-acos (z)
1106 ;; Kahan says we should only compute the parts needed. Thus, the
1107 ;; REALPART's below should only compute the real part, not the whole
1108 ;; complex expression. Doing this can be important because we may get
1109 ;; spurious signals that occur in the part that we are not using.
1111 ;; However, we take a pragmatic approach and just use the whole
1114 ;; NOTE: The formula given by Kahan is somewhat ambiguous in whether
1115 ;; it's the conjugate of the square root or the square root of the
1116 ;; conjugate. This needs to be checked.
1118 ;; I checked. It doesn't matter because (conjugate (sqrt z)) is the
1119 ;; same as (sqrt (conjugate z)) for all z. This follows because
1121 ;; (conjugate (sqrt z)) = exp(0.5*log |z|)*exp(-0.5*j*arg z).
1123 ;; (sqrt (conjugate z)) = exp(0.5*log|z|)*exp(0.5*j*arg conj z)
1125 ;; and these two expressions are equal if and only if arg conj z =
1126 ;; -arg z, which is clearly true for all z.
1127 (declare (type (or rational complex
) z
))
1128 (let ((sqrt-1+z
(complex-sqrt (+ 1 z
)))
1129 (sqrt-1-z (complex-sqrt (- 1 z
))))
1130 (with-float-traps-masked (:divide-by-zero
)
1131 (complex (* 2 (atan (/ (realpart sqrt-1-z
)
1132 (realpart sqrt-1
+z
))))
1133 (asinh (imagpart (* (conjugate sqrt-1
+z
)
1136 ;;; Compute acosh z = 2 * log(sqrt((z+1)/2) + sqrt((z-1)/2))
1138 ;;; Z may be any NUMBER, but the result is always a COMPLEX.
1139 (defun complex-acosh (z)
1140 (declare (type (or rational complex
) z
))
1141 (let ((sqrt-z-1 (complex-sqrt (- z
1)))
1142 (sqrt-z+1 (complex-sqrt (+ z
1))))
1143 (with-float-traps-masked (:divide-by-zero
)
1144 (complex (asinh (realpart (* (conjugate sqrt-z-1
)
1146 (* 2 (atan (/ (imagpart sqrt-z-1
)
1147 (realpart sqrt-z
+1))))))))
1149 ;;; Compute asin z = asinh(i*z)/i.
1151 ;;; Z may be any NUMBER, but the result is always a COMPLEX.
1152 (defun complex-asin (z)
1153 (declare (type (or rational complex
) z
))
1154 (let ((sqrt-1-z (complex-sqrt (- 1 z
)))
1155 (sqrt-1+z
(complex-sqrt (+ 1 z
))))
1156 (with-float-traps-masked (:divide-by-zero
)
1157 (complex (atan (/ (realpart z
)
1158 (realpart (* sqrt-1-z sqrt-1
+z
))))
1159 (asinh (imagpart (* (conjugate sqrt-1-z
)
1162 ;;; Compute asinh z = log(z + sqrt(1 + z*z)).
1164 ;;; Z may be any number, but the result is always a complex.
1165 (defun complex-asinh (z)
1166 (declare (type (or rational complex
) z
))
1167 ;; asinh z = -i * asin (i*z)
1168 (let* ((iz (complex (- (imagpart z
)) (realpart z
)))
1169 (result (complex-asin iz
)))
1170 (complex (imagpart result
)
1171 (- (realpart result
)))))
1173 ;;; Compute atan z = atanh (i*z) / i.
1175 ;;; Z may be any number, but the result is always a complex.
1176 (defun complex-atan (z)
1177 (declare (type (or rational complex
) z
))
1178 ;; atan z = -i * atanh (i*z)
1179 (let* ((iz (complex (- (imagpart z
)) (realpart z
)))
1180 (result (complex-atanh iz
)))
1181 (complex (imagpart result
)
1182 (- (realpart result
)))))