# Diff of /src/code/irrat.lisp

revision 1.19.2.4 by pw, Tue May 23 16:36:34 2000 UTC revision 1.19.2.5 by pw, Sat Mar 23 18:50:03 2002 UTC
# Line 31  Line 31
31                                         "%"                                         "%"
32                                         (string-upcase name)))))                                         (string-upcase name)))))
33      `(progn      `(progn
34         (proclaim '(inline ,function))         (declaim (inline ,function))
35         (export ',function)         (export ',function)
36         (alien:def-alien-routine (,name ,function) double-float         (alien:def-alien-routine (,name ,function) double-float
37           ,@(let ((results nil))           ,@(let ((results nil))
# Line 524  Line 524
524          ((complex)          ((complex)
525           (complex-atan y)))))           (complex-atan y)))))
526
;; It seems that everyone has a C version of sinh, cosh, and
;; tanh. Let's use these for reals because the original
;; implementations based on the definitions lose big in round-off
;; error.  These bad definitions also mean that sin and cos for
;; complex numbers can also lose big.

#+nil
(defun sinh (number)
"Return the hyperbolic sine of NUMBER."
(/ (- (exp number) (exp (- number))) 2))

527  (defun sinh (number)  (defun sinh (number)
528    "Return the hyperbolic sine of NUMBER."    "Return the hyperbolic sine of NUMBER."
529    (number-dispatch ((number number))    (number-dispatch ((number number))
# Line 545  Line 534
534         (complex (* (sinh x) (cos y))         (complex (* (sinh x) (cos y))
535                  (* (cosh x) (sin y)))))))                  (* (cosh x) (sin y)))))))
536
#+nil
(defun cosh (number)
"Return the hyperbolic cosine of NUMBER."
(/ (+ (exp number) (exp (- number))) 2))

