Update of /project/oct/cvsroot/oct In directory clnet:/tmp/cvs-serv22733
Modified Files: qd-io.lisp qd-methods.lisp qd-package.lisp Log Message: qd-io.lisp: o Add RATIONAL-TO-QD, a simple, fast and accurate method to convert rationals to quad-doubles. (From Richard Fateman.) o Use RATIONAL-TO-QD to create a quad-float
qd-methods.lisp: o Use RATIONAL-TO-QD to create a quad-float from a bignum and ratio.
qd-package.lisp: o Export RATIONAL-TO-QD
--- /project/oct/cvsroot/oct/qd-io.lisp 2007/09/24 21:32:15 1.17 +++ /project/oct/cvsroot/oct/qd-io.lisp 2007/10/10 15:24:07 1.18 @@ -381,78 +381,29 @@ (neg-qd (mul-qd xx yy)) (mul-qd xx yy))))))))
+;; This is a slightly modified version of Richard Fateman's code to +;; convert bignums to qd. This supports converting rationals to qd +;; too. +(defun rational-to-qd (rat) + (declare (rational rat)) + (let* ((p (coerce rat 'double-float)) + (ans (make-qd-d p)) + (remainder rat)) + (declare (double-float p) + (rational remainder) + (type %quad-double ans)) + (dotimes (k 3 ans) + (decf remainder (rational p)) + (setf p (coerce remainder 'double-float)) + (setf ans (add-qd-d ans p))))) + (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))))) + (rational-to-qd (* sign + (* (+ int-part (/ frac-part (expt 10 scale))) + (expt 10 exp)))))
;; This seems to work, but really needs to be rewritten! --- /project/oct/cvsroot/oct/qd-methods.lisp 2007/09/24 21:30:07 1.61 +++ /project/oct/cvsroot/oct/qd-methods.lisp 2007/10/10 15:24:07 1.62 @@ -250,11 +250,7 @@
(defun bignum-to-qd (bignum) (make-instance 'qd-real - :value (qdi::make-float (if (minusp bignum) -1 1) - (abs bignum) - 0 - 0 - 0))) + :value (rational-to-qd bignum)))
(defmethod qfloat ((x real) (num-type cl:float)) (cl:float x num-type)) @@ -276,10 +272,7 @@ (qfloat (denominator x) num-type)))
(defmethod qfloat ((x ratio) (num-type qd-real)) - ;; This probably has some issues with roundoff - (let ((top (qd-value (qfloat (numerator x) num-type))) - (bot (qd-value (qfloat (denominator x) num-type)))) - (make-instance 'qd-real :value (div-qd top bot)))) + (make-instance 'qd-real :value (rational-to-qd x)))
#+cmu (defmethod qfloat ((x ext:double-double-float) (num-type qd-real)) @@ -1025,4 +1018,3 @@ ;; and make a real qd-real float, instead of the hackish ;; %qd-real. (set-dispatch-macro-character ## #\Q #'qd-class-reader) - --- /project/oct/cvsroot/oct/qd-package.lisp 2007/09/20 21:04:05 1.41 +++ /project/oct/cvsroot/oct/qd-package.lisp 2007/10/10 15:24:07 1.42 @@ -90,6 +90,7 @@ #:ffloor-qd #:random-qd #:with-qd-parts + #:rational-to-qd ) #+cmu (:export #:add-qd-dd