[oct-cvs] Oct commit: oct qd-io.lisp qd-methods.lisp qd-package.lisp

rtoy rtoy at common-lisp.net
Wed Oct 10 15:24:07 UTC 2007

```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