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

revision 1.1.2.4 by rtoy, Thu Jun 29 14:50:27 2006 UTC revision 1.1.2.5 by rtoy, Thu Jun 29 20:49:12 2006 UTC
# Line 517  Line 517
517    ;; Theoretical peak relative error = 5.3e-37,    ;; Theoretical peak relative error = 5.3e-37,
518    ;; relative peak error spread = 2.3e-14    ;; relative peak error spread = 2.3e-14
519    (defun dd-%log1p (xm1)    (defun dd-%log1p (xm1)
520      (declare (type double-double-float x))      (declare (type double-double-float xm1))
521      (let ((x (+ xm1 1))      (let ((x (+ xm1 1))
522            (z 0w0)            (z 0w0)
523            (y 0w0))            (y 0w0))
524        (declare (type double-double-float x y z))        (declare (type double-double-float x y z)
525                   (optimize (speed 3)))
526        (multiple-value-bind (x e)        (multiple-value-bind (x e)
527            (decode-float x)            (decode-float x)
528            (declare (type double-double-float x)
529                     (type double-float-exponent e))
530          (cond ((or (> e 2)          (cond ((or (> e 2)
531                     (< e -2))                     (< e -2))
532                 ;; Log using log(x) = z + z^3*P(z^2)/Q(z^2)                 ;; Log using log(x) = z + z^3*P(z^2)/Q(z^2)
# Line 1023  pi/4    11001001000011111101101010100010 Line 1026  pi/4    11001001000011111101101010100010
1026  (defun dd-%sin (x)  (defun dd-%sin (x)
1027    (declare (type double-double-float x))    (declare (type double-double-float x))
1028    (when (minusp x)    (when (minusp x)
1029      (return-from dd-%sin (- (dd-%sin (- x)))))      (return-from dd-%sin (- (the double-double-float (dd-%sin (- x))))))
1030    ;; y = integer part of x/(pi/4).    ;; y = integer part of x/(pi/4).
1031    (let* ((y (float (floor (/ x dd-pi/4)) 1w0))    (let* ((y (float (floor (/ x dd-pi/4)) 1w0))
1032           (z (scale-float y -4)))           (z (scale-float y -4)))
# Line 1071  pi/4    11001001000011111101101010100010 Line 1074  pi/4    11001001000011111101101010100010
1074      (let ((i (truncate z))      (let ((i (truncate z))
1075            (j 0)            (j 0)
1076            (sign 1))            (sign 1))
1077          (declare (type (integer 0 7) j)
1078                   (type (integer -1 1) sign))
1079        (unless (zerop (logand i 1))        (unless (zerop (logand i 1))
1080          (incf i)          (incf i)
1081          (incf y))          (incf y))
# Line 1082  pi/4    11001001000011111101101010100010 Line 1087  pi/4    11001001000011111101101010100010
1087        (when (> j 1)        (when (> j 1)
1088          (setf sign (- sign)))          (setf sign (- sign)))
1089
1090        ;; Extended precision modular arithmetic        ;; Extended precision modular arithmetic.  This is basically
1091          ;; computing x - y*(pi/4) accurately so that |z| < pi/4.
1092        (setf z (- (- (- x (* y dp1))        (setf z (- (- (- x (* y dp1))
1093                      (* y dp2))                      (* y dp2))
1094                   (* y dp3)))                   (* y dp3)))
# Line 1091  pi/4    11001001000011111101101010100010 Line 1097  pi/4    11001001000011111101101010100010
1097                  (= j 2))                  (= j 2))
1098              (setf y (+ z (* z (* zz (poly-eval zz sincof)))))              (setf y (+ z (* z (* zz (poly-eval zz sincof)))))
1099              (setf y (+ (- 1 (scale-float zz -1))              (setf y (+ (- 1 (scale-float zz -1))
1100                         (* zz zz (poly-eval zz coscof)))))                         (* zz (poly-eval zz coscof) zz))))
1101          (if (< sign 0)          (if (< sign 0)
1102              (- y)              (- y)
1103              y)))))              y)))))
# Line 1387  pi/4    11001001000011111101101010100010 Line 1393  pi/4    11001001000011111101101010100010
1393                  (t                  (t
1394                   (setf y x)))                   (setf y x)))
1395            (let ((w x))            (let ((w x))
1396                (declare (type double-double-float w))
1397              (setf n (ash n -1))              (setf n (ash n -1))
1398              (loop while (not (zerop n))              (loop while (not (zerop n))
1399                 do                 do
# Line 1426  pi/4    11001001000011111101101010100010 Line 1433  pi/4    11001001000011111101101010100010
1433              (t              (t
1434               (when (/= w y)               (when (/= w y)
1435                 ;; noninteger power of negative number                 ;; noninteger power of negative number
1436                 (let ((p (dd-real-pow (abs x) y))                 (let ((p (the double-double-float (dd-real-pow (abs x) y)))
1437                       (y*pi (* y dd-pi)))                       (y*pi (* y dd-pi)))
1438                   (return-from dd-real-pow (complex (* p (dd-%cos y*pi))                   (return-from dd-real-pow (complex (* p (dd-%cos y*pi))
1439                                                     (* p (dd-%sin y*pi))))))                                                     (* p (dd-%sin y*pi))))))
# Line 1451  pi/4    11001001000011111101101010100010 Line 1458  pi/4    11001001000011111101101010100010
1458    (dd-real-pow x y))    (dd-real-pow x y))
1459
1460
1461    ;; These are essentially the same as in irrat.lisp, but very slightly
1462    ;; modified to operate on double-double-floats.
1463  (defun dd-cssqs (z)  (defun dd-cssqs (z)
1464    ;; 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
1465    ;; result is r + i*k, where k is an integer.    ;; result is r + i*k, where k is an integer.

Legend:
 Removed from v.1.1.2.4 changed lines Added in v.1.1.2.5