X-Git-Url: http://repo.macrolet.net/gitweb/?a=blobdiff_plain;f=src%2Fcode%2Firrat.lisp;h=5e2c8d0bc5ec37c63b435cc3130be3cf79b64342;hb=8eee0d3a30bf39d9f201acff28c92059fe6c3e4e;hp=0f2b7a321a424d09d40a6c5efdbbe79e48209423;hpb=6349f899f37ca9807461e67cb5e2f8b8bd3beba6;p=sbcl.git

diff --git a/src/code/irrat.lisp b/src/code/irrat.lisp
index 0f2b7a3..5e2c8d0 100644
--- a/src/code/irrat.lisp
+++ b/src/code/irrat.lisp
@@ -70,16 +70,48 @@
 #!-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)
+#!-win32
+(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)
+  (def-math-rtn "asinh" 1)
+  (def-math-rtn "acosh" 1)
+  (def-math-rtn "atanh" 1))
+#!+win32
+(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)
@@ -87,8 +119,16 @@
 #!-x86 (def-math-rtn "log10" 1)
 #!-win32(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
+  ;; FIXME: libc hypot "computes the sqrt(x*x+y*y) without undue overflow or underflow"
+  ;; ...we just do the stupid simple thing.
+  (declaim (inline %hypot))
+  (defun %hypot (x y)
+    (sqrt (+ (* x x) (* y y)))))
 
 ;;;; power functions
 
@@ -191,6 +231,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)
@@ -314,11 +355,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)
@@ -586,18 +636,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
 ;;;;