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)