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
17 (defconstant pi 3.14159265358979323846264338327950288419716939937511L0)
18 ;(defconstant e 2.71828182845904523536028747135266249775724709369996L0)
20 ;;; Make these INLINE, since the call to C is at least as compact as a Lisp
21 ;;; call, and saves number consing to boot.
23 ;;; FIXME: This should be (EVAL-WHEN (COMPILE-EVAL) (SB!XC:DEFMACRO ..)),
25 (defmacro def-math-rtn (name num-args)
26 (let ((function (intern (concatenate 'simple-string
28 (string-upcase name)))))
30 (proclaim '(inline ,function))
31 (let ((sb!int::*rogue-export* "DEF-MATH-RTN"))
33 (sb!alien:def-alien-routine (,name ,function) double-float
34 ,@(let ((results nil))
35 (dotimes (i num-args (nreverse results))
36 (push (list (intern (format nil "ARG-~D" i))
40 (eval-when (:compile-toplevel :load-toplevel :execute)
42 (defun handle-reals (function var)
43 `((((foreach fixnum single-float bignum ratio))
44 (coerce (,function (coerce ,var 'double-float)) 'single-float))
50 ;;;; stubs for the Unix math library
52 ;;; Please refer to the Unix man pages for details about these routines.
55 #!-x86 (def-math-rtn "sin" 1)
56 #!-x86 (def-math-rtn "cos" 1)
57 #!-x86 (def-math-rtn "tan" 1)
58 (def-math-rtn "asin" 1)
59 (def-math-rtn "acos" 1)
60 #!-x86 (def-math-rtn "atan" 1)
61 #!-x86 (def-math-rtn "atan2" 2)
62 (def-math-rtn "sinh" 1)
63 (def-math-rtn "cosh" 1)
64 (def-math-rtn "tanh" 1)
65 (def-math-rtn "asinh" 1)
66 (def-math-rtn "acosh" 1)
67 (def-math-rtn "atanh" 1)
69 ;;; Exponential and Logarithmic.
70 #!-x86 (def-math-rtn "exp" 1)
71 #!-x86 (def-math-rtn "log" 1)
72 #!-x86 (def-math-rtn "log10" 1)
73 (def-math-rtn "pow" 2)
74 #!-x86 (def-math-rtn "sqrt" 1)
75 (def-math-rtn "hypot" 2)
76 #!-(or hpux x86) (def-math-rtn "log1p" 1)
78 #!+x86 ;; These are needed for use by byte-compiled files.
81 (declare (double-float x)
82 (values double-float))
85 (declare (double-float x)
86 (values double-float))
89 (declare (double-float x)
90 (values double-float))
93 (declare (double-float x)
94 (values double-float))
97 (declare (double-float x)
98 (values double-float))
100 (defun %tan-quick (x)
101 (declare (double-float x)
102 (values double-float))
105 (declare (double-float x)
106 (values double-float))
109 (declare (double-float x y)
110 (values double-float))
113 (declare (double-float x)
114 (values double-float))
117 (declare (double-float x)
118 (values double-float))
121 (declare (double-float x)
122 (values double-float))
126 (declare (type (double-float 0d0) x)
128 (values (double-float 0d0)))
131 (declare (double-float x)
132 (values double-float))
134 (defun %scalbn (f ex)
135 (declare (double-float f)
136 (type (signed-byte 32) ex)
137 (values double-float))
140 (declare (double-float f ex)
141 (values double-float))
144 (declare (double-float x)
145 (values double-float))
148 (declare (double-float x)
149 (values double-float))
157 "Return e raised to the power NUMBER."
158 (number-dispatch ((number number))
159 (handle-reals %exp number)
161 (* (exp (realpart number))
162 (cis (imagpart number))))))
164 ;;; INTEXP -- Handle the rational base, integer power case.
166 ;;; FIXME: As long as the
167 ;;; system dies on stack overflow or memory exhaustion, it seems reasonable
168 ;;; to have this, but its default should be NIL, and when it's NIL,
169 ;;; anything should be accepted.
170 (defparameter *intexp-maximum-exponent* 10000)
172 ;;; This function precisely calculates base raised to an integral power. It
173 ;;; separates the cases by the sign of power, for efficiency reasons, as powers
174 ;;; can be calculated more efficiently if power is a positive integer. Values
175 ;;; of power are calculated as positive integers, and inverted if negative.
176 (defun intexp (base power)
177 (when (> (abs power) *intexp-maximum-exponent*)
178 ;; FIXME: should be ordinary error, not CERROR. (Once we set the
179 ;; default for the variable to NIL, the un-continuable error will
180 ;; be less obnoxious.)
181 (cerror "Continue with calculation."
182 "The absolute value of ~S exceeds ~S."
183 power '*intexp-maximum-exponent* base power))
184 (cond ((minusp power)
185 (/ (intexp base (- power))))
189 (do ((nextn (ash power -1) (ash power -1))
190 (total (if (oddp power) base 1)
191 (if (oddp power) (* base total) total)))
192 ((zerop nextn) total)
193 (setq base (* base base))
194 (setq power nextn)))))
196 ;;; If an integer power of a rational, use INTEXP above. Otherwise, do
197 ;;; floating point stuff. If both args are real, we try %POW right off,
198 ;;; assuming it will return 0 if the result may be complex. If so, we call
199 ;;; COMPLEX-POW which directly computes the complex result. We also separate
200 ;;; the complex-real and real-complex cases from the general complex case.
201 (defun expt (base power)
203 "Returns BASE raised to the POWER."
206 (labels (;; determine if the double float is an integer.
207 ;; 0 - not an integer
211 (declare (type (unsigned-byte 31) ihi)
212 (type (unsigned-byte 32) lo)
213 (optimize (speed 3) (safety 0)))
215 (declare (type fixnum isint))
216 (cond ((>= ihi #x43400000) ; exponent >= 53
219 (let ((k (- (ash ihi -20) #x3ff))) ; exponent
220 (declare (type (mod 53) k))
222 (let* ((shift (- 52 k))
223 (j (logand (ash lo (- shift))))
225 (declare (type (mod 32) shift)
226 (type (unsigned-byte 32) j j2))
228 (setq isint (- 2 (logand j 1))))))
230 (let* ((shift (- 20 k))
231 (j (ash ihi (- shift)))
233 (declare (type (mod 32) shift)
234 (type (unsigned-byte 31) j j2))
236 (setq isint (- 2 (logand j 1))))))))))
238 (real-expt (x y rtype)
239 (let ((x (coerce x 'double-float))
240 (y (coerce y 'double-float)))
241 (declare (double-float x y))
242 (let* ((x-hi (sb!kernel:double-float-high-bits x))
243 (x-lo (sb!kernel:double-float-low-bits x))
244 (x-ihi (logand x-hi #x7fffffff))
245 (y-hi (sb!kernel:double-float-high-bits y))
246 (y-lo (sb!kernel:double-float-low-bits y))
247 (y-ihi (logand y-hi #x7fffffff)))
248 (declare (type (signed-byte 32) x-hi y-hi)
249 (type (unsigned-byte 31) x-ihi y-ihi)
250 (type (unsigned-byte 32) x-lo y-lo))
252 (when (zerop (logior y-ihi y-lo))
253 (return-from real-expt (coerce 1d0 rtype)))
255 (when (or (> x-ihi #x7ff00000)
256 (and (= x-ihi #x7ff00000) (/= x-lo 0))
258 (and (= y-ihi #x7ff00000) (/= y-lo 0)))
259 (return-from real-expt (coerce (+ x y) rtype)))
260 (let ((yisint (if (< x-hi 0) (isint y-ihi y-lo) 0)))
261 (declare (type fixnum yisint))
262 ;; special value of y
263 (when (and (zerop y-lo) (= y-ihi #x7ff00000))
265 (return-from real-expt
266 (cond ((and (= x-ihi #x3ff00000) (zerop x-lo))
268 (coerce (- y y) rtype))
269 ((>= x-ihi #x3ff00000)
270 ;; (|x|>1)**+-inf = inf,0
275 ;; (|x|<1)**-,+inf = inf,0
278 (coerce 0 rtype))))))
280 (let ((abs-x (abs x)))
281 (declare (double-float abs-x))
282 ;; special value of x
283 (when (and (zerop x-lo)
284 (or (= x-ihi #x7ff00000) (zerop x-ihi)
285 (= x-ihi #x3ff00000)))
286 ;; x is +-0,+-inf,+-1
287 (let ((z (if (< y-hi 0)
288 (/ 1 abs-x) ; z = (1/|x|)
290 (declare (double-float z))
292 (cond ((and (= x-ihi #x3ff00000) (zerop yisint))
294 (let ((y*pi (* y pi)))
295 (declare (double-float y*pi))
296 (return-from real-expt
298 (coerce (%cos y*pi) rtype)
299 (coerce (%sin y*pi) rtype)))))
301 ;; (x<0)**odd = -(|x|**odd)
303 (return-from real-expt (coerce z rtype))))
307 (coerce (sb!kernel::%pow x y) rtype)
309 (let ((pow (sb!kernel::%pow abs-x y)))
310 (declare (double-float pow))
313 (coerce (* -1d0 pow) rtype))
317 (let ((y*pi (* y pi)))
318 (declare (double-float y*pi))
320 (coerce (* pow (%cos y*pi)) rtype)
321 (coerce (* pow (%sin y*pi)) rtype)))))))))))))
322 (declare (inline real-expt))
323 (number-dispatch ((base number) (power number))
324 (((foreach fixnum (or bignum ratio) (complex rational)) integer)
326 (((foreach single-float double-float) rational)
327 (real-expt base power '(dispatch-type base)))
328 (((foreach fixnum (or bignum ratio) single-float)
329 (foreach ratio single-float))
330 (real-expt base power 'single-float))
331 (((foreach fixnum (or bignum ratio) single-float double-float)
333 (real-expt base power 'double-float))
334 ((double-float single-float)
335 (real-expt base power 'double-float))
336 (((foreach (complex rational) (complex float)) rational)
337 (* (expt (abs base) power)
338 (cis (* power (phase base)))))
339 (((foreach fixnum (or bignum ratio) single-float double-float)
341 (if (and (zerop base) (plusp (realpart power)))
343 (exp (* power (log base)))))
344 (((foreach (complex float) (complex rational))
345 (foreach complex double-float single-float))
346 (if (and (zerop base) (plusp (realpart power)))
348 (exp (* power (log base)))))))))
350 (defun log (number &optional (base nil base-p))
352 "Return the logarithm of NUMBER in the base BASE, which defaults to e."
356 (/ (log number) (log base)))
357 (number-dispatch ((number number))
358 (((foreach fixnum bignum ratio))
360 (complex (log (- number)) (coerce pi 'single-float))
361 (coerce (%log (coerce number 'double-float)) 'single-float)))
362 (((foreach single-float double-float))
363 ;; Is (log -0) -infinity (libm.a) or -infinity + i*pi (Kahan)?
364 ;; Since this doesn't seem to be an implementation issue
365 ;; I (pw) take the Kahan result.
366 (if (< (float-sign number)
367 (coerce 0 '(dispatch-type number)))
368 (complex (log (- number)) (coerce pi '(dispatch-type number)))
369 (coerce (%log (coerce number 'double-float))
370 '(dispatch-type number))))
372 (complex-log number)))))
376 "Return the square root of NUMBER."
377 (number-dispatch ((number number))
378 (((foreach fixnum bignum ratio))
380 (complex-sqrt number)
381 (coerce (%sqrt (coerce number 'double-float)) 'single-float)))
382 (((foreach single-float double-float))
384 (complex-sqrt number)
385 (coerce (%sqrt (coerce number 'double-float))
386 '(dispatch-type number))))
388 (complex-sqrt number))))
390 ;;;; trigonometic and related functions
394 "Returns the absolute value of the number."
395 (number-dispatch ((number number))
396 (((foreach single-float double-float fixnum rational))
399 (let ((rx (realpart number))
400 (ix (imagpart number)))
403 (sqrt (+ (* rx rx) (* ix ix))))
405 (coerce (%hypot (coerce rx 'double-float)
406 (coerce ix 'double-float))
411 (defun phase (number)
413 "Returns the angle part of the polar representation of a complex number.
414 For complex numbers, this is (atan (imagpart number) (realpart number)).
415 For non-complex positive numbers, this is 0. For non-complex negative
420 (coerce pi 'single-float)
423 (if (minusp (float-sign number))
424 (coerce pi 'single-float)
427 (if (minusp (float-sign number))
428 (coerce pi 'double-float)
431 (atan (imagpart number) (realpart number)))))
435 "Return the sine of NUMBER."
436 (number-dispatch ((number number))
437 (handle-reals %sin number)
439 (let ((x (realpart number))
440 (y (imagpart number)))
441 (complex (* (sin x) (cosh y))
442 (* (cos x) (sinh y)))))))
446 "Return the cosine of NUMBER."
447 (number-dispatch ((number number))
448 (handle-reals %cos number)
450 (let ((x (realpart number))
451 (y (imagpart number)))
452 (complex (* (cos x) (cosh y))
453 (- (* (sin x) (sinh y))))))))
457 "Return the tangent of NUMBER."
458 (number-dispatch ((number number))
459 (handle-reals %tan number)
461 (complex-tan number))))
465 "Return cos(Theta) + i sin(Theta), AKA exp(i Theta)."
467 (error "Argument to CIS is complex: ~S" theta)
468 (complex (cos theta) (sin theta))))
472 "Return the arc sine of NUMBER."
473 (number-dispatch ((number number))
475 (if (or (> number 1) (< number -1))
476 (complex-asin number)
477 (coerce (%asin (coerce number 'double-float)) 'single-float)))
478 (((foreach single-float double-float))
479 (if (or (> number (coerce 1 '(dispatch-type number)))
480 (< number (coerce -1 '(dispatch-type number))))
481 (complex-asin number)
482 (coerce (%asin (coerce number 'double-float))
483 '(dispatch-type number))))
485 (complex-asin number))))
489 "Return the arc cosine of NUMBER."
490 (number-dispatch ((number number))
492 (if (or (> number 1) (< number -1))
493 (complex-acos number)
494 (coerce (%acos (coerce number 'double-float)) 'single-float)))
495 (((foreach single-float double-float))
496 (if (or (> number (coerce 1 '(dispatch-type number)))
497 (< number (coerce -1 '(dispatch-type number))))
498 (complex-acos number)
499 (coerce (%acos (coerce number 'double-float))
500 '(dispatch-type number))))
502 (complex-acos number))))
504 (defun atan (y &optional (x nil xp))
506 "Return the arc tangent of Y if X is omitted or Y/X if X is supplied."
509 (declare (type double-float y x)
510 (values double-float))
513 (if (plusp (float-sign x))
516 (float-sign y (/ pi 2)))
518 (number-dispatch ((y number) (x number))
520 (foreach double-float single-float fixnum bignum ratio))
521 (atan2 y (coerce x 'double-float)))
522 (((foreach single-float fixnum bignum ratio)
524 (atan2 (coerce y 'double-float) x))
525 (((foreach single-float fixnum bignum ratio)
526 (foreach single-float fixnum bignum ratio))
527 (coerce (atan2 (coerce y 'double-float) (coerce x 'double-float))
529 (number-dispatch ((y number))
530 (handle-reals %atan y)
534 ;; It seems that everyone has a C version of sinh, cosh, and
535 ;; tanh. Let's use these for reals because the original
536 ;; implementations based on the definitions lose big in round-off
537 ;; error. These bad definitions also mean that sin and cos for
538 ;; complex numbers can also lose big.
543 "Return the hyperbolic sine of NUMBER."
544 (/ (- (exp number) (exp (- number))) 2))
548 "Return the hyperbolic sine of NUMBER."
549 (number-dispatch ((number number))
550 (handle-reals %sinh number)
552 (let ((x (realpart number))
553 (y (imagpart number)))
554 (complex (* (sinh x) (cos y))
555 (* (cosh x) (sin y)))))))
560 "Return the hyperbolic cosine of NUMBER."
561 (/ (+ (exp number) (exp (- number))) 2))
565 "Return the hyperbolic cosine of NUMBER."
566 (number-dispatch ((number number))
567 (handle-reals %cosh number)
569 (let ((x (realpart number))
570 (y (imagpart number)))
571 (complex (* (cosh x) (cos y))
572 (* (sinh x) (sin y)))))))
576 "Return the hyperbolic tangent of NUMBER."
577 (number-dispatch ((number number))
578 (handle-reals %tanh number)
580 (complex-tanh number))))
582 (defun asinh (number)
584 "Return the hyperbolic arc sine of NUMBER."
585 (number-dispatch ((number number))
586 (handle-reals %asinh number)
588 (complex-asinh number))))
590 (defun acosh (number)
592 "Return the hyperbolic arc cosine of NUMBER."
593 (number-dispatch ((number number))
595 ;; acosh is complex if number < 1
597 (complex-acosh number)
598 (coerce (%acosh (coerce number 'double-float)) 'single-float)))
599 (((foreach single-float double-float))
600 (if (< number (coerce 1 '(dispatch-type number)))
601 (complex-acosh number)
602 (coerce (%acosh (coerce number 'double-float))
603 '(dispatch-type number))))
605 (complex-acosh number))))
607 (defun atanh (number)
609 "Return the hyperbolic arc tangent of NUMBER."
610 (number-dispatch ((number number))
612 ;; atanh is complex if |number| > 1
613 (if (or (> number 1) (< number -1))
614 (complex-atanh number)
615 (coerce (%atanh (coerce number 'double-float)) 'single-float)))
616 (((foreach single-float double-float))
617 (if (or (> number (coerce 1 '(dispatch-type number)))
618 (< number (coerce -1 '(dispatch-type number))))
619 (complex-atanh number)
620 (coerce (%atanh (coerce number 'double-float))
621 '(dispatch-type number))))
623 (complex-atanh number))))
625 ;;; HP-UX does not supply a C version of log1p, so
626 ;;; use the definition.
629 #!-sb-fluid (declaim (inline %log1p))
631 (defun %log1p (number)
632 (declare (double-float number)
633 (optimize (speed 3) (safety 0)))
634 (the double-float (log (the (double-float 0d0) (+ number 1d0)))))