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

revision 1.2 by rtoy, Fri Jun 30 18:41:22 2006 UTC revision 1.3 by rtoy, Fri Jul 7 18:26:43 2006 UTC
# Line 79  Line 79
79  ;; Evaluate polynomial  ;; Evaluate polynomial
80  (defun poly-eval (x coef)  (defun poly-eval (x coef)
81    (declare (type double-double-float x)    (declare (type double-double-float x)
82             (type (simple-array double-double-float (*)) coef))             (type (simple-array double-double-float (*)) coef)
83               (optimize (speed 3) (space 0)))
84    ;; y = C0 + C1*x + C2*x^2 + ...    ;; y = C0 + C1*x + C2*x^2 + ...
85    ;;    ;;
86    ;; But coefficients are stored in reverse (descending powers) order:    ;; But coefficients are stored in reverse (descending powers) order:
# Line 93  Line 94
94
95  (defun poly-eval-1 (x coef)  (defun poly-eval-1 (x coef)
96    (declare (type double-double-float x)    (declare (type double-double-float x)
97             (type (simple-array double-double-float (*)) coef))             (type (simple-array double-double-float (*)) coef)
98               (optimize (speed 3) (space 0)))
99    ;; Like poly-eval, except it assumes coef[N] = 1 and is omitted.    ;; Like poly-eval, except it assumes coef[N] = 1 and is omitted.
100    (let ((y 1w0))    (let ((y 1w0))
101      (declare (type double-double-float y))      (declare (type double-double-float y))
# Line 153  Line 155
155    (defun dd-%expm1 (x)    (defun dd-%expm1 (x)
156      "exp(x) - 1"      "exp(x) - 1"
157      (declare (type double-double-float x)      (declare (type double-double-float x)
158               (optimize (speed 3)))               (optimize (speed 3) (space 0)))
159      ;; Range reduction is accomplished by separating the argument      ;; Range reduction is accomplished by separating the argument
160      ;; into an integer k and fraction f such that      ;; into an integer k and fraction f such that
161      ;;      ;;
# Line 282  Line 284
284    ;; relative error = 2.2e-37, relative peak error spread = 9.2e-38.    ;; relative error = 2.2e-37, relative peak error spread = 9.2e-38.
285    (defun dd-%exp (x)    (defun dd-%exp (x)
286      (declare (type double-double-float x)      (declare (type double-double-float x)
287               (optimize (speed 3)))               (optimize (speed 3) (space 0)))
288      (when (> x max-log)      (when (> x max-log)
289        (return-from dd-%exp        (return-from dd-%exp
290          (kernel:%make-double-double-float #.ext:double-float-positive-infinity          (kernel:%make-double-double-float #.ext:double-float-positive-infinity
# Line 394  Line 396
396    ;;    ;;
397    (defun dd-%log (x)    (defun dd-%log (x)
398      (declare (type double-double-float x)      (declare (type double-double-float x)
399               (optimize (speed 3)))               (optimize (speed 3) (space 0)))
400      ;; Separate mantissa from exponent      ;; Separate mantissa from exponent
401      (multiple-value-bind (x e)      (multiple-value-bind (x e)
402          (decode-float x)          (decode-float x)
# Line 517  Line 519
519    ;; Theoretical peak relative error = 5.3e-37,    ;; Theoretical peak relative error = 5.3e-37,
520    ;; relative peak error spread = 2.3e-14    ;; relative peak error spread = 2.3e-14
521    (defun dd-%log1p (xm1)    (defun dd-%log1p (xm1)
522      (declare (type double-double-float xm1))      (declare (type double-double-float xm1)
523                 (optimize (space 0)))
524      (let ((x (+ xm1 1))      (let ((x (+ xm1 1))
525            (z 0w0)            (z 0w0)
526            (y 0w0))            (y 0w0))
# Line 584  Line 587
587                                           3.889701475261939961851358705515223019890w14))))                                           3.889701475261939961851358705515223019890w14))))
588    (defun dd-%sinh (x)    (defun dd-%sinh (x)
589      (declare (type double-double-float x)      (declare (type double-double-float x)
590               (optimize (speed 3)))               (optimize (speed 3) (space 0)))
591      (let ((a (abs x)))      (let ((a (abs x)))
592        (declare (type double-double-float a))        (declare (type double-double-float a))
593        (cond ((> a 1)        (cond ((> a 1)
# Line 618  Line 621
621               (* 0.5w0 (+ d (/ d (+ 1 d)))))))))               (* 0.5w0 (+ d (/ d (+ 1 d)))))))))
622
623  (defun dd-%cosh (x)  (defun dd-%cosh (x)
624    (declare (type double-double-float x))    (declare (type double-double-float x)
625               (optimize (speed 3) (space 0)))
626    (let ((y (dd-%exp x)))    (let ((y (dd-%exp x)))
627      (scale-float (+ y (/ y)) -1)))      (scale-float (+ y (/ y)) -1)))
628
# Line 640  Line 644
644                                           1.365413143842835040443257277862054198329w8))))                                           1.365413143842835040443257277862054198329w8))))
645    (defun dd-%tanh (x)    (defun dd-%tanh (x)
646      (declare (type double-double-float x)      (declare (type double-double-float x)
647               (optimize (speed 3)))               (optimize (speed 3) (space 0)))
648      (let ((z (abs x)))      (let ((z (abs x)))
649        (declare (type double-double-float z))        (declare (type double-double-float z))
650        (cond ((> z (* 0.5w0 max-log))        (cond ((> z (* 0.5w0 max-log))
# Line 690  Line 694
694                                           ))))                                           ))))
695    (defun dd-%atanh (x)    (defun dd-%atanh (x)
696      (declare (type double-double-float x)      (declare (type double-double-float x)
697               (optimize (speed 3)))               (optimize (speed 3) (space 0)))
698      (cond ((minusp x)      (cond ((minusp x)
699             (- (the double-double-float (dd-%atanh (- x)))))             (- (the double-double-float (dd-%atanh (- x)))))
700            ((< x 1w-12)            ((< x 1w-12)
# Line 730  Line 734
734                                           6.226145049170155830806967678679167550122w4))))                                           6.226145049170155830806967678679167550122w4))))
735    (defun dd-%asinh (x)    (defun dd-%asinh (x)
736      (declare (type double-double-float x)      (declare (type double-double-float x)
737               (optimize (speed 3)))               (optimize (speed 3) (space 0)))
738      (cond ((minusp x)      (cond ((minusp x)
739             (- (the double-double-float (dd-%asinh (- x)))))             (- (the double-double-float (dd-%asinh (- x)))))
740            #+nil            #+nil
# Line 772  Line 776
776                                           ))))                                           ))))
777    (defun dd-%acosh (x)    (defun dd-%acosh (x)
778      (declare (type double-double-float x)      (declare (type double-double-float x)
779               (optimize (speed 3)))               (optimize (speed 3) (space 0)))
780      (cond ((> x 1w17)      (cond ((> x 1w17)
781             (+ loge2 (dd-%log x)))             (+ loge2 (dd-%log x)))
782            (t            (t
# Line 813  Line 817
817                                           ))))                                           ))))
818    (defun dd-%asin (x)    (defun dd-%asin (x)
819      (declare (type double-double-float x)      (declare (type double-double-float x)
820               (optimize (speed 3)))               (optimize (speed 3) (space 0)))
821      (cond ((minusp x)      (cond ((minusp x)
822             (- (the double-double-float (dd-%asin (- x)))))             (- (the double-double-float (dd-%asin (- x)))))
823            #+nil            #+nil
# Line 843  Line 847
847
848  (defun dd-%acos (x)  (defun dd-%acos (x)
849    (declare (type double-double-float x)    (declare (type double-double-float x)
850             (optimize (speed 3)))             (optimize (speed 3) (space 0)))
851    (cond ((< x -0.5w0)    (cond ((< x -0.5w0)
852           (- dd-pi           (- dd-pi
853              (* 2 (dd-%asin (sqrt (* 0.5w0 (+ 1 x)))))))              (* 2 (dd-%asin (sqrt (* 0.5w0 (+ 1 x)))))))
# Line 884  Line 888
888    (declare (type double-double-float t3p8 tp8))    (declare (type double-double-float t3p8 tp8))
889    (defun dd-%atan (x)    (defun dd-%atan (x)
890      (declare (type double-double-float x)      (declare (type double-double-float x)
891               (optimize (speed 3)))               (optimize (speed 3) (space 0)))
892      (when (minusp x)      (when (minusp x)
893        (return-from dd-%atan (- (the double-double-float (dd-%atan (- x))))))        (return-from dd-%atan (- (the double-double-float (dd-%atan (- x))))))
894      ;; range reduction      ;; range reduction
# Line 910  Line 914
914
915  (defun dd-%atan2 (y x)  (defun dd-%atan2 (y x)
916    (declare (type double-double-float x y)    (declare (type double-double-float x y)
917             (optimize (speed 3)))             (optimize (speed 3) (space 0)))
918    (let ((code 0)    (let ((code 0)
919          (w 0w0))          (w 0w0))
920      (declare (type (integer 0 3) code)      (declare (type (integer 0 3) code)
# Line 1024  pi/4    11001001000011111101101010100010 Line 1028  pi/4    11001001000011111101101010100010
1028    (scale-float (float #b0111101111100101010001100110110011110011010011101001000011000110110011000000101011000010100110110111110010 1w0) (* 4 -106)))    (scale-float (float #b0111101111100101010001100110110011110011010011101001000011000110110011000000101011000010100110110111110010 1w0) (* 4 -106)))
1029
1030  (defun dd-%sin (x)  (defun dd-%sin (x)
1031    (declare (type double-double-float x))    (declare (type double-double-float x)
1032               (optimize (speed 3) (space 0)))
1033    (when (minusp x)    (when (minusp x)
1034      (return-from dd-%sin (- (the double-double-float (dd-%sin (- x))))))      (return-from dd-%sin (- (the double-double-float (dd-%sin (- x))))))
1035    ;; y = integer part of x/(pi/4).    ;; y = integer part of x/(pi/4).
# Line 1061  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 3)))             (optimize (speed 3) (space 0)))
1070    (when (minusp x)    (when (minusp x)
1071      (return-from dd-%cos (dd-%cos (- x))))      (return-from dd-%cos (dd-%cos (- x))))
1072    ;; y = integer part of x/(pi/4).    ;; y = integer part of x/(pi/4).
# Line 1125  pi/4    11001001000011111101101010100010 Line 1130  pi/4    11001001000011111101101010100010
1130                         ))))                         ))))
1131    (defun dd-tancot (xx cotflag)    (defun dd-tancot (xx cotflag)
1132      (declare (type double-double-float xx)      (declare (type double-double-float xx)
1133               (optimize (speed 3)))               (optimize (speed 3) (space 0)))
1134      (let ((x 0w0)      (let ((x 0w0)
1135            (sign 1))            (sign 1))
1136        (declare (type double-double-float x)        (declare (type double-double-float x)
# Line 1237  pi/4    11001001000011111101101010100010 Line 1242  pi/4    11001001000011111101101010100010
1242                         ))))                         ))))
1243    (defun dd-%log2 (x)    (defun dd-%log2 (x)
1244      (declare (type double-double-float x)      (declare (type double-double-float x)
1245               (optimize (speed 3)))               (optimize (speed 3) (space 0)))
1246      (multiple-value-bind (x e)      (multiple-value-bind (x e)
1247          (decode-float x)          (decode-float x)
1248        (declare (type double-double-float x)        (declare (type double-double-float x)
# Line 1305  pi/4    11001001000011111101101010100010 Line 1310  pi/4    11001001000011111101101010100010
1310                         ))))                         ))))
1311    (defun dd-%exp2 (x)    (defun dd-%exp2 (x)
1312      (declare (type double-double-float x)      (declare (type double-double-float x)
1313               (optimize (speed 3)))               (optimize (speed 3) (space 0)))
1314      (when (>= x 1024w0)      (when (>= x 1024w0)
1315        (return-from dd-%exp2        (return-from dd-%exp2
1316          (%make-double-double-float ext:double-float-positive-infinity          (%make-double-double-float ext:double-float-positive-infinity
# Line 1327  pi/4    11001001000011111101101010100010 Line 1332  pi/4    11001001000011111101101010100010
1332  (defun dd-%powil (x nn)  (defun dd-%powil (x nn)
1333    (declare (type double-double-float x)    (declare (type double-double-float x)
1334             (fixnum nn)             (fixnum nn)
1335             (optimize (speed 3)))             (optimize (speed 3) (space 0)))
1336    (when (zerop x)    (when (zerop x)
1337      (return-from dd-%powil      (return-from dd-%powil
1338        (cond ((zerop nn)        (cond ((zerop nn)
# Line 1415  pi/4    11001001000011111101101010100010 Line 1420  pi/4    11001001000011111101101010100010
1420
1421  (defun dd-real-pow (x y)  (defun dd-real-pow (x y)
1422    (declare (type double-double-float x y)    (declare (type double-double-float x y)
1423             (optimize (speed 3)))             (optimize (speed 3) (space 0)))
1424    (let ((nflg 0)    (let ((nflg 0)
1425          (w (floor y)))          (w (floor y)))
1426      ;; nflg = 1 if x < 0 raised to integer power      ;; nflg = 1 if x < 0 raised to integer power

Legend:
 Removed from v.1.2 changed lines Added in v.1.3