September 2007 - oct-cvs - mailman3.common-lisp.net

[oct-cvs] Oct commit: oct qd-io.lisp
by rtoy 24 Sep '07

24 Sep '07

Update of /project/oct/cvsroot/oct In directory clnet:/tmp/cvs-serv16369 Modified Files: qd-io.lisp Log Message: A new version of MAKE-FLOAT that converts the number to rational before converting to a quad-double. This reduces round-off errors. This still needs work, I think. --- /project/oct/cvsroot/oct/qd-io.lisp 2007/09/24 02:37:31 1.16 +++ /project/oct/cvsroot/oct/qd-io.lisp 2007/09/24 21:32:15 1.17 @@ -305,11 +305,13 @@ (t (qd-output-aux arg stream)))) +;; This version has problems with roundoff. +#+nil (defun make-float (sign int-part frac-part scale exp) (declare (type (member -1 1) sign) (type unsigned-byte int-part frac-part) (fixnum scale exp)) - #+(or) + ;;#+(or) (progn (format t "sign = ~A~%" sign) (format t "int-part = ~A~%" int-part) @@ -319,7 +321,7 @@ (let ((int (cl:+ (cl:* int-part (expt 10 scale)) frac-part)) (power (cl:- exp scale))) - #+(or) + ;;#+(or) (format t "~A * ~A * 10^(~A)~%" sign int power) (let* ((len (integer-length int))) #+(or) @@ -379,6 +381,80 @@ (neg-qd (mul-qd xx yy)) (mul-qd xx yy)))))))) +(defun make-float (sign int-part frac-part scale exp) + (declare (type (member -1 1) sign) + (type unsigned-byte int-part frac-part) + (fixnum scale exp)) + (flet ((convert-int (int) + ;; Convert integer INT to a quad-double. + (let ((len (integer-length int))) + (cond ((<= len 53) + ;; The simple case that fits in a double-float + (let ((xx (make-qd-d (float int 1d0)))) + xx)) + (t + ;; The complicated case. We look at the top 5*53 + ;; bits and create doubles from them (properly + ;; scaled) and then combine into a quad-double. + ;; Looking at only 4*53 (212 bits) isn't accurate + ;; enough. In particulare 10^100 isn't converted + ;; as accurately as desired if only 212 bits are + ;; used. + ;; + ;; If the integer doesn't have 5*53 bits, left + ;; shift it until it does, and set the length to + ;; 5*53, so that the scaling is done properly. + (let* + ((new-int (if (< len (* 5 53)) + (progn + (setf len (* 5 53)) + (ash int (- (* 5 53) len))) + int)) + (q0 (ldb (byte 53 (cl:- len 53)) new-int)) + (q1 (ldb (byte 53 (cl:- len (* 2 53))) new-int)) + (q2 (ldb (byte 53 (cl:- len (* 3 53))) new-int)) + (q3 (ldb (byte 53 (cl:- len (* 4 53))) new-int)) + (q4 (ldb (byte 53 (cl:- len (* 5 53))) new-int)) + (xx (multiple-value-call #'%make-qd-d + (renorm-5 (scale-float (float q0 1d0) + (cl:- len 53)) + (scale-float (float q1 1d0) + (cl:- len (* 2 53))) + (scale-float (float q2 1d0) + (cl:- len (* 3 53))) + (scale-float (float q3 1d0) + (cl:- len (* 4 53))) + (scale-float (float q4 1d0) + (cl:- len (* 5 53))))))) + #+(or) + (progn + (format t "xx = ~A~%" xx) + #+(or) + (format t " = ~/qd::qd-format/~%" xx) + (format t "yy = ~A~%" yy) + #+(or) + (format t " = ~/qd::qd-format/~%" yy) + (format t "q0 = ~X (~A)~%" q0 + (scale-float (float q0 1d0) + (cl:- len 53))) + (format t "q1 = ~X (~A)~%" q1 + (scale-float (float q1 1d0) + (cl:- len (* 2 53)))) + #+(or) + (format t "~/qdi::qd-format/~%" (mul-qd xx yy))) + xx)))))) + (let* ((rat (* (cl:+ (cl:* int-part (expt 10 scale)) + frac-part) + (expt 10 (cl:- exp scale)))) + (top (numerator rat)) + (bot (denominator rat))) + (div-qd (let ((top-qd (convert-int top))) + (if (minusp sign) + (neg-qd top-qd) + top-qd)) + (convert-int bot))))) + + ;; This seems to work, but really needs to be rewritten! (defun read-qd (stream) (labels ((read-digits (s)

1 0

[oct-cvs] Oct commit: oct qd-methods.lisp
by rtoy 24 Sep '07

24 Sep '07