537  (defun cosh (number)  (defun cosh (number)
538    "Return the hyperbolic cosine of NUMBER."    "Return the hyperbolic cosine of NUMBER."
539    (number-dispatch ((number number))    (number-dispatch ((number number))
# Line 607  Line 591
591      ((complex)      ((complex)
592       (complex-atanh number))))       (complex-atanh number))))
593
594  ;;; HP-UX does not supply a C version of log1p, so  ;;; HP-UX does not supply a C version of log1p, so use the definition.
595  ;;; use the definition.  ;;; We really need to fix this.  The definition really loses big-time
596    ;;; in roundoff as x gets small.
597
598  #+hpux  #+hpux
599  (declaim (inline %log1p))  (declaim (inline %log1p))
# Line 619  Line 604
604    (the double-float (log (the (double-float 0d0) (+ number 1d0)))))    (the double-float (log (the (double-float 0d0) (+ number 1d0)))))
605
606
#+old-elfun
(progn
;;; Here are the old definitions of the special functions, for
;;; complex-valued arguments.  Some of these functions suffer from
;;; severe round-off error or unnecessary overflow.

(proclaim '(inline mult-by-i))
(defun mult-by-i (number)
(complex (- (imagpart number))
(realpart number)))

(defun complex-sqrt (x)
(exp (/ (log x) 2)))

(defun complex-log (x)
(complex (log (abs x))
(phase x)))

(defun complex-atanh (number)
(/ (- (log (1+ number)) (log (- 1 number))) 2))

(defun complex-tanh (number)
(/ (- (exp number) (exp (- number)))
(+ (exp number) (exp (- number)))))

(defun complex-acos (number)
(* -2 (mult-by-i (log (+ (sqrt (/ (1+ number) 2))
(mult-by-i (sqrt (/ (- 1 number) 2))))))))

(defun complex-acosh (number)
(* 2 (log (+ (sqrt (/ (1+ number) 2)) (sqrt (/ (1- number) 2))))))

(defun complex-asin (number)
(- (mult-by-i (log (+ (mult-by-i number) (sqrt (- 1 (* number number))))))))

(defun complex-asinh (number)
(log (+ number (sqrt (1+ (* number number))))))

(defun complex-atan (y)
(let ((im (imagpart y))
(re (realpart y)))
(/ (- (log (complex (- 1 im) re))
(log (complex (+ 1 im) (- re))))
(complex 0 2))))

(defun complex-tan (number)
(let* ((num (sin number))
(denom (cos number)))
(if (zerop denom) (error "~S undefined tangent." number)
(/ num denom))))
)

#-old-specfun
(progn
607  ;;;;  ;;;;
608  ;;;; This is a set of routines that implement many elementary  ;;;; This is a set of routines that implement many elementary
609  ;;;; transcendental functions as specified by ANSI Common Lisp.  The  ;;;; transcendental functions as specified by ANSI Common Lisp.  The
# Line 723  Line 654
654  ;;;; State of the Art in Numerical Analysis", pp. 165-211, Clarendon  ;;;; State of the Art in Numerical Analysis", pp. 165-211, Clarendon
655  ;;;; Press, 1987  ;;;; Press, 1987
656  ;;;;  ;;;;
657
658  (declaim (inline square))  (declaim (inline square))
(declaim (ftype (function (double-float) (double-float 0d0)) square))
659  (defun square (x)  (defun square (x)
660    (declare (double-float x)    (declare (double-float x))
(values (double-float 0d0)))
661    (* x x))    (* x x))
662
663  ;; If you have these functions in libm, perhaps they should be used  ;; If you have these functions in libm, perhaps they should be used
# Line 742  underlying floating-point format" Line 672  underlying floating-point format"
672             (type double-float-exponent n))             (type double-float-exponent n))
673    (scale-float x n))    (scale-float x n))
674
675    (declaim (inline logb-finite))
676    (defun logb-finite (x)
677      "Same as logb but X is not infinity and non-zero and not a NaN, so
678    that we can always return an integer"
679      (declare (type double-float x))
680      (multiple-value-bind (signif expon sign)
681          (decode-float x)
682        (declare (ignore signif sign))
683        ;; decode-float is almost right, except that the exponent
684        ;; is off by one
685        (1- expon)))
686
687  (defun logb (x)  (defun logb (x)
688    "Compute an integer N such that 1 <= |2^N * x| < 2.    "Compute an integer N such that 1 <= |2^(-N) * x| < 2.
689  For the special cases, the following values are used:  For the special cases, the following values are used:
690
691      x             logb      x             logb
# Line 758  For the special cases, the following val Line 700  For the special cases, the following val
700           #.ext:double-float-positive-infinity)           #.ext:double-float-positive-infinity)
701          ((zerop x)          ((zerop x)
702           ;; The answer is negative infinity, but we are supposed to           ;; The answer is negative infinity, but we are supposed to
703           ;; signal divide-by-zero.           ;; signal divide-by-zero, so do the actual division
;; (error 'division-by-zero :operation 'logb :operands (list x))
704           (/ -1.0d0 x)           (/ -1.0d0 x)
705           )           )
706          (t          (t
707           (multiple-value-bind (signif expon sign)           (logb-finite x))))
708               (decode-float x)
709             (declare (ignore signif sign))
;; decode-float is almost right, except that the exponent
;; is off by one
(1- expon)))))
710
711  ;; This function is used to create a complex number of the appropriate  ;; This function is used to create a complex number of the appropriate
712  ;; type.  ;; type.
# Line 779  For the special cases, the following val Line 717  For the special cases, the following val
717  it has the same type as Z.  If Z has type (complex rational), the X  it has the same type as Z.  If Z has type (complex rational), the X
718  and Y are coerced to single-float."  and Y are coerced to single-float."
719    (declare (double-float x y)    (declare (double-float x y)
720             (number z))             (number z)
721    (if (subtypep (type-of (realpart z)) 'double-float)             (optimize (extensions:inhibit-warnings 3)))
722      (if (typep (realpart z) 'double-float)
723        (complex x y)        (complex x y)
724        ;; Convert anything that's not a double-float to a single-float.        ;; Convert anything that's not a double-float to a single-float.
725        (complex (float x 1.0)        (complex (float x 1f0)
726                 (float y 1.0))))                 (float y 1f0))))
727
728  (defun cssqs (z)  (defun cssqs (z)
729    ;; Compute |(x+i*y)/2^k|^2 scaled to avoid over/underflow. The    ;; Compute |(x+i*y)/2^k|^2 scaled to avoid over/underflow. The
# Line 792  and Y are coerced to single-float." Line 731  and Y are coerced to single-float."
731
732    ;; Save all FP flags    ;; Save all FP flags
733    (let ((x (float (realpart z) 1d0))    (let ((x (float (realpart z) 1d0))
734          (y (float (imagpart z) 1d0))          (y (float (imagpart z) 1d0)))
(k 0)
(rho 0d0))
(declare (double-float x y)
(type (double-float 0d0) rho)
(fixnum k))
735      ;; Would this be better handled using an exception handler to      ;; Would this be better handled using an exception handler to
736      ;; catch the overflow or underflow signal?  For now, we turn all      ;; catch the overflow or underflow signal?  For now, we turn all
737      ;; traps off and look at the accrued exceptions to see if any      ;; traps off and look at the accrued exceptions to see if any
738      ;; signal would have been raised.      ;; signal would have been raised.
740        (setf rho (+ (square x) (square y)))        (let ((rho (+ (square x) (square y))))
741        (cond ((and (or (float-nan-p rho)          (declare (optimize (speed 3) (space 0)))
742                        (float-infinity-p rho))          (cond ((and (or (float-nan-p rho)
743                    (or (float-infinity-p (abs x))                          (float-infinity-p rho))
744                        (float-infinity-p (abs y))))                      (or (float-infinity-p (abs x))
745               (setf rho #.ext:double-float-positive-infinity))                          (float-infinity-p (abs y))))
746              ((let ((threshold #.(/ least-positive-double-float                 (values ext:double-float-positive-infinity 0))
747                                     double-float-epsilon))                ((let ((threshold #.(/ least-positive-double-float
748                     (traps (ldb vm::float-sticky-bits                                       double-float-epsilon))
749                                 (vm:floating-point-modes))))                       (traps (ldb vm::float-sticky-bits
750                 ;; Overflow raised or (underflow raised and rho <                                   (vm:floating-point-modes))))
751                 ;; lambda/eps)                   ;; Overflow raised or (underflow raised and rho <
752                 (or (not (zerop (logand vm:float-overflow-trap-bit traps)))                   ;; lambda/eps)
753                     (and (not (zerop (logand vm:float-underflow-trap-bit traps)))                   (or (not (zerop (logand vm:float-overflow-trap-bit traps)))
754                          (< rho threshold))))                       (and (not (zerop (logand vm:float-underflow-trap-bit traps)))
755               (setf k (logb (max (abs x) (abs y))))                            (< rho threshold))))
756               (setf rho (+ (square (scalb x (- k)))                 ;; If we're here, neither x nor y are infinity and at
757                            (square (scalb y (- k))))))))                 ;; least one is non-zero.. Thus logb returns a nice
758      (values rho k)))                 ;; integer.
759                   (let ((k (- (logb-finite (max (abs x) (abs y))))))
760                     (values (+ (square (scalb x k))
761                                (square (scalb y k)))
762                             (- k))))
763                  (t
764                   (values rho 0)))))))
765
766  (defun complex-sqrt (z)  (defun complex-sqrt (z)
767    "Principle square root of Z    "Principle square root of Z
# Line 830  Z may be any number, but the result is a Line 770  Z may be any number, but the result is a
770    (declare (number z))    (declare (number z))
771    (multiple-value-bind (rho k)    (multiple-value-bind (rho k)
772        (cssqs z)        (cssqs z)
773      (declare (type (double-float 0d0) rho)      (declare (type (or (member 0d0) (double-float 0d0)) rho)
774               (fixnum k))               (type fixnum k))
775      (let ((x (float (realpart z) 1.0d0))      (let ((x (float (realpart z) 1.0d0))
776            (y (float (imagpart z) 1.0d0))            (y (float (imagpart z) 1.0d0))
777            (eta 0d0)            (eta 0d0)
778            (nu 0d0))            (nu 0d0))
779        (declare (double-float x y eta nu))        (declare (double-float x y eta nu))
780
781        (if (not (float-nan-p x))        (locally
782            (setf rho (+ (scalb (abs x) (- k)) (sqrt rho))))            ;; space 0 to get maybe-inline functions inlined.
783              (declare (optimize (speed 3) (space 0)))
784        (cond ((oddp k)
785               (setf k (ash k -1)))          (if (not (float-nan-p x))
786              (t              (setf rho (+ (scalb (abs x) (- k)) (sqrt rho))))
787               (setf k (1- (ash k -1)))
788               (setf rho (+ rho rho))))          (cond ((oddp k)
789                   (setf k (ash k -1)))
790        (setf rho (scalb (sqrt rho) k))                (t
791                   (setf k (1- (ash k -1)))
792                   (setf rho (+ rho rho))))
793
794            (setf rho (scalb (sqrt rho) k))
795
796            (setf eta rho)
797            (setf nu y)
798
799            (when (/= rho 0d0)
800              (when (not (float-infinity-p (abs nu)))
801                (setf nu (/ (/ nu rho) 2d0)))
802              (when (< x 0d0)
803                (setf eta (abs nu))
804                (setf nu (float-sign y rho))))
805            (coerce-to-complex-type eta nu z)))))
806
(setf eta rho)
(setf nu y)

(when (/= rho 0d0)
(when (not (float-infinity-p (abs nu)))
(setf nu (/ (/ nu rho) 2d0)))
(when (< x 0d0)
(setf eta (abs nu))
(setf nu (float-sign y rho))))
(coerce-to-complex-type eta nu z))))

807  (defun complex-log-scaled (z j)  (defun complex-log-scaled (z j)
808    "Compute log(2^j*z).    "Compute log(2^j*z).
809
# Line 879  This is for use with J /= 0 only when |z Line 823  This is for use with J /= 0 only when |z
823          (y (float (imagpart z) 1.0d0)))          (y (float (imagpart z) 1.0d0)))
824      (multiple-value-bind (rho k)      (multiple-value-bind (rho k)
825          (cssqs z)          (cssqs z)
826        (declare (type (double-float 0d0) rho)        (declare (optimize (speed 3)))
(fixnum k))
827        (let ((beta (max (abs x) (abs y)))        (let ((beta (max (abs x) (abs y)))
828              (theta (min (abs x) (abs y))))              (theta (min (abs x) (abs y))))
829          (declare (type (double-float 0d0) beta theta))          (coerce-to-complex-type (if (and (zerop k)
830          (if (and (zerop k)                                           (< t0 beta)
831                   (< t0 beta)                                           (or (<= beta t1)
832                   (or (<= beta t1)                                               (< rho t2)))
833                       (< rho t2)))                                      (/ (%log1p (+ (* (- beta 1.0d0)
834              (setf rho (/ (%log1p (+ (* (- beta 1.0d0)                                                       (+ beta 1.0d0))
835                                         (+ beta 1.0d0))                                                    (* theta theta)))
836                                      (* theta theta)))                                         2d0)
837                           2d0))                                      (+ (/ (log rho) 2d0)
838              (setf rho (+ (/ (log rho) 2d0)                                         (* (+ k j) ln2)))
839                           (* (+ k j) ln2))))                                  (atan y x)
840          (setf theta (atan y x))                                  z)))))
(coerce-to-complex-type rho theta z)))))
841
842  (defun complex-log (z)  (defun complex-log (z)
843    "Log of Z = log |Z| + i * arg Z    "Log of Z = log |Z| + i * arg Z
# Line 913  Z may be any number, but the result is a Line 855  Z may be any number, but the result is a
855    "Compute atanh z = (log(1+z) - log(1-z))/2"    "Compute atanh z = (log(1+z) - log(1-z))/2"
856    (declare (number z))    (declare (number z))
857    (let* (;; Constants    (let* (;; Constants
858           (theta #.(/ (sqrt most-positive-double-float) 4.0d0))           (theta (/ (sqrt most-positive-double-float) 4.0d0))
859           (rho #.(/ 4.0d0 (sqrt most-positive-double-float)))           (rho (/ 4.0d0 (sqrt most-positive-double-float)))
860           (half-pi #.(/ pi 2.0d0))           (half-pi (/ pi 2.0d0))
861           (rp (float (realpart z) 1.0d0))           (rp (float (realpart z) 1.0d0))
862           (beta (float-sign rp 1.0d0))           (beta (float-sign rp 1.0d0))
863           (x (* beta rp))           (x (* beta rp))
864           (y (* beta (- (float (imagpart z) 1.0d0))))           (y (* beta (- (float (imagpart z) 1.0d0))))
865           (eta 0.0d0)           (eta 0.0d0)
866           (nu 0.0d0))           (nu 0.0d0))
867      (declare (double-float theta rho half-pi rp beta y eta nu)      ;; Shouldn't need this declare.
868               (type (double-float 0d0) x))      (declare (double-float x y))
869      (cond ((or (> x theta)      (locally
870                 (> (abs y) theta))          (declare (optimize (speed 3)))
871             ;; To avoid overflow...        (cond ((or (> x theta)
872             (setf eta (float-sign y half-pi))                   (> (abs y) theta))
873             ;; nu is real part of 1/(x + iy).  This is x/(x^2+y^2),               ;; To avoid overflow...
874             ;; which can cause overflow.  Arrange this computation so               (setf eta (float-sign y half-pi))
875             ;; that it won't overflow.               ;; nu is real part of 1/(x + iy).  This is x/(x^2+y^2),
876             (setf nu (let* ((x-bigger (> x (abs y)))               ;; which can cause overflow.  Arrange this computation so
877                             (r (if x-bigger (/ y x) (/ x y)))               ;; that it won't overflow.
878                             (d (+ 1.0d0 (* r r))))               (setf nu (let* ((x-bigger (> x (abs y)))
879                        (declare (double-float r d))                               (r (if x-bigger (/ y x) (/ x y)))
880                        (if x-bigger                               (d (+ 1.0d0 (* r r))))
881                            (/ (/ x) d)                          (if x-bigger
882                            (/ (/ r y) d)))))                              (/ (/ x) d)
883            ((= x 1.0d0)                              (/ (/ r y) d)))))
884             ;; Should this be changed so that if y is zero, eta is set              ((= x 1.0d0)
885             ;; to +infinity instead of approx 176?  In any case               ;; Should this be changed so that if y is zero, eta is set
886             ;; tanh(176) is 1.0d0 within working precision.               ;; to +infinity instead of approx 176?  In any case
887             (let ((t1 (+ 4d0 (square y)))               ;; tanh(176) is 1.0d0 within working precision.
888                   (t2 (+ (abs y) rho)))               (let ((t1 (+ 4d0 (square y)))
889               (declare (type (double-float 0d0) t1 t2))                     (t2 (+ (abs y) rho)))
890               #+nil                 (setf eta (log (/ (sqrt (sqrt t1)))
891               (setf eta (log (/ (sqrt (sqrt t1)))                                (sqrt t2)))
892                              (sqrt t2)))                 (setf nu (* 0.5d0
893               (setf eta (* 0.5d0 (log (the (double-float 0.0d0)                             (float-sign y
894                                            (/ (sqrt t1) t2)))))                                         (+ half-pi (atan (* 0.5d0 t2))))))))
895               (setf nu (* 0.5d0              (t
896                           (float-sign y               (let ((t1 (+ (abs y) rho)))
897                                       (+ half-pi (atan (* 0.5d0 t2))))))))                 ;; Normal case using log1p(x) = log(1 + x)
898            (t                 (setf eta (* 0.25d0
899             (let ((t1 (+ (abs y) rho)))                              (%log1p (/ (* 4.0d0 x)
900               (declare (double-float t1))                                         (+ (square (- 1.0d0 x))
901               ;; Normal case using log1p(x) = log(1 + x)                                            (square t1))))))
902               (setf eta (* 0.25d0                 (setf nu (* 0.5d0
903                            (%log1p (/ (* 4.0d0 x)                             (atan (* 2.0d0 y)
904                                       (+ (square (- 1.0d0 x))                                   (- (* (- 1.0d0 x)
905                                          (square t1))))))                                         (+ 1.0d0 x))
906               (setf nu (* 0.5d0                                      (square t1))))))))
907                           (atan (* 2.0d0 y)        (coerce-to-complex-type (* beta eta)
908                                 (- (* (- 1.0d0 x)                                (- (* beta nu))
909                                       (+ 1.0d0 x))                                z))))
(square t1))))))))
(coerce-to-complex-type (* beta eta)
(- (* beta nu))
z)))
910
911  (defun complex-tanh (z)  (defun complex-tanh (z)
912    "Compute tanh z = sinh z / cosh z"    "Compute tanh z = sinh z / cosh z"
913    (declare (number z))    (declare (number z))
914    (let ((x (float (realpart z) 1.0d0))    (let ((x (float (realpart z) 1.0d0))
915          (y (float (imagpart z) 1.0d0)))          (y (float (imagpart z) 1.0d0)))
916      (declare (double-float x y))      (locally
917      (cond ((> (abs x)          ;; space 0 to get maybe-inline functions inlined
918                #-(or linux hpux) #.(/ (%asinh most-positive-double-float) 4d0)          (declare (optimize (speed 3) (space 0)))
919                ;; This is more accurate under linux.        (cond ((> (abs x)
920                #+(or linux hpux) #.(/ (+ (%log 2.0d0)                  #-(or linux hpux) #.(/ (%asinh most-positive-double-float) 4d0)
921                                          (%log most-positive-double-float)) 4d0))                  ;; This is more accurate under linux.
922             (complex (float-sign x)                  #+(or linux hpux) #.(/ (+ (%log 2.0d0)
923                      (float-sign y 0.0d0)))                                            (%log most-positive-double-float)) 4d0))
924            (t               (coerce-to-complex-type (float-sign x)
925             (let* ((tv (%tan y))                                       (float-sign y) z))
926                    (beta (+ 1.0d0 (* tv tv)))              (t
927                    (s (sinh x))               (let* ((tv (%tan y))
928                    (rho (sqrt (+ 1.0d0 (* s s)))))                      (beta (+ 1.0d0 (* tv tv)))
929               (declare (double-float tv s)                      (s (sinh x))
930                        (type (double-float 0.0d0) beta rho))                      (rho (sqrt (+ 1.0d0 (* s s)))))
931               (if (float-infinity-p (abs tv))                 (if (float-infinity-p (abs tv))
932                   (coerce-to-complex-type (/ rho s)                     (coerce-to-complex-type (/ rho s)
933                                           (/ tv)                                             (/ tv)
934                                           z)                                             z)
935                   (let ((den (+ 1.0d0 (* beta s s))))                     (let ((den (+ 1.0d0 (* beta s s))))
936                     (coerce-to-complex-type (/ (* beta rho s)                       (coerce-to-complex-type (/ (* beta rho s)
937                                                den)                                                  den)
938                                             (/ tv den)                                               (/ tv den)
939                                             z))))))))                                               z)))))))))
940
941  ;; Kahan says we should only compute the parts needed.  Thus, the  ;; Kahan says we should only compute the parts needed.  Thus, the
942  ;; realpart's below should only compute the real part, not the whole  ;; realpart's below should only compute the real part, not the whole
# Line 1094  Z may be any number, but the result is a Line 1032  Z may be any number, but the result is a
1032           (result (complex-tanh iz)))           (result (complex-tanh iz)))
1033      (complex (imagpart result)      (complex (imagpart result)
1034               (- (realpart result)))))               (- (realpart result)))))
1035  )

Legend:
 Removed from v.1.19.2.4 changed lines Added in v.1.19.2.5