/[cmucl]/src/code/irrat-dd.lisp
ViewVC logotype

Diff of /src/code/irrat-dd.lisp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.16 by rtoy, Wed Jun 27 16:38:54 2007 UTC revision 1.17 by rtoy, Thu Aug 2 18:18:18 2007 UTC
# Line 145  Line 145 
145    (defun dd-%expm1 (x)    (defun dd-%expm1 (x)
146      "exp(x) - 1"      "exp(x) - 1"
147      (declare (type double-double-float x)      (declare (type double-double-float x)
148               (optimize (speed 3) (space 0)))               (optimize (speed 3) (space 0)
149                           (inhibit-warnings 3)))
150      ;; Range reduction is accomplished by separating the argument      ;; Range reduction is accomplished by separating the argument
151      ;; into an integer k and fraction f such that      ;; into an integer k and fraction f such that
152      ;;      ;;
# Line 242  Line 243 
243    ;; relative error = 2.2e-37, relative peak error spread = 9.2e-38.    ;; relative error = 2.2e-37, relative peak error spread = 9.2e-38.
244    (defun dd-%exp (x)    (defun dd-%exp (x)
245      (declare (type double-double-float x)      (declare (type double-double-float x)
246               (optimize (speed 3) (space 0)))               (optimize (speed 3) (space 0)
247                           (inhibit-warnings 3)))
248      (when (> x max-log)      (when (> x max-log)
249        (return-from dd-%exp        (return-from dd-%exp
250          (kernel:%make-double-double-float #.ext:double-float-positive-infinity          (kernel:%make-double-double-float #.ext:double-float-positive-infinity
# Line 354  Line 356 
356    ;;    ;;
357    (defun dd-%log (x)    (defun dd-%log (x)
358      (declare (type double-double-float x)      (declare (type double-double-float x)
359               (optimize (speed 3) (space 0)))               (optimize (speed 3) (space 0)
360                           (inhibit-warnings 3)))
361      ;; Separate mantissa from exponent      ;; Separate mantissa from exponent
362      (multiple-value-bind (x e)      (multiple-value-bind (x e)
363          (decode-float x)          (decode-float x)
# Line 483  Line 486 
486            (z 0w0)            (z 0w0)
487            (y 0w0))            (y 0w0))
488        (declare (type double-double-float x y z)        (declare (type double-double-float x y z)
489                 (optimize (speed 3)))                 (optimize (speed 3)
490                             (inhibit-warnings 3)))
491        (multiple-value-bind (x e)        (multiple-value-bind (x e)
492            (decode-float x)            (decode-float x)
493          (declare (type double-double-float x)          (declare (type double-double-float x)
# Line 545  Line 549 
549                                           3.889701475261939961851358705515223019890w14))))                                           3.889701475261939961851358705515223019890w14))))
550    (defun dd-%sinh (x)    (defun dd-%sinh (x)
551      (declare (type double-double-float x)      (declare (type double-double-float x)
552               (optimize (speed 3) (space 0)))               (optimize (speed 3) (space 0)
553                           (inhibit-warnings 3)))
554      (let ((a (abs x)))      (let ((a (abs x)))
555        (declare (type double-double-float a))        (declare (type double-double-float a))
556        (cond ((> a 1)        (cond ((> a 1)
# Line 580  Line 585 
585    
586  (defun dd-%cosh (x)  (defun dd-%cosh (x)
587    (declare (type double-double-float x)    (declare (type double-double-float x)
588             (optimize (speed 3) (space 0)))             (optimize (speed 3) (space 0)
589                         (inhibit-warnings 3)))
590    (let ((y (dd-%exp x)))    (let ((y (dd-%exp x)))
591      (scale-float (+ y (/ y)) -1)))      (scale-float (+ y (/ y)) -1)))
592    
# Line 602  Line 608 
608                                           1.365413143842835040443257277862054198329w8))))                                           1.365413143842835040443257277862054198329w8))))
609    (defun dd-%tanh (x)    (defun dd-%tanh (x)
610      (declare (type double-double-float x)      (declare (type double-double-float x)
611               (optimize (speed 3) (space 0)))               (optimize (speed 3) (space 0)
612                           (inhibit-warnings 3)))
613      (let ((z (abs x)))      (let ((z (abs x)))
614        (declare (type double-double-float z))        (declare (type double-double-float z))
615        (cond ((> z (* 0.5w0 max-log))        (cond ((> z (* 0.5w0 max-log))
# Line 652  Line 659 
659                                           ))))                                           ))))
660    (defun dd-%atanh (x)    (defun dd-%atanh (x)
661      (declare (type double-double-float x)      (declare (type double-double-float x)
662               (optimize (speed 3) (space 0)))               (optimize (speed 3) (space 0)
663                           (inhibit-warnings 3)))
664      (cond ((minusp x)      (cond ((minusp x)
665             (- (the double-double-float (dd-%atanh (- x)))))             (- (the double-double-float (dd-%atanh (- x)))))
666            ((< x 1w-12)            ((< x 1w-12)
# Line 692  Line 700 
700                                           6.226145049170155830806967678679167550122w4))))                                           6.226145049170155830806967678679167550122w4))))
701    (defun dd-%asinh (x)    (defun dd-%asinh (x)
702      (declare (type double-double-float x)      (declare (type double-double-float x)
703               (optimize (speed 3) (space 0)))               (optimize (speed 3) (space 0)
704                           (inhibit-warnings 3)))
705      (cond ((minusp x)      (cond ((minusp x)
706             (- (the double-double-float (dd-%asinh (- x)))))             (- (the double-double-float (dd-%asinh (- x)))))
707            #+nil            #+nil
# Line 734  Line 743 
743                                           ))))                                           ))))
744    (defun dd-%acosh (x)    (defun dd-%acosh (x)
745      (declare (type (double-double-float 1w0) x)      (declare (type (double-double-float 1w0) x)
746               (optimize (speed 3) (space 0)))               (optimize (speed 3) (space 0)
747                           (inhibit-warnings 3)))
748      (cond ((> x 1w17)      (cond ((> x 1w17)
749             (+ loge2 (dd-%log x)))             (+ loge2 (dd-%log x)))
750            (t            (t
# Line 776  Line 786 
786                                           ))))                                           ))))
787    (defun dd-%asin (x)    (defun dd-%asin (x)
788      (declare (type (double-double-float -1w0 1w0) x)      (declare (type (double-double-float -1w0 1w0) x)
789               (optimize (speed 3) (space 0)))               (optimize (speed 3) (space 0)
790                           (inhibit-warnings 3)))
791      (cond ((minusp x)      (cond ((minusp x)
792             (- (the double-double-float (dd-%asin (- x)))))             (- (the double-double-float (dd-%asin (- x)))))
793            #+nil            #+nil
# Line 802  Line 813 
813    
814  (defun dd-%acos (x)  (defun dd-%acos (x)
815    (declare (type (double-double-float -1w0 1w0) x)    (declare (type (double-double-float -1w0 1w0) x)
816             (optimize (speed 3) (space 0)))             (optimize (speed 3) (space 0)
817                         (inhibit-warnings 3)))
818    (cond ((< x -0.5w0)    (cond ((< x -0.5w0)
819           (- dd-pi           (- dd-pi
820              (* 2 (dd-%asin (sqrt (* 0.5w0 (+ 1 x)))))))              (* 2 (dd-%asin (sqrt (* 0.5w0 (+ 1 x)))))))
# Line 843  Line 855 
855    (declare (type double-double-float t3p8 tp8))    (declare (type double-double-float t3p8 tp8))
856    (defun dd-%atan (x)    (defun dd-%atan (x)
857      (declare (type double-double-float x)      (declare (type double-double-float x)
858               (optimize (speed 3) (space 0)))               (optimize (speed 3) (space 0)
859                           (inhibit-warnings 3)))
860      (when (minusp x)      (when (minusp x)
861        (return-from dd-%atan (- (the double-double-float (dd-%atan (- x))))))        (return-from dd-%atan (- (the double-double-float (dd-%atan (- x))))))
862      ;; range reduction      ;; range reduction
# Line 869  Line 882 
882    
883  (defun dd-%atan2 (y x)  (defun dd-%atan2 (y x)
884    (declare (type double-double-float x y)    (declare (type double-double-float x y)
885             (optimize (speed 3) (space 0)))             (optimize (speed 3) (space 0)
886                         (inhibit-warnings 3)))
887    (let ((code 0)    (let ((code 0)
888          (neg-x (minusp (float-sign x)))          (neg-x (minusp (float-sign x)))
889          (neg-y (minusp (float-sign y))))          (neg-y (minusp (float-sign y))))
# Line 1013  pi/4    11001001000011111101101010100010 Line 1027  pi/4    11001001000011111101101010100010
1027    
1028  (defun dd-%%sin (x)  (defun dd-%%sin (x)
1029    (declare (type double-double-float x)    (declare (type double-double-float x)
1030             (optimize (speed 2) (space 0)))             (optimize (speed 2) (space 0)
1031                         (inhibit-warnings 3)))
1032    (when (minusp x)    (when (minusp x)
1033      (return-from dd-%%sin (- (the double-double-float (dd-%%sin (- x))))))      (return-from dd-%%sin (- (the double-double-float (dd-%%sin (- x))))))
1034    ;; y = integer part of x/(pi/4).    ;; y = integer part of x/(pi/4).
# Line 1051  pi/4    11001001000011111101101010100010 Line 1066  pi/4    11001001000011111101101010100010
1066    
1067  (defun dd-%%cos (x)  (defun dd-%%cos (x)
1068    (declare (type double-double-float x)    (declare (type double-double-float x)
1069             (optimize (speed 2) (space 0)))             (optimize (speed 2) (space 0)
1070                         (inhibit-warnings 3)))
1071    (when (minusp x)    (when (minusp x)
1072      (return-from dd-%%cos (dd-%%cos (- x))))      (return-from dd-%%cos (dd-%%cos (- x))))
1073    ;; y = integer part of x/(pi/4).    ;; y = integer part of x/(pi/4).
# Line 1201  pi/4    11001001000011111101101010100010 Line 1217  pi/4    11001001000011111101101010100010
1217    ;; double-float and stored in PARTS.    ;; double-float and stored in PARTS.
1218    (defun dd-expand (x)    (defun dd-expand (x)
1219      (declare (double-double-float x)      (declare (double-double-float x)
1220               (optimize (speed 3) (space 0)))               (optimize (speed 3) (space 0)
1221                           (inhibit-warnings 3)))
1222      (multiple-value-bind (frac exp)      (multiple-value-bind (frac exp)
1223          (decode-float x)          (decode-float x)
1224        (declare (double-double-float frac)        (declare (double-double-float frac)
# Line 1214  pi/4    11001001000011111101101010100010 Line 1231  pi/4    11001001000011111101101010100010
1231        exp))        exp))
1232    (defun reduce-arg (x)    (defun reduce-arg (x)
1233      (declare (double-double-float x)      (declare (double-double-float x)
1234               (optimize (speed 3)))               (optimize (speed 3)
1235                           (inhibit-warnings 3)))
1236      (let* ((e0 (dd-expand x))      (let* ((e0 (dd-expand x))
1237             (n (sys:without-gcing             (n (sys:without-gcing
1238                 (%kernel-rem-pi/2 (vector-sap parts)                 (%kernel-rem-pi/2 (vector-sap parts)
# Line 1346  pi/4    11001001000011111101101010100010 Line 1364  pi/4    11001001000011111101101010100010
1364                         ))))                         ))))
1365    (defun dd-%log2 (x)    (defun dd-%log2 (x)
1366      (declare (type double-double-float x)      (declare (type double-double-float x)
1367               (optimize (speed 3) (space 0)))               (optimize (speed 3) (space 0)
1368                           (inhibit-warnings 3)))
1369      (multiple-value-bind (x e)      (multiple-value-bind (x e)
1370          (decode-float x)          (decode-float x)
1371        (declare (type double-double-float x)        (declare (type double-double-float x)
# Line 1414  pi/4    11001001000011111101101010100010 Line 1433  pi/4    11001001000011111101101010100010
1433                         ))))                         ))))
1434    (defun dd-%exp2 (x)    (defun dd-%exp2 (x)
1435      (declare (type double-double-float x)      (declare (type double-double-float x)
1436               (optimize (speed 3) (space 0)))               (optimize (speed 3) (space 0)
1437                           (inhibit-warnings 3)))
1438      (when (>= x 1024w0)      (when (>= x 1024w0)
1439        (return-from dd-%exp2        (return-from dd-%exp2
1440          (%make-double-double-float ext:double-float-positive-infinity          (%make-double-double-float ext:double-float-positive-infinity
# Line 1436  pi/4    11001001000011111101101010100010 Line 1456  pi/4    11001001000011111101101010100010
1456  (defun dd-%powil (x nn)  (defun dd-%powil (x nn)
1457    (declare (type double-double-float x)    (declare (type double-double-float x)
1458             (fixnum nn)             (fixnum nn)
1459             (optimize (speed 3) (space 0)))             (optimize (speed 3) (space 0)
1460                         (inhibit-warnings 3)))
1461    (when (zerop x)    (when (zerop x)
1462      (return-from dd-%powil      (return-from dd-%powil
1463        (cond ((zerop nn)        (cond ((zerop nn)
# Line 1524  pi/4    11001001000011111101101010100010 Line 1545  pi/4    11001001000011111101101010100010
1545    
1546  (defun dd-real-pow (x y)  (defun dd-real-pow (x y)
1547    (declare (type double-double-float x y)    (declare (type double-double-float x y)
1548             (optimize (speed 3) (space 0)))             (optimize (speed 3) (space 0)
1549                         (inhibit-warnings 3)))
1550    (let ((nflg 0)    (let ((nflg 0)
1551          (w (floor y)))          (w (floor y)))
1552      ;; nflg = 1 if x < 0 raised to integer power      ;; nflg = 1 if x < 0 raised to integer power
# Line 1587  pi/4    11001001000011111101101010100010 Line 1609  pi/4    11001001000011111101101010100010
1609      ;; invalid operation.      ;; invalid operation.
1610      (with-float-traps-masked (:underflow :overflow :invalid)      (with-float-traps-masked (:underflow :overflow :invalid)
1611        (let ((rho (+ (square x) (square y))))        (let ((rho (+ (square x) (square y))))
1612          (declare (optimize (speed 3) (space 0)))          (declare (optimize (speed 3) (space 0)
1613                               (inhibit-warnings 3)))
1614          (cond ((and (or (float-nan-p rho)          (cond ((and (or (float-nan-p rho)
1615                          (float-infinity-p rho))                          (float-infinity-p rho))
1616                      (or (float-infinity-p (abs x))                      (or (float-infinity-p (abs x))
# Line 1630  Z may be any number, but the result is a Line 1653  Z may be any number, but the result is a
1653    
1654        (locally        (locally
1655            ;; space 0 to get maybe-inline functions inlined.            ;; space 0 to get maybe-inline functions inlined.
1656            (declare (optimize (speed 3) (space 0)))            (declare (optimize (speed 3) (space 0)
1657                                 (inhibit-warnings 3)))
1658    
1659          (if (not (float-nan-p x))          (if (not (float-nan-p x))
1660              (setf rho (+ (scalb (abs x) (- k)) (sqrt rho))))              (setf rho (+ (scalb (abs x) (- k)) (sqrt rho))))
# Line 1673  This is for use with J /= 0 only when |z Line 1697  This is for use with J /= 0 only when |z
1697          (y (float (imagpart z) 1.0w0)))          (y (float (imagpart z) 1.0w0)))
1698      (multiple-value-bind (rho k)      (multiple-value-bind (rho k)
1699          (dd-cssqs z)          (dd-cssqs z)
1700        (declare (optimize (speed 3)))        (declare (optimize (speed 3)
1701                             (inhibit-warnings 3)))
1702        (let ((beta (max (abs x) (abs y)))        (let ((beta (max (abs x) (abs y)))
1703              (theta (min (abs x) (abs y))))              (theta (min (abs x) (abs y))))
1704          (complex (if (and (zerop k)          (complex (if (and (zerop k)
# Line 1720  Z may be any number, but the result is a Line 1745  Z may be any number, but the result is a
1745             ;; Shouldn't need this declare.             ;; Shouldn't need this declare.
1746             (declare (double-double-float x y))             (declare (double-double-float x y))
1747             (locally             (locally
1748                 (declare (optimize (speed 3)))                 (declare (optimize (speed 3)
1749                                      (inhibit-warnings 3)))
1750               (cond ((or (> x theta)               (cond ((or (> x theta)
1751                          (> (abs y) theta))                          (> (abs y) theta))
1752                      ;; To avoid overflow...                      ;; To avoid overflow...
# Line 1767  Z may be any number, but the result is a Line 1793  Z may be any number, but the result is a
1793          (y (float (imagpart z) 1.0w0)))          (y (float (imagpart z) 1.0w0)))
1794      (locally      (locally
1795          ;; space 0 to get maybe-inline functions inlined          ;; space 0 to get maybe-inline functions inlined
1796          (declare (optimize (speed 3) (space 0)))          (declare (optimize (speed 3) (space 0)
1797                               (inhibit-warnings 3)))
1798        (cond ((> (abs x)        (cond ((> (abs x)
1799                  #-(or linux hpux) #.(/ (%asinh most-positive-double-float) 4d0)                  #-(or linux hpux) #.(/ (%asinh most-positive-double-float) 4d0)
1800                  ;; This is more accurate under linux.                  ;; This is more accurate under linux.

Legend:
Removed from v.1.16  
changed lines
  Added in v.1.17

  ViewVC Help
Powered by ViewVC 1.1.5