Update of /project/oct/cvsroot/oct In directory clnet:/tmp/cvs-serv16245 Modified Files: qd-methods.lisp Log Message: Minor optimization converting a rational to a qd-real. If the denominator is 1, skip the qd-real division. --- /project/oct/cvsroot/oct/qd-methods.lisp 2007/09/20 21:04:05 1.60 +++ /project/oct/cvsroot/oct/qd-methods.lisp 2007/09/24 21:30:07 1.61 @@ -81,9 +81,11 @@ ;; We should do something better than this. (let ((top (numerator x)) (bot (denominator x))) - (make-instance 'qd-real - :value (div-qd (qdi::make-float (signum top) (abs top) 0 0 0) - (qdi::make-float (signum bot) (abs bot) 0 0 0))))) + (if (= bot 1) + (make-instance 'qd-real :value (qdi::make-float (signum top) (abs top) 0 0 0)) + (make-instance 'qd-real + :value (div-qd (qdi::make-float (signum top) (abs top) 0 0 0) + (qdi::make-float (signum bot) (abs bot) 0 0 0)))))) (defmethod add1 ((a number))

1 0

[oct-cvs] Oct commit: oct qd-io.lisp
by rtoy 23 Sep '07

23 Sep '07

Update of /project/oct/cvsroot/oct In directory clnet:/tmp/cvs-serv23811 Modified Files: qd-io.lisp Log Message: Use more bits (265 instead of 212) when converting a rational to a quad-double. This fixes the issue that converting 10^100 to a quad-double isn't as accurate as it could be. --- /project/oct/cvsroot/oct/qd-io.lisp 2007/09/12 02:31:14 1.15 +++ /project/oct/cvsroot/oct/qd-io.lisp 2007/09/24 02:37:31 1.16 @@ -340,29 +340,25 @@ (neg-qd (mul-qd xx yy)) (mul-qd xx yy)))) (t - (let* #+nil - ((hi (ldb (byte 106 (cl:- len 106)) int)) - (lo (ldb (byte 106 (cl:- len 212)) int)) - (xx (make-qd-dd (cl:* sign (scale-float (float hi 1w0) - (cl:- len 106))) - (cl:* sign (scale-float (float lo 1w0) - (cl:- len 106 106))))) - (yy (npow (make-qd-d 10d0) - power))) - ((q0 (ldb (byte 53 (cl:- len 53)) int)) - (q1 (ldb (byte 53 (cl:- len (* 2 53))) int)) - (q2 (ldb (byte 53 (cl:- len (* 3 53))) int)) - (q3 (ldb (byte 53 (cl:- len (* 4 53))) int)) - (xx (make-qd-d (scale-float (float q0 1d0) - (cl:- len 53)) - (scale-float (float q1 1d0) - (cl:- len (* 2 53))) - (scale-float (float q2 1d0) - (cl:- len (* 3 53))) - (scale-float (float q3 1d0) - (cl:- len (* 4 53))))) - (yy (npow (make-qd-d 10d0) - power))) + (let* + ((q0 (ldb (byte 53 (cl:- len 53)) int)) + (q1 (ldb (byte 53 (cl:- len (* 2 53))) int)) + (q2 (ldb (byte 53 (cl:- len (* 3 53))) int)) + (q3 (ldb (byte 53 (cl:- len (* 4 53))) int)) + (q4 (ldb (byte 53 (cl:- len (* 5 53))) int)) + (xx (multiple-value-call #'%make-qd-d + (renorm-5 (scale-float (float q0 1d0) + (cl:- len 53)) + (scale-float (float q1 1d0) + (cl:- len (* 2 53))) + (scale-float (float q2 1d0) + (cl:- len (* 3 53))) + (scale-float (float q3 1d0) + (cl:- len (* 4 53))) + (scale-float (float q4 1d0) + (cl:- len (* 5 53)))))) + (yy (npow (make-qd-d 10d0) + power))) #+(or) (progn (format t "xx = ~A~%" xx)

1 0

[oct-cvs] Oct commit: oct qd-methods.lisp qd-package.lisp
by rtoy 20 Sep '07

20 Sep '07

