0.6.11.24:
[sbcl.git] / src / code / irrat.lisp
1 ;;;; This file contains all the irrational functions. (Actually, most
2 ;;;; of the work is done by calling out to C.)
3
4 ;;;; This software is part of the SBCL system. See the README file for
5 ;;;; more information.
6 ;;;;
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.
12
13 (in-package "SB!KERNEL")
14 \f
15 ;;;; miscellaneous constants, utility functions, and macros
16
17 (defconstant pi 3.14159265358979323846264338327950288419716939937511L0)
18 ;(defconstant e 2.71828182845904523536028747135266249775724709369996L0)
19
20 ;;; Make these INLINE, since the call to C is at least as compact as a
21 ;;; Lisp call, and saves number consing to boot.
22 (eval-when (:compile-toplevel :execute)
23
24 (sb!xc:defmacro def-math-rtn (name num-args)
25   (let ((function (symbolicate "%" (string-upcase name))))
26     `(progn
27        (proclaim '(inline ,function))
28        (sb!alien:def-alien-routine (,name ,function) double-float
29          ,@(let ((results nil))
30              (dotimes (i num-args (nreverse results))
31                (push (list (intern (format nil "ARG-~D" i))
32                            'double-float)
33                      results)))))))
34
35 (defun handle-reals (function var)
36   `((((foreach fixnum single-float bignum ratio))
37      (coerce (,function (coerce ,var 'double-float)) 'single-float))
38     ((double-float)
39      (,function ,var))))
40
41 ) ; EVAL-WHEN
42 \f
43 ;;;; stubs for the Unix math library
44
45 ;;; Please refer to the Unix man pages for details about these routines.
46
47 ;;; trigonometric
48 #!-x86 (def-math-rtn "sin" 1)
49 #!-x86 (def-math-rtn "cos" 1)
50 #!-x86 (def-math-rtn "tan" 1)
51 (def-math-rtn "asin" 1)
52 (def-math-rtn "acos" 1)
53 #!-x86 (def-math-rtn "atan" 1)
54 #!-x86 (def-math-rtn "atan2" 2)
55 (def-math-rtn "sinh" 1)
56 (def-math-rtn "cosh" 1)
57 (def-math-rtn "tanh" 1)
58 (def-math-rtn "asinh" 1)
59 (def-math-rtn "acosh" 1)
60 (def-math-rtn "atanh" 1)
61
62 ;;; exponential and logarithmic
63 #!-x86 (def-math-rtn "exp" 1)
64 #!-x86 (def-math-rtn "log" 1)
65 #!-x86 (def-math-rtn "log10" 1)
66 (def-math-rtn "pow" 2)
67 #!-x86 (def-math-rtn "sqrt" 1)
68 (def-math-rtn "hypot" 2)
69 #!-(or hpux x86) (def-math-rtn "log1p" 1)
70
71 #!+x86 ;; These are needed for use by byte-compiled files.
72 (progn
73   (defun %sin (x)
74     (declare (double-float x)
75              (values double-float))
76     (%sin x))
77   (defun %sin-quick (x)
78     (declare (double-float x)
79              (values double-float))
80     (%sin-quick x))
81   (defun %cos (x)
82     (declare (double-float x)
83              (values double-float))
84     (%cos x))
85   (defun %cos-quick (x)
86     (declare (double-float x)
87              (values double-float))
88     (%cos-quick x))
89   (defun %tan (x)
90     (declare (double-float x)
91              (values double-float))
92     (%tan x))
93   (defun %tan-quick (x)
94     (declare (double-float x)
95              (values double-float))
96     (%tan-quick x))
97   (defun %atan (x)
98     (declare (double-float x)
99              (values double-float))
100     (%atan x))
101   (defun %atan2 (x y)
102     (declare (double-float x y)
103              (values double-float))
104     (%atan2 x y))
105   (defun %exp (x)
106     (declare (double-float x)
107              (values double-float))
108     (%exp x))
109   (defun %log (x)
110     (declare (double-float x)
111              (values double-float))
112     (%log x))
113   (defun %log10 (x)
114     (declare (double-float x)
115              (values double-float))
116     (%log10 x))
117   #+nil ;; notyet
118   (defun %pow (x y)
119     (declare (type (double-float 0d0) x)
120              (double-float y)
121              (values (double-float 0d0)))
122     (%pow x y))
123   (defun %sqrt (x)
124     (declare (double-float x)
125              (values double-float))
126     (%sqrt x))
127   (defun %scalbn (f ex)
128     (declare (double-float f)
129              (type (signed-byte 32) ex)
130              (values double-float))
131     (%scalbn f ex))
132   (defun %scalb (f ex)
133     (declare (double-float f ex)
134              (values double-float))
135     (%scalb f ex))
136   (defun %logb (x)
137     (declare (double-float x)
138              (values double-float))
139     (%logb x))
140   (defun %log1p (x)
141     (declare (double-float x)
142              (values double-float))
143     (%log1p x))
144   ) ; progn
145 \f
146 ;;;; power functions
147
148 (defun exp (number)
149   #!+sb-doc
150   "Return e raised to the power NUMBER."
151   (number-dispatch ((number number))
152     (handle-reals %exp number)
153     ((complex)
154      (* (exp (realpart number))
155         (cis (imagpart number))))))
156
157 ;;; INTEXP -- Handle the rational base, integer power case.
158
159 ;;; FIXME: As long as the system dies on stack overflow or memory
160 ;;; exhaustion, it seems reasonable to have this, but its default
161 ;;; should be NIL, and when it's NIL, anything should be accepted.
162 (defparameter *intexp-maximum-exponent* 10000)
163
164 ;;; This function precisely calculates base raised to an integral
165 ;;; power. It separates the cases by the sign of power, for efficiency
166 ;;; reasons, as powers can be calculated more efficiently if power is
167 ;;; a positive integer. Values of power are calculated as positive
168 ;;; integers, and inverted if negative.
169 (defun intexp (base power)
170   (when (> (abs power) *intexp-maximum-exponent*)
171     ;; FIXME: should be ordinary error, not CERROR. (Once we set the
172     ;; default for the variable to NIL, the un-continuable error will
173     ;; be less obnoxious.)
174     (cerror "Continue with calculation."
175             "The absolute value of ~S exceeds ~S."
176             power '*intexp-maximum-exponent* base power))
177   (cond ((minusp power)
178          (/ (intexp base (- power))))
179         ((eql base 2)
180          (ash 1 power))
181         (t
182          (do ((nextn (ash power -1) (ash power -1))
183               (total (if (oddp power) base 1)
184                      (if (oddp power) (* base total) total)))
185              ((zerop nextn) total)
186            (setq base (* base base))
187            (setq power nextn)))))
188
189 ;;; If an integer power of a rational, use INTEXP above. Otherwise, do
190 ;;; floating point stuff. If both args are real, we try %POW right
191 ;;; off, assuming it will return 0 if the result may be complex. If
192 ;;; so, we call COMPLEX-POW which directly computes the complex
193 ;;; result. We also separate the complex-real and real-complex cases
194 ;;; from the general complex case.
195 (defun expt (base power)
196   #!+sb-doc
197   "Returns BASE raised to the POWER."
198   (if (zerop power)
199       (1+ (* base power))
200     (labels (;; determine if the double float is an integer.
201              ;;  0 - not an integer
202              ;;  1 - an odd int
203              ;;  2 - an even int
204              (isint (ihi lo)
205                (declare (type (unsigned-byte 31) ihi)
206                         (type (unsigned-byte 32) lo)
207                         (optimize (speed 3) (safety 0)))
208                (let ((isint 0))
209                  (declare (type fixnum isint))
210                  (cond ((>= ihi #x43400000)     ; exponent >= 53
211                         (setq isint 2))
212                        ((>= ihi #x3ff00000)
213                         (let ((k (- (ash ihi -20) #x3ff)))      ; exponent
214                           (declare (type (mod 53) k))
215                           (cond ((> k 20)
216                                  (let* ((shift (- 52 k))
217                                         (j (logand (ash lo (- shift))))
218                                         (j2 (ash j shift)))
219                                    (declare (type (mod 32) shift)
220                                             (type (unsigned-byte 32) j j2))
221                                    (when (= j2 lo)
222                                      (setq isint (- 2 (logand j 1))))))
223                                 ((= lo 0)
224                                  (let* ((shift (- 20 k))
225                                         (j (ash ihi (- shift)))
226                                         (j2 (ash j shift)))
227                                    (declare (type (mod 32) shift)
228                                             (type (unsigned-byte 31) j j2))
229                                    (when (= j2 ihi)
230                                      (setq isint (- 2 (logand j 1))))))))))
231                  isint))
232              (real-expt (x y rtype)
233                (let ((x (coerce x 'double-float))
234                      (y (coerce y 'double-float)))
235                  (declare (double-float x y))
236                  (let* ((x-hi (sb!kernel:double-float-high-bits x))
237                         (x-lo (sb!kernel:double-float-low-bits x))
238                         (x-ihi (logand x-hi #x7fffffff))
239                         (y-hi (sb!kernel:double-float-high-bits y))
240                         (y-lo (sb!kernel:double-float-low-bits y))
241                         (y-ihi (logand y-hi #x7fffffff)))
242                    (declare (type (signed-byte 32) x-hi y-hi)
243                             (type (unsigned-byte 31) x-ihi y-ihi)
244                             (type (unsigned-byte 32) x-lo y-lo))
245                    ;; y==zero: x**0 = 1
246                    (when (zerop (logior y-ihi y-lo))
247                      (return-from real-expt (coerce 1d0 rtype)))
248                    ;; +-NaN return x+y
249                    (when (or (> x-ihi #x7ff00000)
250                              (and (= x-ihi #x7ff00000) (/= x-lo 0))
251                              (> y-ihi #x7ff00000)
252                              (and (= y-ihi #x7ff00000) (/= y-lo 0)))
253                      (return-from real-expt (coerce (+ x y) rtype)))
254                    (let ((yisint (if (< x-hi 0) (isint y-ihi y-lo) 0)))
255                      (declare (type fixnum yisint))
256                      ;; special value of y
257                      (when (and (zerop y-lo) (= y-ihi #x7ff00000))
258                        ;; y is +-inf
259                        (return-from real-expt
260                          (cond ((and (= x-ihi #x3ff00000) (zerop x-lo))
261                                 ;; +-1**inf is NaN
262                                 (coerce (- y y) rtype))
263                                ((>= x-ihi #x3ff00000)
264                                 ;; (|x|>1)**+-inf = inf,0
265                                 (if (>= y-hi 0)
266                                     (coerce y rtype)
267                                     (coerce 0 rtype)))
268                                (t
269                                 ;; (|x|<1)**-,+inf = inf,0
270                                 (if (< y-hi 0)
271                                     (coerce (- y) rtype)
272                                     (coerce 0 rtype))))))
273
274                      (let ((abs-x (abs x)))
275                        (declare (double-float abs-x))
276                        ;; special value of x
277                        (when (and (zerop x-lo)
278                                   (or (= x-ihi #x7ff00000) (zerop x-ihi)
279                                       (= x-ihi #x3ff00000)))
280                          ;; x is +-0,+-inf,+-1
281                          (let ((z (if (< y-hi 0)
282                                       (/ 1 abs-x)       ; z = (1/|x|)
283                                       abs-x)))
284                            (declare (double-float z))
285                            (when (< x-hi 0)
286                              (cond ((and (= x-ihi #x3ff00000) (zerop yisint))
287                                     ;; (-1)**non-int
288                                     (let ((y*pi (* y pi)))
289                                       (declare (double-float y*pi))
290                                       (return-from real-expt
291                                         (complex
292                                          (coerce (%cos y*pi) rtype)
293                                          (coerce (%sin y*pi) rtype)))))
294                                    ((= yisint 1)
295                                     ;; (x<0)**odd = -(|x|**odd)
296                                     (setq z (- z)))))
297                            (return-from real-expt (coerce z rtype))))
298
299                        (if (>= x-hi 0)
300                            ;; x>0
301                            (coerce (sb!kernel::%pow x y) rtype)
302                            ;; x<0
303                            (let ((pow (sb!kernel::%pow abs-x y)))
304                              (declare (double-float pow))
305                              (case yisint
306                                (1 ; odd
307                                 (coerce (* -1d0 pow) rtype))
308                                (2 ; even
309                                 (coerce pow rtype))
310                                (t ; non-integer
311                                 (let ((y*pi (* y pi)))
312                                   (declare (double-float y*pi))
313                                   (complex
314                                    (coerce (* pow (%cos y*pi))
315                                            rtype)
316                                    (coerce (* pow (%sin y*pi))
317                                            rtype)))))))))))))
318       (declare (inline real-expt))
319       (number-dispatch ((base number) (power number))
320         (((foreach fixnum (or bignum ratio) (complex rational)) integer)
321          (intexp base power))
322         (((foreach single-float double-float) rational)
323          (real-expt base power '(dispatch-type base)))
324         (((foreach fixnum (or bignum ratio) single-float)
325           (foreach ratio single-float))
326          (real-expt base power 'single-float))
327         (((foreach fixnum (or bignum ratio) single-float double-float)
328           double-float)
329          (real-expt base power 'double-float))
330         ((double-float single-float)
331          (real-expt base power 'double-float))
332         (((foreach (complex rational) (complex float)) rational)
333          (* (expt (abs base) power)
334             (cis (* power (phase base)))))
335         (((foreach fixnum (or bignum ratio) single-float double-float)
336           complex)
337          (if (and (zerop base) (plusp (realpart power)))
338              (* base power)
339              (exp (* power (log base)))))
340         (((foreach (complex float) (complex rational))
341           (foreach complex double-float single-float))
342          (if (and (zerop base) (plusp (realpart power)))
343              (* base power)
344              (exp (* power (log base)))))))))
345
346 (defun log (number &optional (base nil base-p))
347   #!+sb-doc
348   "Return the logarithm of NUMBER in the base BASE, which defaults to e."
349   (if base-p
350       (if (zerop base)
351           base                          ; ANSI spec
352           (/ (log number) (log base)))
353       (number-dispatch ((number number))
354         (((foreach fixnum bignum ratio))
355          (if (minusp number)
356              (complex (log (- number)) (coerce pi 'single-float))
357              (coerce (%log (coerce number 'double-float)) 'single-float)))
358         (((foreach single-float double-float))
359          ;; Is (log -0) -infinity (libm.a) or -infinity + i*pi (Kahan)?
360          ;; Since this doesn't seem to be an implementation issue
361          ;; I (pw) take the Kahan result.
362          (if (< (float-sign number)
363                 (coerce 0 '(dispatch-type number)))
364              (complex (log (- number)) (coerce pi '(dispatch-type number)))
365              (coerce (%log (coerce number 'double-float))
366                      '(dispatch-type number))))
367         ((complex)
368          (complex-log number)))))
369
370 (defun sqrt (number)
371   #!+sb-doc
372   "Return the square root of NUMBER."
373   (number-dispatch ((number number))
374     (((foreach fixnum bignum ratio))
375      (if (minusp number)
376          (complex-sqrt number)
377          (coerce (%sqrt (coerce number 'double-float)) 'single-float)))
378     (((foreach single-float double-float))
379      (if (minusp number)
380          (complex-sqrt number)
381          (coerce (%sqrt (coerce number 'double-float))
382                  '(dispatch-type number))))
383      ((complex)
384       (complex-sqrt number))))
385 \f
386 ;;;; trigonometic and related functions
387
388 (defun abs (number)
389   #!+sb-doc
390   "Returns the absolute value of the number."
391   (number-dispatch ((number number))
392     (((foreach single-float double-float fixnum rational))
393      (abs number))
394     ((complex)
395      (let ((rx (realpart number))
396            (ix (imagpart number)))
397        (etypecase rx
398          (rational
399           (sqrt (+ (* rx rx) (* ix ix))))
400          (single-float
401           (coerce (%hypot (coerce rx 'double-float)
402                           (coerce ix 'double-float))
403                   'single-float))
404          (double-float
405           (%hypot rx ix)))))))
406
407 (defun phase (number)
408   #!+sb-doc
409   "Return the angle part of the polar representation of a complex number.
410   For complex numbers, this is (atan (imagpart number) (realpart number)).
411   For non-complex positive numbers, this is 0. For non-complex negative
412   numbers this is PI."
413   (etypecase number
414     (rational
415      (if (minusp number)
416          (coerce pi 'single-float)
417          0.0f0))
418     (single-float
419      (if (minusp (float-sign number))
420          (coerce pi 'single-float)
421          0.0f0))
422     (double-float
423      (if (minusp (float-sign number))
424          (coerce pi 'double-float)
425          0.0d0))
426     (complex
427      (atan (imagpart number) (realpart number)))))
428
429 (defun sin (number)
430   #!+sb-doc
431   "Return the sine of NUMBER."
432   (number-dispatch ((number number))
433     (handle-reals %sin number)
434     ((complex)
435      (let ((x (realpart number))
436            (y (imagpart number)))
437        (complex (* (sin x) (cosh y))
438                 (* (cos x) (sinh y)))))))
439
440 (defun cos (number)
441   #!+sb-doc
442   "Return the cosine of NUMBER."
443   (number-dispatch ((number number))
444     (handle-reals %cos number)
445     ((complex)
446      (let ((x (realpart number))
447            (y (imagpart number)))
448        (complex (* (cos x) (cosh y))
449                 (- (* (sin x) (sinh y))))))))
450
451 (defun tan (number)
452   #!+sb-doc
453   "Return the tangent of NUMBER."
454   (number-dispatch ((number number))
455     (handle-reals %tan number)
456     ((complex)
457      (complex-tan number))))
458
459 (defun cis (theta)
460   #!+sb-doc
461   "Return cos(Theta) + i sin(Theta), AKA exp(i Theta)."
462   (declare (type real theta))
463   (complex (cos theta) (sin theta)))
464
465 (defun asin (number)
466   #!+sb-doc
467   "Return the arc sine of NUMBER."
468   (number-dispatch ((number number))
469     ((rational)
470      (if (or (> number 1) (< number -1))
471          (complex-asin number)
472          (coerce (%asin (coerce number 'double-float)) 'single-float)))
473     (((foreach single-float double-float))
474      (if (or (> number (coerce 1 '(dispatch-type number)))
475              (< number (coerce -1 '(dispatch-type number))))
476          (complex-asin number)
477          (coerce (%asin (coerce number 'double-float))
478                  '(dispatch-type number))))
479     ((complex)
480      (complex-asin number))))
481
482 (defun acos (number)
483   #!+sb-doc
484   "Return the arc cosine of NUMBER."
485   (number-dispatch ((number number))
486     ((rational)
487      (if (or (> number 1) (< number -1))
488          (complex-acos number)
489          (coerce (%acos (coerce number 'double-float)) 'single-float)))
490     (((foreach single-float double-float))
491      (if (or (> number (coerce 1 '(dispatch-type number)))
492              (< number (coerce -1 '(dispatch-type number))))
493          (complex-acos number)
494          (coerce (%acos (coerce number 'double-float))
495                  '(dispatch-type number))))
496     ((complex)
497      (complex-acos number))))
498
499 (defun atan (y &optional (x nil xp))
500   #!+sb-doc
501   "Return the arc tangent of Y if X is omitted or Y/X if X is supplied."
502   (if xp
503       (flet ((atan2 (y x)
504                (declare (type double-float y x)
505                         (values double-float))
506                (if (zerop x)
507                    (if (zerop y)
508                        (if (plusp (float-sign x))
509                            y
510                            (float-sign y pi))
511                        (float-sign y (/ pi 2)))
512                    (%atan2 y x))))
513         (number-dispatch ((y number) (x number))
514           ((double-float
515             (foreach double-float single-float fixnum bignum ratio))
516            (atan2 y (coerce x 'double-float)))
517           (((foreach single-float fixnum bignum ratio)
518             double-float)
519            (atan2 (coerce y 'double-float) x))
520           (((foreach single-float fixnum bignum ratio)
521             (foreach single-float fixnum bignum ratio))
522            (coerce (atan2 (coerce y 'double-float) (coerce x 'double-float))
523                    'single-float))))
524       (number-dispatch ((y number))
525         (handle-reals %atan y)
526         ((complex)
527          (complex-atan y)))))
528
529 ;; It seems that everyone has a C version of sinh, cosh, and
530 ;; tanh. Let's use these for reals because the original
531 ;; implementations based on the definitions lose big in round-off
532 ;; error. These bad definitions also mean that sin and cos for
533 ;; complex numbers can also lose big.
534
535 #+nil
536 (defun sinh (number)
537   #!+sb-doc
538   "Return the hyperbolic sine of NUMBER."
539   (/ (- (exp number) (exp (- number))) 2))
540
541 (defun sinh (number)
542   #!+sb-doc
543   "Return the hyperbolic sine of NUMBER."
544   (number-dispatch ((number number))
545     (handle-reals %sinh number)
546     ((complex)
547      (let ((x (realpart number))
548            (y (imagpart number)))
549        (complex (* (sinh x) (cos y))
550                 (* (cosh x) (sin y)))))))
551
552 #+nil
553 (defun cosh (number)
554   #!+sb-doc
555   "Return the hyperbolic cosine of NUMBER."
556   (/ (+ (exp number) (exp (- number))) 2))
557
558 (defun cosh (number)
559   #!+sb-doc
560   "Return the hyperbolic cosine of NUMBER."
561   (number-dispatch ((number number))
562     (handle-reals %cosh number)
563     ((complex)
564      (let ((x (realpart number))
565            (y (imagpart number)))
566        (complex (* (cosh x) (cos y))
567                 (* (sinh x) (sin y)))))))
568
569 (defun tanh (number)
570   #!+sb-doc
571   "Return the hyperbolic tangent of NUMBER."
572   (number-dispatch ((number number))
573     (handle-reals %tanh number)
574     ((complex)
575      (complex-tanh number))))
576
577 (defun asinh (number)
578   #!+sb-doc
579   "Return the hyperbolic arc sine of NUMBER."
580   (number-dispatch ((number number))
581     (handle-reals %asinh number)
582     ((complex)
583      (complex-asinh number))))
584
585 (defun acosh (number)
586   #!+sb-doc
587   "Return the hyperbolic arc cosine of NUMBER."
588   (number-dispatch ((number number))
589     ((rational)
590      ;; acosh is complex if number < 1
591      (if (< number 1)
592          (complex-acosh number)
593          (coerce (%acosh (coerce number 'double-float)) 'single-float)))
594     (((foreach single-float double-float))
595      (if (< number (coerce 1 '(dispatch-type number)))
596          (complex-acosh number)
597          (coerce (%acosh (coerce number 'double-float))
598                  '(dispatch-type number))))
599     ((complex)
600      (complex-acosh number))))
601
602 (defun atanh (number)
603   #!+sb-doc
604   "Return the hyperbolic arc tangent of NUMBER."
605   (number-dispatch ((number number))
606     ((rational)
607      ;; atanh is complex if |number| > 1
608      (if (or (> number 1) (< number -1))
609          (complex-atanh number)
610          (coerce (%atanh (coerce number 'double-float)) 'single-float)))
611     (((foreach single-float double-float))
612      (if (or (> number (coerce 1 '(dispatch-type number)))
613              (< number (coerce -1 '(dispatch-type number))))
614          (complex-atanh number)
615          (coerce (%atanh (coerce number 'double-float))
616                  '(dispatch-type number))))
617     ((complex)
618      (complex-atanh number))))
619
620 ;;; HP-UX does not supply a C version of log1p, so
621 ;;; use the definition.
622
623 #!+hpux
624 #!-sb-fluid (declaim (inline %log1p))
625 #!+hpux
626 (defun %log1p (number)
627   (declare (double-float number)
628            (optimize (speed 3) (safety 0)))
629   (the double-float (log (the (double-float 0d0) (+ number 1d0)))))
630