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

[rtoy]

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