Update of /project/oct/cvsroot/oct In directory clnet:/tmp/cvs-serv8223 Modified Files: qd-methods.lisp qd-package.lisp Log Message: Add RATIONALIZE methods. Algorithm graciously provided by Bruno Haible. qd-package.lisp: o Shadow RATIONALIZE qd-methods.lisp: o Add RATIONALIZE methods for CL:REAL's and QD-REAL's.. --- /project/oct/cvsroot/oct/qd-methods.lisp 2007/09/19 17:30:04 1.59 +++ /project/oct/cvsroot/oct/qd-methods.lisp 2007/09/20 21:04:05 1.60 @@ -835,6 +835,102 @@ (cl:rational x2) (cl:rational x3)))) +(defmethod rationalize ((x cl:real)) + (cl:rationalize x)) + +;;; The algorithm here is the method described in CLISP. Bruno Haible has +;;; graciously given permission to use this algorithm. He says, "You can use +;;; it, if you present the following explanation of the algorithm." +;;; +;;; Algorithm (recursively presented): +;;; If x is a rational number, return x. +;;; If x = 0.0, return 0. +;;; If x < 0.0, return (- (rationalize (- x))). +;;; If x > 0.0: +;;; Call (integer-decode-float x). It returns a m,e,s=1 (mantissa, +;;; exponent, sign). +;;; If m = 0 or e >= 0: return x = m*2^e. +;;; Search a rational number between a = (m-1/2)*2^e and b = (m+1/2)*2^e +;;; with smallest possible numerator and denominator. +;;; Note 1: If m is a power of 2, we ought to take a = (m-1/4)*2^e. +;;; But in this case the result will be x itself anyway, regardless of +;;; the choice of a. Therefore we can simply ignore this case. +;;; Note 2: At first, we need to consider the closed interval [a,b]. +;;; but since a and b have the denominator 2^(|e|+1) whereas x itself +;;; has a denominator <= 2^|e|, we can restrict the seach to the open +;;; interval (a,b). +;;; So, for given a and b (0 < a < b) we are searching a rational number +;;; y with a <= y <= b. +;;; Recursive algorithm fraction_between(a,b): +;;; c := (ceiling a) +;;; if c < b +;;; then return c ; because a <= c < b, c integer +;;; else +;;; ; a is not integer (otherwise we would have had c = a < b) +;;; k := c-1 ; k = floor(a), k < a < b <= k+1 +;;; return y = k + 1/fraction_between(1/(b-k), 1/(a-k)) +;;; ; note 1 <= 1/(b-k) < 1/(a-k) +;;; +;;; You can see that we are actually computing a continued fraction expansion. +;;; +;;; Algorithm (iterative): +;;; If x is rational, return x. +;;; Call (integer-decode-float x). It returns a m,e,s (mantissa, +;;; exponent, sign). +;;; If m = 0 or e >= 0, return m*2^e*s. (This includes the case x = 0.0.) +;;; Create rational numbers a := (2*m-1)*2^(e-1) and b := (2*m+1)*2^(e-1) +;;; (positive and already in lowest terms because the denominator is a +;;; power of two and the numerator is odd). +;;; Start a continued fraction expansion +;;; p[-1] := 0, p[0] := 1, q[-1] := 1, q[0] := 0, i := 0. +;;; Loop +;;; c := (ceiling a) +;;; if c >= b +;;; then k := c-1, partial_quotient(k), (a,b) := (1/(b-k),1/(a-k)), +;;; goto Loop +;;; finally partial_quotient(c). +;;; Here partial_quotient(c) denotes the iteration +;;; i := i+1, p[i] := c*p[i-1]+p[i-2], q[i] := c*q[i-1]+q[i-2]. +;;; At the end, return s * (p[i]/q[i]). +;;; This rational number is already in lowest terms because +;;; p[i]*q[i-1]-p[i-1]*q[i] = (-1)^i. +;;; +(defmethod rationalize ((x qd-real)) + ;; This is a fairly straigtforward implementation of the iterative + ;; algorithm above. + (multiple-value-bind (frac expo sign) + (integer-decode-float x) + (cond ((or (zerop frac) (>= expo 0)) + (if (minusp sign) + (- (ash frac expo)) + (ash frac expo))) + (t + ;; expo < 0 and (2*m-1) and (2*m+1) are coprime to 2^(1-e), + ;; so build the fraction up immediately, without having to do + ;; a gcd. + (let ((a (/ (- (* 2 frac) 1) (ash 1 (- 1 expo)))) + (b (/ (+ (* 2 frac) 1) (ash 1 (- 1 expo)))) + (p0 0) + (q0 1) + (p1 1) + (q1 0)) + (do ((c (ceiling a) (ceiling a))) + ((< c b) + (let ((top (+ (* c p1) p0)) + (bot (+ (* c q1) q0))) + (/ (if (minusp sign) + (- top) + top) + bot))) + (let* ((k (- c 1)) + (p2 (+ (* k p1) p0)) + (q2 (+ (* k q1) q0))) + (psetf a (/ (- b k)) + b (/ (- a k))) + (setf p0 p1 + q0 q1 + p1 p2 + q1 q2)))))))) (define-compiler-macro + (&whole form &rest args) (declare (ignore form)) --- /project/oct/cvsroot/oct/qd-package.lisp 2007/09/18 12:46:36 1.40 +++ /project/oct/cvsroot/oct/qd-package.lisp 2007/09/20 21:04:05 1.41 @@ -168,6 +168,7 @@ #:decf #:float-digits #:rational + #:rationalize ) ;; Export types (:export #:qd-real @@ -237,6 +238,7 @@ #:decf #:float-digits #:rational + #:rationalize ) ;; Constants (:export #:+pi+

1 0

[oct-cvs] Oct commit: oct qd-class.lisp qd-methods.lisp
by rtoy 19 Sep '07

19 Sep '07

Update of /project/oct/cvsroot/oct In directory clnet:/tmp/cvs-serv18474 Modified Files: qd-class.lisp qd-methods.lisp Log Message: MAKE-QD should handle rationals better instead of converting them to doubles and then converting the qd-real. Convert the numerator and denominator to qd-real, and the divide. (This should be done better.) qd-class.lisp: o Change method to work on floats, instead of reals. qd-methods.lisp: o Add method to handle rationals. --- /project/oct/cvsroot/oct/qd-class.lisp 2007/09/06 02:58:38 1.25 +++ /project/oct/cvsroot/oct/qd-class.lisp 2007/09/19 17:30:04 1.26 @@ -60,7 +60,7 @@ (defmethod print-object ((qd qd-real) stream) (print-qd (qd-value qd) stream)) -(defmethod make-qd ((x real)) +(defmethod make-qd ((x cl:float)) (make-instance 'qd-real :value (make-qd-d (float x 1d0)))) (defmethod make-qd ((x qd-real)) --- /project/oct/cvsroot/oct/qd-methods.lisp 2007/09/18 12:46:36 1.58 +++ /project/oct/cvsroot/oct/qd-methods.lisp 2007/09/19 17:30:04 1.59 @@ -76,6 +76,16 @@ (defconstant +qd-real-one+ (make-instance 'qd-real :value (make-qd-d 1d0))) + +(defmethod make-qd ((x cl:rational)) + ;; We should do something better than this. + (let ((top (numerator x)) + (bot (denominator x))) + (make-instance 'qd-real + :value (div-qd (qdi::make-float (signum top) (abs top) 0 0 0) + (qdi::make-float (signum bot) (abs bot) 0 0 0))))) + + (defmethod add1 ((a number)) (cl::1+ a)) @@ -824,8 +834,10 @@ (cl:rational x1) (cl:rational x2) (cl:rational x3)))) + (define-compiler-macro + (&whole form &rest args) + (declare (ignore form)) (if (null args) 0 (do ((args (cdr args) (cdr args)) @@ -834,6 +846,7 @@ ((null args) res)))) (define-compiler-macro - (&whole form number &rest more-numbers) + (declare (ignore form)) (if more-numbers (do ((nlist more-numbers (cdr nlist)) (result number)) @@ -843,6 +856,7 @@ `(unary-minus ,number))) (define-compiler-macro * (&whole form &rest args) + (declare (ignore form)) (if (null args) 1 (do ((args (cdr args) (cdr args))

1 0

[oct-cvs] Oct commit: oct qd-complex.lisp qd-methods.lisp qd-package.lisp
by rtoy 18 Sep '07

18 Sep '07

Update of /project/oct/cvsroot/oct In directory clnet:/tmp/cvs-serv7806 Modified Files: qd-complex.lisp qd-methods.lisp qd-package.lisp Log Message: Add method RATIONAL to convert a quad-double to a rational. qd-package.lisp: o Appropriately shadow and export RATIONAL. o Need to export WITH-QD-PARTS from QDI. qd-methods.lisp: o Define methods for RATIONAL for reals and qd-reals. qd-complex.lisp: o Use CL:RATIONAL as appropriate for the CL rational type. --- /project/oct/cvsroot/oct/qd-complex.lisp 2007/09/12 21:01:13 1.36 +++ /project/oct/cvsroot/oct/qd-complex.lisp 2007/09/18 12:46:36 1.37 @@ -35,27 +35,27 @@ :real (sub-qd-d (qd-value (realpart a)) 1d0) :imag (qd-value (imagpart a)))) -(defmethod two-arg-/ ((a qd-real) (b rational)) +(defmethod two-arg-/ ((a qd-real) (b cl:rational)) (make-instance 'qd-real :value (div-qd (qd-value a) (qd-value (float b #q0))))) -(defmethod two-arg-/ ((a rational) (b qd-real)) +(defmethod two-arg-/ ((a cl:rational) (b qd-real)) (make-instance 'qd-real :value (div-qd (qd-value (float a #q0)) (qd-value b)))) -(defmethod two-arg-* ((a qd-real) (b rational)) +(defmethod two-arg-* ((a qd-real) (b cl:rational)) (make-instance 'qd-real :value (mul-qd (qd-value a) (qd-value (float b #q0))))) -(defmethod two-arg-+ ((a qd-real) (b rational)) +(defmethod two-arg-+ ((a qd-real) (b cl:rational)) (make-instance 'qd-real :value (add-qd (qd-value a) (qd-value (float b #q0))))) -(defmethod two-arg-+ ((a rational) (b qd-real)) +(defmethod two-arg-+ ((a cl:rational) (b qd-real)) (make-instance 'qd-real :value (add-qd (qd-value b) (qd-value (float a #q0))))) -(defmethod two-arg-- ((a qd-real) (b rational)) +(defmethod two-arg-- ((a qd-real) (b cl:rational)) (make-instance 'qd-real :value (sub-qd (qd-value a) (qd-value (float b #q0))))) -(defmethod two-arg-- ((a rational) (b qd-real)) +(defmethod two-arg-- ((a cl:rational) (b qd-real)) (make-instance 'qd-real :value (sub-qd (qd-value (float a #q0)) (qd-value b)))) (defmethod unary-minus ((z qd-complex)) --- /project/oct/cvsroot/oct/qd-methods.lisp 2007/09/16 02:39:29 1.57 +++ /project/oct/cvsroot/oct/qd-methods.lisp 2007/09/18 12:46:36 1.58 @@ -814,6 +814,16 @@ (defmethod float-digits ((x qd-real)) (* 4 (float-digits 1d0))) +(defmethod rational ((x real)) + (cl:rational x)) + +(defmethod rational ((x qd-real)) + (with-qd-parts (x0 x1 x2 x3) + (qd-value x) + (+ (cl:rational x0) + (cl:rational x1) + (cl:rational x2) + (cl:rational x3)))) (define-compiler-macro + (&whole form &rest args) (if (null args) --- /project/oct/cvsroot/oct/qd-package.lisp 2007/09/16 05:04:05 1.39 +++ /project/oct/cvsroot/oct/qd-package.lisp 2007/09/18 12:46:36 1.40 @@ -89,6 +89,7 @@ #:scale-float-qd #:ffloor-qd #:random-qd + #:with-qd-parts ) #+cmu (:export #:add-qd-dd @@ -166,6 +167,7 @@ #:incf #:decf #:float-digits + #:rational ) ;; Export types (:export #:qd-real @@ -234,6 +236,7 @@ #:incf #:decf #:float-digits + #:rational ) ;; Constants (:export #:+pi+

1 0

[oct-cvs] Oct commit: oct qd-rep.lisp qd.lisp
by rtoy 18 Sep '07

18 Sep '07

Update of /project/oct/cvsroot/oct In directory clnet:/tmp/cvs-serv19559 Modified Files: qd-rep.lisp qd.lisp Log Message: qd-rep.lisp: o Add macro WITH-QD-PARTS to extract the components of a quad-double. qd.lisp: o Use the macro as needed. --- /project/oct/cvsroot/oct/qd-rep.lisp 2007/09/17 03:07:27 1.6 +++ /project/oct/cvsroot/oct/qd-rep.lisp 2007/09/18 11:20:16 1.7 @@ -179,3 +179,13 @@ (aref qd 3))) ) ; end progn + +(defmacro with-qd-parts ((a0 a1 a2 a3) qd &body body) + (let ((q (gensym))) + `(let* ((,q ,qd) + (,a0 (qd-0 ,q)) + (,a1 (qd-1 ,q)) + (,a2 (qd-2 ,q)) + (,a3 (qd-3 ,q))) + ,@body))) + --- /project/oct/cvsroot/oct/qd.lisp 2007/09/17 19:04:23 1.53 +++ /project/oct/cvsroot/oct/qd.lisp 2007/09/18 11:20:16 1.54 @@ -393,10 +393,12 @@ ;; version? It's quite a bit more complicated. ;; ;; In addition, this is reorganized to minimize data dependency. - (multiple-value-bind (a0 a1 a2 a3) - (qd-parts a) - (multiple-value-bind (b0 b1 b2 b3) - (qd-parts b) + (with-qd-parts (a0 a1 a2 a3) + a + (declare (double-float a0 a1 a2 a3)) + (with-qd-parts (b0 b1 b2 b3) + b + (declare (double-float b0 b1 b2 b3)) (let ((s0 (cl:+ a0 b0)) (s1 (cl:+ a1 b1)) (s2 (cl:+ a2 b2)) @@ -411,18 +413,22 @@ (v1 (cl:- s1 a1)) (v2 (cl:- s2 a2)) (v3 (cl:- s3 a3))) + (declare (double-float v0 v1 v2 v3)) (let ((u0 (cl:- s0 v0)) (u1 (cl:- s1 v1)) (u2 (cl:- s2 v2)) (u3 (cl:- s3 v3))) + (declare (double-float u0 u1 u2 u3)) (let ((w0 (cl:- a0 u0)) (w1 (cl:- a1 u1)) (w2 (cl:- a2 u2)) (w3 (cl:- a3 u3))) + (declare (double-float w0 w1 w2 w3)) (let ((u0 (cl:- b0 v0)) (u1 (cl:- b1 v1)) (u2 (cl:- b2 v2)) (u3 (cl:- b3 v3))) + (declare (double-float u0 u1 u2 u3)) (let ((t0 (cl:+ w0 u0)) (t1 (cl:+ w1 u1)) (t2 (cl:+ w2 u2)) @@ -432,17 +438,18 @@ (three-sum s2 t0 t1 s2 t0 t1) (three-sum2 s3 t0 s3 t0 t2) (setf t0 (cl:+ t0 t1 t3)) - ;; Renormalize - (multiple-value-setq (s0 s1 s2 s3) - (renorm-5 s0 s1 s2 s3 t0)) - (if (and (zerop a0) (zerop b0)) - (%make-qd-d (+ a0 b0) 0d0 0d0 0d0) - (%make-qd-d s0 s1 s2 s3))))))))))) + ;; Renormalize + (multiple-value-setq (s0 s1 s2 s3) + (renorm-5 s0 s1 s2 s3 t0)) + (if (and (zerop a0) (zerop b0)) + (%make-qd-d (+ a0 b0) 0d0 0d0 0d0) + (%make-qd-d s0 s1 s2 s3))))))))))) (defun neg-qd (a) (declare (type %quad-double a)) - (multiple-value-bind (a0 a1 a2 a3) - (qd-parts a) + (with-qd-parts (a0 a1 a2 a3) + a + (declare (double-float a0 a1 a2 a3)) (%make-qd-d (cl:- a0) (cl:- a1) (cl:- a2) (cl:- a3)))) (defun sub-qd (a b) @@ -603,10 +610,12 @@ (space 0)) #+cmu (inline ext:float-infinity-p)) - (multiple-value-bind (a0 a1 a2 a3) - (qd-parts a) - (multiple-value-bind (b0 b1 b2 b3) - (qd-parts b) + (with-qd-parts (a0 a1 a2 a3) + a + (declare (double-float a0 a1 a2 a3)) + (with-qd-parts (b0 b1 b2 b3) + b + (declare (double-float b0 b1 b2 b3)) (multiple-value-bind (p0 q0) (two-prod a0 b0) #+cmu @@ -1124,8 +1133,9 @@ (k2 (- k k1)) (s1 (scale-float 1d0 k1)) (s2 (scale-float 1d0 k2))) - (multiple-value-bind (q0 q1 q2 q3) - (qd-parts qd) + (with-qd-parts (q0 q1 q2 q3) + qd + (declare (double-float q0 q1 q2 q3)) (%make-qd-d (cl:* (cl:* q0 s1) s2) (cl:* (cl:* q1 s1) s2) (cl:* (cl:* q2 s1) s2)

1 0

[oct-cvs] Oct commit: oct qd-test.lisp rt-tests.lisp
by rtoy 17 Sep '07

17 Sep '07

Update of /project/oct/cvsroot/oct In directory clnet:/tmp/cvs-serv24380 Modified Files: qd-test.lisp rt-tests.lisp Log Message: qd-test.lisp: o Add optional arg to enable/disable printing of results. Default is on. rt-tests.lisp: o Don't print results. --- /project/oct/cvsroot/oct/qd-test.lisp 2007/08/27 17:49:19 1.18 +++ /project/oct/cvsroot/oct/qd-test.lisp 2007/09/18 03:05:56 1.19 @@ -67,7 +67,7 @@ ;; pi = ;; 3.1415926535897932384626433832795028841971693993751058209749445923078L0 -(defun test2 () +(defun test2 (&optional (printp t)) ;; pi/4 = 4 * arctan(1/5) - arctan(1/239) ;; ;; Arctan is computed using the Taylor series @@ -100,11 +100,12 @@ (p (mul-qd-d (sub-qd (mul-qd-d s1 4d0) s2) 4d0))) - (format t "~2&pi via Machin's atan formula~%") - (print-result p +qd-pi+) + (when printp + (format t "~2&pi via Machin's atan formula~%") + (print-result p +qd-pi+)) p))) -(defun test3 () +(defun test3 (&optional (printp t)) (declare (optimize (speed 3))) ;; Salamin-Brent Quadratic formula for pi (let* ((a +qd-one+) @@ -127,11 +128,12 @@ (setf s s-new) (setf p (div-qd (mul-qd-d (sqr-qd a) 2d0) s)))) - (format t "~2&Salamin-Brent Quadratic formula for pi~%") - (print-result p +qd-pi+) + (when printp + (format t "~2&Salamin-Brent Quadratic formula for pi~%") + (print-result p +qd-pi+)) p)) -(defun test4 () +(defun test4 (&optional (printp t)) (declare (optimize (speed 3))) ;; Borwein Quartic formula for pi (let* ((a (sub-qd (make-qd-d 6d0) @@ -160,13 +162,14 @@ m))) (setf p (div-qd +qd-one+ a)))) - (format t "~2&Borwein's Quartic formula for pi~%") - (print-result p +qd-pi+) + (when printp + (format t "~2&Borwein's Quartic formula for pi~%") + (print-result p +qd-pi+)) p)) ;; e = ;; 2.718281828459045235360287471352662497757247093699959574966967627724L0 -(defun test5 () +(defun test5 (&optional (printp t)) ;; Taylor series for e (let ((s (make-qd-d 2d0)) (tmp +qd-one+) @@ -179,13 +182,14 @@ (setf tmp (div-qd tmp (make-qd-d (float n 1d0)))) (setf s (add-qd s tmp))) - (format t "~2&e via Taylor series~%") - (print-result s +qd-e+) + (when printp + (format t "~2&e via Taylor series~%") + (print-result s +qd-e+)) s)) ;; log(2) = ;; 0.6931471805599453094172321214581765680755001343602552541206800094934L0 -(defun test6 () +(defun test6 (&optional (printp t)) ;; Taylor series for log 2 ;; ;; -log(1-x) = x + x^2/2 + x^3/3 + x^4/4 + ... @@ -201,11 +205,12 @@ (setf tt (mul-qd-d tt .5d0)) (setf s (add-qd s (div-qd tt (make-qd-d (float n 1d0)))))) - (format t "~2&log(2) via Taylor series~%") - (print-result s +qd-log2+) + (when printp + (format t "~2&log(2) via Taylor series~%") + (print-result s +qd-log2+)) s)) -(defun test-atan (&optional (fun #'atan-qd)) +(defun test-atan (&optional (fun #'atan-qd) (printp t)) ;; Compute atan for known values (format t "~2&atan via ~S~%" fun) @@ -234,7 +239,7 @@ (format t "bits = ~,1f~%" (bit-accuracy y true)))) -(defun test-sin () +(defun test-sin (&optional (printp t)) (format t "~2&sin~%") (let* ((arg (div-qd +qd-pi+ (make-qd-d 6d0))) (y (sin-qd arg)) @@ -259,7 +264,7 @@ (bit-accuracy y true))) ) -(defun test-tan (&optional (f #'tan-qd)) +(defun test-tan (&optional (f #'tan-qd) (printp t)) (format t "~2&tan via ~S~%" f) (let* ((arg (div-qd +qd-pi+ (make-qd-d 6d0))) (y (funcall f arg)) @@ -284,7 +289,7 @@ (bit-accuracy y true))) ) -(defun test-asin () +(defun test-asin (&optional (printp t)) (format t "~2&asin~%") (let* ((arg (make-qd-d 0.5d0)) (y (asin-qd arg)) @@ -309,7 +314,7 @@ (bit-accuracy y true))) ) -(defun test-log (f) +(defun test-log (f &optional (printp t)) (format t "~2&Log via ~A~%" f) (let* ((arg (make-qd-d 2d0)) (y (funcall f arg)) @@ -334,7 +339,7 @@ (bit-accuracy y true))) ) -(defun test-log1p (f) +(defun test-log1p (f &optional (printp t)) (format t "~2&Log1p via ~A~%" f) (let* ((arg (make-qd-d 1d0)) (y (funcall f arg)) @@ -359,7 +364,7 @@ (bit-accuracy y true))) ) -(defun test-exp (f) +(defun test-exp (f &optional (printp t)) (format t "~2&Exp via ~A~%" f) (let* ((arg +qd-zero+) (y (funcall f arg)) @@ -385,6 +390,7 @@ (bit-accuracy y true))) ) + (defun all-tests () (test2) (test3) --- /project/oct/cvsroot/oct/rt-tests.lisp 2007/08/27 18:05:12 1.1 +++ /project/oct/cvsroot/oct/rt-tests.lisp 2007/09/18 03:05:56 1.2 @@ -50,7 +50,7 @@ ;; Pi via Machin's formula (rt:deftest oct.pi.machin (let* ((*standard-output* *null*) - (val (make-instance 'qd-real :value (qdi::test2))) + (val (make-instance 'qd-real :value (qdi::test2 nil))) (true qd:+pi+)) (check-accuracy 213 val true)) nil) @@ -58,7 +58,7 @@ ;; Pi via Salamin-Brent algorithm (rt:deftest oct.pi.salamin-brent (let* ((*standard-output* *null*) - (val (make-instance 'qd-real :value (qdi::test3))) + (val (make-instance 'qd-real :value (qdi::test3 nil))) (true qd:+pi+)) (check-accuracy 202 val true)) nil) @@ -66,7 +66,7 @@ ;; Pi via Borweign's Quartic formula (rt:deftest oct.pi.borweign (let* ((*standard-output* *null*) - (val (make-instance 'qd-real :value (qdi::test4))) + (val (make-instance 'qd-real :value (qdi::test4 nil))) (true qd:+pi+)) (check-accuracy 211 val true)) nil) @@ -74,7 +74,7 @@ ;; e via Taylor series (rt:deftest oct.e.taylor (let* ((*standard-output* *null*) - (val (make-instance 'qd-real :value (qdi::test5))) + (val (make-instance 'qd-real :value (qdi::test5 nil))) (true (make-instance 'qd-real :value qdi::+qd-e+))) (check-accuracy 212 val true)) nil) @@ -82,7 +82,7 @@ ;; log(2) via Taylor series (rt:deftest oct.log2.taylor (let* ((*standard-output* *null*) - (val (make-instance 'qd-real :value (qdi::test6))) + (val (make-instance 'qd-real :value (qdi::test6 nil))) (true (make-instance 'qd-real :value qdi::+qd-log2+))) (check-accuracy 212 val true)) nil)

1 0

[oct-cvs] Oct commit: oct qd.lisp
by rtoy 17 Sep '07

17 Sep '07

Update of /project/oct/cvsroot/oct In directory clnet:/tmp/cvs-serv8007 Modified Files: qd.lisp Log Message: New SCALE-FLOAT-QD implementation that shouldn't suffer from premature overflow/underflow. (Still issue if the exponent is very large or very small, though, but not if the exponent is < 2000 or so.) --- /project/oct/cvsroot/oct/qd.lisp 2007/09/17 17:15:04 1.52 +++ /project/oct/cvsroot/oct/qd.lisp 2007/09/17 19:04:23 1.53 @@ -1085,6 +1085,7 @@ q0-sign)))))) (declaim (inline scale-float-qd)) +#+(or) (defun scale-float-qd (qd k) (declare (type %quad-double qd) (type fixnum k) @@ -1112,6 +1113,29 @@ (cl:* (qd-2 qd) scale) (cl:* (qd-3 qd) scale)))) +(defun scale-float-qd (qd k) + (declare (type %quad-double qd) + (type (integer -2000 2000) k) + (optimize (speed 3) (space 0))) + ;; Split the exponent in half and do the scaling in two parts. + ;; Requires 2 multiplications, but should not prematurely return 0, + ;; and should be faster than the original version above. + (let* ((k1 (floor k 2)) + (k2 (- k k1)) + (s1 (scale-float 1d0 k1)) + (s2 (scale-float 1d0 k2))) + (multiple-value-bind (q0 q1 q2 q3) + (qd-parts qd) + (%make-qd-d (cl:* (cl:* q0 s1) s2) + (cl:* (cl:* q1 s1) s2) + (cl:* (cl:* q2 s1) s2) + (cl:* (cl:* q3 s1) s2))))) + + + + + + (defun decode-float-qd (q) (declare (type %quad-double q)) (multiple-value-bind (frac exp sign)

1 0

[oct-cvs] Oct commit: oct qd.lisp
by rtoy 17 Sep '07

17 Sep '07

Update of /project/oct/cvsroot/oct In directory clnet:/tmp/cvs-serv21126 Modified Files: qd.lisp Log Message: o Convert THREE-SUM2 to a macro instead of a function (to speed things up for Allegro and other Lisps that don't inline). o Update code for the THREE-SUM2 macro. --- /project/oct/cvsroot/oct/qd.lisp 2007/09/17 14:06:20 1.51 +++ /project/oct/cvsroot/oct/qd.lisp 2007/09/17 17:15:04 1.52 @@ -43,7 +43,7 @@ (setf *inline-expansion-limit* 1600)) ;; All of the following functions should be inline. -(declaim (inline three-sum2)) +;;(declaim (inline three-sum2)) ;; Internal routines for implementing quad-double. @@ -73,6 +73,7 @@ +#+nil (defun three-sum2 (a b c) (declare (double-float a b c) (optimize (speed 3))) @@ -81,8 +82,25 @@ (t3 t1)) (two-sum t1 t2 a b) (two-sum a t3 c t1) - (values a (cl:+ t2 t3) c))) + (values a (cl:+ t2 t3)))) +(defmacro three-sum2 (s0 s1 a b c) + (let ((aa (gensym)) + (bb (gensym)) + (cc (gensym)) + (t1 (gensym)) + (t2 (gensym)) + (t3 (gensym))) + `(let* ((,aa ,a) + (,bb ,b) + (,cc ,c) + (,t1 0d0) + (,t2 ,t1) + (,t3 ,t1)) + (declare (double-float ,t1 ,t2 ,t3)) + (two-sum ,t1 ,t2 ,aa ,bb) + (two-sum ,s0 ,t3 ,cc ,t1) + (setf ,s1 (+ ,t2 ,t3))))) ;; Not needed???? #+nil @@ -409,18 +427,17 @@ (t1 (cl:+ w1 u1)) (t2 (cl:+ w2 u2)) (t3 (cl:+ w3 u3))) + (declare (double-float t0 t1 t2 t3)) (two-sum s1 t0 s1 t0) (three-sum s2 t0 t1 s2 t0 t1) - (multiple-value-bind (s3 t0) - (three-sum2 s3 t0 t2) - (declare (double-float t0)) - (setf t0 (cl:+ t0 t1 t3)) + (three-sum2 s3 t0 s3 t0 t2) + (setf t0 (cl:+ t0 t1 t3)) ;; Renormalize (multiple-value-setq (s0 s1 s2 s3) (renorm-5 s0 s1 s2 s3 t0)) (if (and (zerop a0) (zerop b0)) (%make-qd-d (+ a0 b0) 0d0 0d0 0d0) - (%make-qd-d s0 s1 s2 s3)))))))))))) + (%make-qd-d s0 s1 s2 s3))))))))))) (defun neg-qd (a) (declare (type %quad-double a)) @@ -472,8 +489,10 @@ (return-from mul-qd-d (%make-qd-d p0 0d0 0d0 0d0))) (multiple-value-bind (p1 q1) (two-prod (qd-1 a) b) + (declare (double-float p1 q1)) (multiple-value-bind (p2 q2) (two-prod (qd-2 a) b) + (declare (double-float p2 q2)) (let* ((p3 (cl:* (qd-3 a) b)) (s0 p0) (s1 p0) @@ -481,15 +500,14 @@ (declare (double-float s0 s1 s2 p3)) (two-sum s1 s2 q0 p1) (three-sum s2 q1 p2 s2 q1 p2) - (multiple-value-bind (q1 q2) - (three-sum2 q1 q2 p3) - (let ((s3 q1) - (s4 (cl:+ q2 p2))) - (multiple-value-bind (s0 s1 s2 s3) - (renorm-5 s0 s1 s2 s3 s4) - (if (zerop s0) - (%make-qd-d (float-sign p0 0d0) 0d0 0d0 0d0) - (%make-qd-d s0 s1 s2 s3)))))))))) + (three-sum2 q1 q2 q1 q2 p3) + (let ((s3 q1) + (s4 (cl:+ q2 p2))) + (multiple-value-bind (s0 s1 s2 s3) + (renorm-5 s0 s1 s2 s3 s4) + (if (zerop s0) + (%make-qd-d (float-sign p0 0d0) 0d0 0d0 0d0) + (%make-qd-d s0 s1 s2 s3))))))))) ;; a0 * b0 0 ;; a0 * b1 1 @@ -556,7 +574,7 @@ (cl:* (qd-3 a) (kernel:double-double-lo b)) q3 q4))) - (multiple-value-setq (p3 q0 s1) + (multiple-value-setq (p3 q0) (three-sum2 p3 q0 s1)) (let ((p4 (cl:+ q0 s2))) (multiple-value-call #'%make-qd-d

1 0