X-Git-Url: http://repo.macrolet.net/gitweb/?a=blobdiff_plain;f=src%2Fcode%2Firrat.lisp;h=86940c3f41e02c06996774893e8208010772862c;hb=0e3c4b4db102bd204a30402d7e5a0de44aea57ce;hp=1a38ed99a6cafedc6f76de7fae6900dbf99cfd91;hpb=7fb597b585fc715537ea644f7d84440eca217ca1;p=sbcl.git

diff --git a/src/code/irrat.lisp b/src/code/irrat.lisp
index 1a38ed9..86940c3 100644
--- a/src/code/irrat.lisp
+++ b/src/code/irrat.lisp
@@ -23,15 +23,18 @@
 (eval-when (:compile-toplevel :execute)
 
 (sb!xc:defmacro def-math-rtn (name num-args)
-  (let ((function (symbolicate "%" (string-upcase name))))
+  (let ((function (symbolicate "%" (string-upcase name)))
+        (args (loop for i below num-args
+                    collect (intern (format nil "ARG~D" i)))))
     `(progn
        (declaim (inline ,function))
-       (sb!alien:define-alien-routine (,name ,function) double-float
-         ,@(let ((results nil))
-             (dotimes (i num-args (nreverse results))
-               (push (list (intern (format nil "ARG-~D" i))
-                           'double-float)
-                     results)))))))
+       (defun ,function ,args
+         (alien-funcall
+          (extern-alien ,name
+                        (function double-float
+                                  ,@(loop repeat num-args
+                                          collect 'double-float)))
+          ,@args)))))
 
 (defun handle-reals (function var)
   `((((foreach fixnum single-float bignum ratio))
@@ -46,8 +49,11 @@
              `(defun ,name ,ll (,name ,@ll))))
   (def %atan2 (x y))
   (def %atan (x))
+  (def %tan (x))
   (def %tan-quick (x))
+  (def %cos (x))
   (def %cos-quick (x))
+  (def %sin (x))
   (def %sin-quick (x))
   (def %sqrt (x))
   (def %log (x))
@@ -67,25 +73,84 @@
 #!-x86 (def-math-rtn "sin" 1)
 #!-x86 (def-math-rtn "cos" 1)
 #!-x86 (def-math-rtn "tan" 1)
-(def-math-rtn "asin" 1)
-(def-math-rtn "acos" 1)
 #!-x86 (def-math-rtn "atan" 1)
 #!-x86 (def-math-rtn "atan2" 2)
-(def-math-rtn "sinh" 1)
-(def-math-rtn "cosh" 1)
-#!-win32(def-math-rtn "tanh" 1)
-#!-win32(def-math-rtn "asinh" 1)
-#!-win32(def-math-rtn "acosh" 1)
-#!-win32(def-math-rtn "atanh" 1)
+#!-(and win32 x86)
+(progn
+  (def-math-rtn "acos" 1)
+  (def-math-rtn "asin" 1)
+  (def-math-rtn "cosh" 1)
+  (def-math-rtn "sinh" 1)
+  (def-math-rtn "tanh" 1)
+  #!-win32
+  (progn
+    (def-math-rtn "asinh" 1)
+    (def-math-rtn "acosh" 1)
+    (def-math-rtn "atanh" 1)))
+#!+win32
+(progn
+  #!-x86-64
+  (progn
+    (declaim (inline %asin))
+    (defun %asin (number)
+      (%atan (/ number (sqrt (- 1 (* number number))))))
+    (declaim (inline %acos))
+    (defun %acos (number)
+      (- (/ pi 2) (%asin number)))
+    (declaim (inline %cosh))
+    (defun %cosh (number)
+      (/ (+ (exp number) (exp (- number))) 2))
+    (declaim (inline %sinh))
+    (defun %sinh (number)
+      (/ (- (exp number) (exp (- number))) 2))
+    (declaim (inline %tanh))
+    (defun %tanh (number)
+      (/ (%sinh number) (%cosh number))))
+  (declaim (inline %asinh))
+  (defun %asinh (number)
+    (log (+ number (sqrt (+ (* number number) 1.0d0))) #.(exp 1.0d0)))
+  (declaim (inline %acosh))
+  (defun %acosh (number)
+    (log (+ number (sqrt (- (* number number) 1.0d0))) #.(exp 1.0d0)))
+  (declaim (inline %atanh))
+  (defun %atanh (number)
+    (let ((ratio (/ (+ 1 number) (- 1 number))))
+      ;; Were we effectively zero?
+      (if (= ratio -1.0d0)
+          0.0d0
+          (/ (log ratio #.(exp 1.0d0)) 2.0d0)))))
 
 ;;; exponential and logarithmic
 #!-x86 (def-math-rtn "exp" 1)
 #!-x86 (def-math-rtn "log" 1)
 #!-x86 (def-math-rtn "log10" 1)
-#!-win32(def-math-rtn "pow" 2)
+#!-(and win32 x86) (def-math-rtn "pow" 2)
 #!-(or x86 x86-64) (def-math-rtn "sqrt" 1)
-(def-math-rtn "hypot" 2)
-#!-(or hpux x86) (def-math-rtn "log1p" 1)
+#!-win32 (def-math-rtn "hypot" 2)
+#!-x86 (def-math-rtn "log1p" 1)
+
+#!+win32
+(progn
+  ;; This is written in a peculiar way to avoid overflow. Note that in
+  ;; sqrt(x^2 + y^2), either square or the sum can overflow.
+  ;;
+  ;; Factoring x^2 out of sqrt(x^2 + y^2) gives us the expression
+  ;; |x|sqrt(1 + (y/x)^2), which, assuming |x| >= |y|, can only overflow
+  ;; if |x| is sufficiently large.
+  ;;
+  ;; The ZEROP test suffices (y is non-negative) to guard against
+  ;; divisions by zero: x >= y > 0.
+  (declaim (inline %hypot))
+  (defun %hypot (x y)
+    (declare (type double-float x y))
+    (let ((x (abs x))
+          (y (abs y)))
+      (when (> y x)
+        (rotatef x y))
+      (if (zerop y)
+          x
+          (let ((y/x (/ y x)))
+            (* x (sqrt (1+ (* y/x y/x)))))))))
 
 ;;;; power functions
 
@@ -135,10 +200,15 @@
   #!+sb-doc
   "Return BASE raised to the POWER."
   (if (zerop power)
-      (let ((result (1+ (* base power))))
-        (if (and (floatp result) (float-nan-p result))
-            (float 1 result)
-            result))
+    (if (and (zerop base) (floatp power))
+        (error 'arguments-out-of-domain-error
+               :operands (list base power)
+               :operation 'expt
+               :references (list '(:ansi-cl :function expt)))
+        (let ((result (1+ (* base power))))
+          (if (and (floatp result) (float-nan-p result))
+              (float 1 result)
+              result)))
     (labels (;; determine if the double float is an integer.
              ;;  0 - not an integer
              ;;  1 - an odd int
@@ -188,6 +258,7 @@
                    (when (zerop (logior y-ihi y-lo))
                      (return-from real-expt (coerce 1d0 rtype)))
                    ;; +-NaN return x+y
+                   ;; FIXME: Hardcoded qNaN/sNaN values are not portable.
                    (when (or (> x-ihi #x7ff00000)
                              (and (= x-ihi #x7ff00000) (/= x-lo 0))
                              (> y-ihi #x7ff00000)
@@ -256,8 +327,12 @@
                                    (coerce (* pow (%cos y*pi))
                                            rtype)
                                    (coerce (* pow (%sin y*pi))
-                                           rtype)))))))))))))
-      (declare (inline real-expt))
+                                           rtype))))))))))))
+             (complex-expt (base power)
+               (if (and (zerop base) (plusp (realpart power)))
+                   (* base power)
+                   (exp (* power (log base))))))
+      (declare (inline real-expt complex-expt))
       (number-dispatch ((base number) (power number))
         (((foreach fixnum (or bignum ratio) (complex rational)) integer)
          (intexp base power))
@@ -271,19 +346,33 @@
          (real-expt base power 'double-float))
         ((double-float single-float)
          (real-expt base power 'double-float))
-        (((foreach (complex rational) (complex float)) rational)
+        ;; Handle (expt <complex> <rational>), except the case dealt with
+        ;; in the first clause above, (expt <(complex rational)> <integer>).
+        (((foreach (complex rational) (complex single-float)
+                   (complex double-float))
+          rational)
          (* (expt (abs base) power)
             (cis (* power (phase base)))))
-        (((foreach fixnum (or bignum ratio) single-float double-float)
-          complex)
-         (if (and (zerop base) (plusp (realpart power)))
-             (* base power)
-             (exp (* power (log base)))))
-        (((foreach (complex float) (complex rational))
+        ;; The next three clauses handle (expt <real> <complex>).
+        (((foreach fixnum (or bignum ratio) single-float)
+          (foreach (complex single-float) (complex rational)))
+         (complex-expt base power))
+        (((foreach fixnum (or bignum ratio) single-float)
+          (complex double-float))
+         (complex-expt (coerce base 'double-float) power))
+        ((double-float complex)
+         (complex-expt base power))
+        ;; The next three clauses handle (expt <complex> <float>) and
+        ;; (expt <complex> <complex>).
+        (((foreach (complex single-float) (complex rational))
+          (foreach (complex single-float) (complex rational) single-float))
+         (complex-expt base power))
+        (((foreach (complex single-float) (complex rational))
+          (foreach (complex double-float) double-float))
+         (complex-expt (coerce base '(complex double-float)) power))
+        (((complex double-float)
           (foreach complex double-float single-float))
-         (if (and (zerop base) (plusp (realpart power)))
-             (* base power)
-             (exp (* power (log base)))))))))
+         (complex-expt base power))))))
 
 ;;; FIXME: Maybe rename this so that it's clearer that it only works
 ;;; on integers?
@@ -311,11 +400,20 @@
   "Return the logarithm of NUMBER in the base BASE, which defaults to e."
   (if base-p
       (cond
-        ((zerop base) 0f0) ; FIXME: type
+        ((zerop base)
+         (if (or (typep number 'double-float) (typep base 'double-float))
+             0.0d0
+             0.0f0))
         ((and (typep number '(integer (0) *))
               (typep base '(integer (0) *)))
          (coerce (/ (log2 number) (log2 base)) 'single-float))
-        (t (/ (log number) (log base))))
+        ((and (typep number 'integer) (typep base 'double-float))
+         ;; No single float intermediate result
+         (/ (log2 number) (log base 2.0d0)))
+        ((and (typep number 'double-float) (typep base 'integer))
+         (/ (log number 2.0d0) (log2 base)))
+        (t
+         (/ (log number) (log base))))
       (number-dispatch ((number number))
         (((foreach fixnum bignum))
          (if (minusp number)
@@ -583,18 +681,6 @@
     ((complex)
      (complex-atanh number))))
 
-;;; HP-UX does not supply a C version of log1p, so use the definition.
-;;;
-;;; FIXME: This is really not a good definition. As per Raymond Toy
-;;; working on CMU CL, "The definition really loses big-time in
-;;; roundoff as x gets small."
-#!+hpux
-#!-sb-fluid (declaim (inline %log1p))
-#!+hpux
-(defun %log1p (number)
-  (declare (double-float number)
-           (optimize (speed 3) (safety 0)))
-  (the double-float (log (the (double-float 0d0) (+ number 1d0)))))
 
 ;;;; not-OLD-SPECFUN stuff
 ;;;;
@@ -767,8 +853,13 @@
              (values
               (double-from-bits 0 (1+ sb!vm:double-float-normal-exponent-max) 0)
               0))
-            ((let ((threshold #.(/ least-positive-double-float
-                                   double-float-epsilon))
+            ((let ((threshold
+                    ;; (/ least-positive-double-float double-float-epsilon)
+                    (load-time-value
+                     #!-long-float
+                     (sb!kernel:make-double-float #x1fffff #xfffffffe)
+                     #!+long-float
+                     (error "(/ least-positive-long-float long-float-epsilon)")))
                    (traps (ldb sb!vm::float-sticky-bits
                                (sb!vm:floating-point-modes))))
                 ;; Overflow raised or (underflow raised and rho <
@@ -845,10 +936,22 @@
   ;; influences the choices of these constants but doesn't say how to
   ;; choose them.  We'll just assume his choices matches our
   ;; implementation of log1p.
-  (let ((t0 #.(/ 1 (sqrt 2.0d0)))
+  (let ((t0 (load-time-value
+             #!-long-float
+             (sb!kernel:make-double-float #x3fe6a09e #x667f3bcd)
+             #!+long-float
+             (error "(/ (sqrt 2l0))")))
+        ;; KLUDGE: if repeatable fasls start failing under some weird
+        ;; xc host, this 1.2d0 might be a good place to examine: while
+        ;; it _should_ be the same in all vaguely-IEEE754 hosts, 1.2
+        ;; is not exactly representable, so something could go wrong.
         (t1 1.2d0)
         (t2 3d0)
-        (ln2 #.(log 2d0))
+        (ln2 (load-time-value
+              #!-long-float
+              (sb!kernel:make-double-float #x3fe62e42 #xfefa39ef)
+              #!+long-float
+              (error "(log 2l0)")))
         (x (float (realpart z) 1.0d0))
         (y (float (imagpart z) 1.0d0)))
     (multiple-value-bind (rho k)
@@ -946,16 +1049,11 @@
       ;; space 0 to get maybe-inline functions inlined
       (declare (optimize (speed 3) (space 0)))
     (cond ((> (abs x)
-              ;; FIXME: this form is hideously broken wrt
-              ;; cross-compilation portability.  Much else in this
-              ;; file is too, of course, sometimes hidden by
-              ;; constant-folding, but this one in particular clearly
-              ;; depends on host and target
-              ;; MOST-POSITIVE-DOUBLE-FLOATs being equal.  -- CSR,
-              ;; 2003-04-20
-              #.(/ (+ (log 2.0d0)
-                      (log most-positive-double-float))
-                   4d0))
+              (load-time-value
+               #!-long-float
+               (sb!kernel:make-double-float #x406633ce #x8fb9f87e)
+               #!+long-float
+               (error "(/ (+ (log 2l0) (log most-positive-long-float)) 4l0)")))
            (coerce-to-complex-type (float-sign x)
                                    (float-sign y) z))
           (t