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

revision 1.15 by rtoy, Thu Jun 21 15:50:50 2007 UTC revision 1.16 by rtoy, Wed Jun 27 16:38:54 2007 UTC
# Line 70  Line 70
70  (defun poly-eval (x coef)  (defun poly-eval (x coef)
71    (declare (type double-double-float x)    (declare (type double-double-float x)
72             (type (simple-array double-double-float (*)) coef)             (type (simple-array double-double-float (*)) coef)
73             (optimize (speed 3) (space 0)))             (optimize (speed 3) (space 0) (inhibit-warnings 3)))
74    ;; y = C0 + C1*x + C2*x^2 + ...    ;; y = C0 + C1*x + C2*x^2 + ...
75    ;;    ;;
76    ;; But coefficients are stored in reverse (descending powers) order:    ;; But coefficients are stored in reverse (descending powers) order:
# Line 85  Line 85
85  (defun poly-eval-1 (x coef)  (defun poly-eval-1 (x coef)
86    (declare (type double-double-float x)    (declare (type double-double-float x)
87             (type (simple-array double-double-float (*)) coef)             (type (simple-array double-double-float (*)) coef)
88             (optimize (speed 3) (space 0)))             (optimize (speed 3) (space 0) (inhibit-warnings 3)))
89    ;; Like poly-eval, except it assumes coef[N] = 1 and is omitted.    ;; Like poly-eval, except it assumes coef[N] = 1 and is omitted.
90    (let ((y 1w0))    (let ((y 1w0))
91      (declare (type double-double-float y))      (declare (type double-double-float y))
# Line 733  Line 733
733                                           3.190482951215438078279772140481195226621w7                                           3.190482951215438078279772140481195226621w7
734                                           ))))                                           ))))
735    (defun dd-%acosh (x)    (defun dd-%acosh (x)
736      (declare (type double-double-float x)      (declare (type (double-double-float 1w0) x)
737               (optimize (speed 3) (space 0)))               (optimize (speed 3) (space 0)))
738      (cond ((> x 1w17)      (cond ((> x 1w17)
739             (+ loge2 (dd-%log x)))             (+ loge2 (dd-%log x)))
# Line 746  Line 746
746                     (t                     (t
747                      (dd-%log (+ x (sqrt (* z (+ 1 x))))))))))))                      (dd-%log (+ x (sqrt (* z (+ 1 x))))))))))))
748
749
750  (let ((P (make-array 10 :element-type 'double-double-float  (let ((P (make-array 10 :element-type 'double-double-float
751                       :initial-contents '(                       :initial-contents '(
752                                           -8.067112765482705313585175280952515549833w-1                                           -8.067112765482705313585175280952515549833w-1
# Line 774  Line 775
775                                           7.122795760168575261226395089432959614179w4                                           7.122795760168575261226395089432959614179w4
776                                           ))))                                           ))))
777    (defun dd-%asin (x)    (defun dd-%asin (x)
778      (declare (type double-double-float x)      (declare (type (double-double-float -1w0 1w0) x)
779               (optimize (speed 3) (space 0)))               (optimize (speed 3) (space 0)))
780      (cond ((minusp x)      (cond ((minusp x)
781             (- (the double-double-float (dd-%asin (- x)))))             (- (the double-double-float (dd-%asin (- x)))))
# Line 782  Line 783
783            ((< x 1w-8)            ((< x 1w-8)
784             x)             x)
785            (t            (t
786             (let ((flag nil)             (multiple-value-bind (flag z zz)
787                   (z 0w0)                 (cond ((> x 0.5w0)
788                   (zz 0w0))                        (let* ((zz1 (- 0.5w0 x))
789                                 (zz (scale-float (+ zz1 0.5w0) -1)))
790                            (values t (sqrt zz) zz)))
791                         (t
792                          (values nil x (* x x))))
793               (declare (type double-double-float z zz))               (declare (type double-double-float z zz))
794               (cond ((> x 0.5w0)
(setf zz (- 0.5w0 x))
(setf zz (scale-float (+ zz 0.5w0) -1))
(setf z (sqrt zz))
(setf flag t))
(t
(setf z x)
(setf zz (* z z))
(setf flag nil)))
795               (let ((p (* zz (/ (poly-eval zz p)               (let ((p (* zz (/ (poly-eval zz p)
796                                 (poly-eval-1 zz q)))))                                 (poly-eval-1 zz q)))))
797                 (setf z (+ (* z p) z))                 (setf z (+ (* z p) z))
# Line 804  Line 801
801                 z))))))                 z))))))
802
803  (defun dd-%acos (x)  (defun dd-%acos (x)
804    (declare (type double-double-float x)    (declare (type (double-double-float -1w0 1w0) x)
805             (optimize (speed 3) (space 0)))             (optimize (speed 3) (space 0)))
806    (cond ((< x -0.5w0)    (cond ((< x -0.5w0)
807           (- dd-pi           (- dd-pi
# Line 1016  pi/4    11001001000011111101101010100010 Line 1013  pi/4    11001001000011111101101010100010
1013
1014  (defun dd-%%sin (x)  (defun dd-%%sin (x)
1015    (declare (type double-double-float x)    (declare (type double-double-float x)
1016             (optimize (speed 3) (space 0)))             (optimize (speed 2) (space 0)))
1017    (when (minusp x)    (when (minusp x)
1018      (return-from dd-%%sin (- (the double-double-float (dd-%%sin (- x))))))      (return-from dd-%%sin (- (the double-double-float (dd-%%sin (- x))))))
1019    ;; y = integer part of x/(pi/4).    ;; y = integer part of x/(pi/4).
# Line 1028  pi/4    11001001000011111101101010100010 Line 1025  pi/4    11001001000011111101101010100010
1025
1026      (let ((j (truncate z))      (let ((j (truncate z))
1027            (sign 1))            (sign 1))
1028          (declare (type (integer -1 1) sign))
1029        (unless (zerop (logand j 1))        (unless (zerop (logand j 1))
1030          (incf j)          (incf j)
1031          (incf y))          (incf y))
# Line 1053  pi/4    11001001000011111101101010100010 Line 1051  pi/4    11001001000011111101101010100010
1051
1052  (defun dd-%%cos (x)  (defun dd-%%cos (x)
1053    (declare (type double-double-float x)    (declare (type double-double-float x)
1054             (optimize (speed 3) (space 0)))             (optimize (speed 2) (space 0)))
1055    (when (minusp x)    (when (minusp x)
1056      (return-from dd-%%cos (dd-%%cos (- x))))      (return-from dd-%%cos (dd-%%cos (- x))))
1057    ;; y = integer part of x/(pi/4).    ;; y = integer part of x/(pi/4).
# Line 1117  pi/4    11001001000011111101101010100010 Line 1115  pi/4    11001001000011111101101010100010
1115                         ))))                         ))))
1116    (defun dd-tancot (xx cotflag)    (defun dd-tancot (xx cotflag)
1117      (declare (type double-double-float xx)      (declare (type double-double-float xx)
1118               (optimize (speed 3) (space 0)))               (optimize (speed 2) (space 0)))
1119      (let ((x 0w0)      (let ((x 0w0)
1120            (sign 1))            (sign 1))
1121        (declare (type double-double-float x)        (declare (type double-double-float x)
# Line 1231  pi/4    11001001000011111101101010100010 Line 1229  pi/4    11001001000011111101101010100010
1229        (values n sum))))        (values n sum))))
1230
1231
1232    (declaim (ftype (function (double-double-float) double-double-float)
1233                    dd-%sin))
1234  (defun dd-%sin (x)  (defun dd-%sin (x)
1235    (declare (double-double-float x))    (declare (double-double-float x))
1236    (cond ((minusp (float-sign x))    (cond ((minusp (float-sign x))
# Line 1241  pi/4    11001001000011111101101010100010 Line 1241  pi/4    11001001000011111101101010100010
1241           ;; Argument reduction needed           ;; Argument reduction needed
1242           (multiple-value-bind (n reduced)           (multiple-value-bind (n reduced)
1243               (reduce-arg x)               (reduce-arg x)
1244             (case (logand n 3)             (ecase (logand n 3)
1245               (0 (dd-%%sin reduced))               (0 (dd-%%sin reduced))
1246               (1 (dd-%%cos reduced))               (1 (dd-%%cos reduced))
1247               (2 (- (dd-%%sin reduced)))               (2 (- (dd-%%sin reduced)))
1248               (3 (- (dd-%%cos reduced))))))))               (3 (- (dd-%%cos reduced))))))))
1249
1250    (declaim (ftype (function (double-double-float) double-double-float)
1251                    dd-%cos))
1252  (defun dd-%cos (x)  (defun dd-%cos (x)
1253    (declare (double-double-float x))    (declare (double-double-float x))
1254    (cond ((minusp x)    (cond ((minusp x)
# Line 1257  pi/4    11001001000011111101101010100010 Line 1259  pi/4    11001001000011111101101010100010
1259           ;; Argument reduction needed           ;; Argument reduction needed
1260           (multiple-value-bind (n reduced)           (multiple-value-bind (n reduced)
1261               (reduce-arg x)               (reduce-arg x)
1262             (case (logand n 3)             (ecase (logand n 3)
1263               (0 (dd-%%cos reduced))               (0 (dd-%%cos reduced))
1264               (1 (- (dd-%%sin reduced)))               (1 (- (dd-%%sin reduced)))
1265               (2 (- (dd-%%cos reduced)))               (2 (- (dd-%%cos reduced)))
1266               (3 (dd-%%sin reduced)))))))               (3 (dd-%%sin reduced)))))))
1267
1268    (declaim (ftype (function (double-double-float) double-double-float)
1269                    dd-%tan))
1270  (defun dd-%tan (x)  (defun dd-%tan (x)
1271    (declare (double-double-float x))    (declare (double-double-float x))
1272    (cond ((minusp (float-sign x))    (cond ((minusp (float-sign x))

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