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

[oct-cvs] Oct commit: oct timing2.lisp
by rtoy 11 Oct '07

11 Oct '07

Update of /project/oct/cvsroot/oct In directory clnet:/tmp/cvs-serv20797 Added Files: timing2.lisp Log Message: Initial version. --- /project/oct/cvsroot/oct/timing2.lisp 2007/10/11 16:29:08 NONE +++ /project/oct/cvsroot/oct/timing2.lisp 2007/10/11 16:29:08 1.1 ;;;; -*- Mode: lisp -*- ;;;; ;;;; Copyright (c) 2007 Raymond Toy ;;;; ;;;; Permission is hereby granted, free of charge, to any person ;;;; obtaining a copy of this software and associated documentation ;;;; files (the "Software"), to deal in the Software without ;;;; restriction, including without limitation the rights to use, ;;;; copy, modify, merge, publish, distribute, sublicense, and/or sell ;;;; copies of the Software, and to permit persons to whom the ;;;; Software is furnished to do so, subject to the following ;;;; conditions: ;;;; ;;;; The above copyright notice and this permission notice shall be ;;;; included in all copies or substantial portions of the Software. ;;;; ;;;; THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, ;;;; EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES ;;;; OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND ;;;; NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT ;;;; HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, ;;;; WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING ;;;; FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR ;;;; OTHER DEALINGS IN THE SOFTWARE. ;;; Here are some simple timing tests, based on Yozo Hida's qd_timer ;;; test code. I've tried to make these versions time the same ;;; operations as Yozo's. (in-package #:qd) (defun time-add (&optional (n 100000)) (declare (fixnum n) (optimize (speed 3))) (let ((a (/ #q7)) (b #q0)) (declare (type qd-real a b)) (time (dotimes (k n) (declare (fixnum k)) (incf b a))) (format t "n = ~D~%" n) (format t "b = ~W~%" b))) (defun time-mul (&optional (n 100000)) (declare (fixnum n) (optimize (speed 3))) (let ((a (+ 1 (/ (float n #q1)))) (b #q1)) (declare (type qd-real a b)) (time (dotimes (k n) (declare (fixnum k)) (setf b (* a b)))) (format t "n = ~D~%" n) (format t "b = ~W~%" b))) (defun time-mul (&optional (n 100000)) (declare (fixnum n) (optimize (speed 3))) (let ((a (+ 1 (/ (float n #q1)))) (b #q1)) (declare (type qd-real a b)) (time (dotimes (k n) (declare (fixnum k)) (setf b (* a b)))) (format t "n = ~D~%" n) (format t "b = ~W~%" b))) (defun time-div (&optional (n 100000)) (declare (fixnum n) (optimize (speed 3))) (let ((a (+ 1 (/ (float n #q1)))) (b #q1)) (declare (type qd-real a b)) (time (dotimes (k n) (declare (fixnum k)) (setf b (/ b a)))) (format t "n = ~D~%" n) (format t "b = ~W~%" b))) (defun time-sqrt (&optional (n 100000)) (declare (fixnum n) (optimize (speed 3))) (let ((a #q0) (b (+ 2 +pi+))) (declare (type qd-real a b)) (time (dotimes (k n) (declare (fixnum k)) (setf a (sqrt (+ a b))))) (format t "n = ~D~%" n) (format t "a = ~W~%" a))) (defun time-sin (&optional (n 2000)) (declare (fixnum n) (optimize (speed 3))) (let ((a #q0) (b (/ +pi+ n)) (c #q0)) (declare (type qd-real a b c)) (time (dotimes (k n) (declare (fixnum k)) (incf a b) (incf c (sin a)))) (format t "n = ~D~%" n) (format t "a = ~W~%" a))) (defun time-log (&optional (n 1000)) (declare (fixnum n) (optimize (speed 3))) (let ((a #q0) (c (exp #q-50.1)) (d (exp (/ #q100.2 n)))) (declare (type qd-real a c d)) (time (dotimes (k n) (declare (fixnum k)) (incf a (log c)) (setf c (* c d)))) (format t "n = ~D~%" n) (format t "a = ~W~%" a))) #|| Some test results. These were all run on a Sun Blade 1500 using a 1.5 GHz Ultrasparc III. I used the default configuration when compiling qd, and used Sun's C++ compiler. For the Lisp timing, I used CMUCL. Executive summary: Test Time qd oct ---- ----------- add 0.036 0.09 mul 0.117 0.13 div 0.388 0.29 sqrt 0.142 0.11 sin 0.136 0.14 log 0.231 0.12 Times are in sec for the test. The default number of iterations were used. Some of the results are a bit surprising. I was expecting the C++ code to be faster, and that's the case for add and mul. But oct with CMUCL was faster in div, sqrt, and log. Note that oct uses a different algorithm for log than qd, so that could explain the difference. Note also that CLOS dispatch is included in these timing results. CMUCL's CLOS dispatch is ok, but not great. ------------------------------------------------------------------------------- The raw data: The output from qd_timer -qd -v: Timing qd_real -------------- Timing addition... n = 100000 t = 0.0362288 b = 1.428571e+04 100000 operations in 0.0362288 s. 0.362288 us Timing multiplication ... n = 100000 t = 0.11686 b = 2.718268e+00 100000 operations in 0.11686 s. 1.168602 us Timing division ... n = 100000 t = 0.388279 b = 0.367881 100000 operations in 0.388279 s. 3.882788 us Timing square root ... n = 10000 t = 0.141866 a = 2.821980 10000 operations in 0.141866 s. 14.186575 us Timing sin ... n = 2000 t = 0.136080 a = 3.141593 2000 operations in 0.136080 s. 68.039791 us Timing log ... n = 1000 t = 0.230506 a = -50.100000 1000 operations in 0.230506 s. 230.506166 us The output from CMUCL: QD> (time-add) ; Evaluation took: ; 0.09 seconds of real time ; 0.1 seconds of user run time ; 0.0 seconds of system run time ; 147,285,856 CPU cycles ; 0 page faults and ; 7,200,016 bytes consed. ; n = 100000 b = #q14285.7142857142857142857142857142857142857142857142857142857142855q0 NIL QD> (time-mul) ; Evaluation took: ; 0.13 seconds of real time ; 0.1 seconds of user run time ; 0.02 seconds of system run time ; 203,790,588 CPU cycles ; 0 page faults and ; 7,200,824 bytes consed. ; n = 100000 b = #q2.71826823717448966803506482442604644797444693267782286300915989397q0 NIL QD> (time-div) ; Evaluation took: ; 0.29 seconds of real time ; 0.28 seconds of user run time ; 0.01 seconds of system run time ; 460,956,912 CPU cycles ; 0 page faults and ; 7,200,016 bytes consed. ; n = 100000 b = #q0.36788128056098406210328658773118942247132502490133718973918140856q0 NIL QD> (time-sqrt 10000) ; Evaluation took: ; 0.11 seconds of real time ; 0.1 seconds of user run time ; 0.0 seconds of system run time ; 173,209,708 CPU cycles ; 0 page faults and ; 2,402,560 bytes consed. ; n = 10000 a = #q2.82198033014704783016853125515542796898998765943212617578596649019q0 NIL QD> (time-sin) ; Evaluation took: ; 0.14 seconds of real time ; 0.14 seconds of user run time ; 0.0 seconds of system run time ; 213,378,476 CPU cycles ; 0 page faults and ; 3,105,800 bytes consed. ; n = 2000 a = #q3.14159265358979323846264338327950288419716939937510582097494459409q0 NIL QD> (time-log) ; Evaluation took: ; 0.12 seconds of real time ; 0.12 seconds of user run time ; 0.01 seconds of system run time ; 192,187,304 CPU cycles ; 0 page faults and ; 1,621,792 bytes consed. ; n = 1000 a = #q-50.100000000000000000000000000000000000000000000000000000000208796q0 NIL QD> --------------------------------------------- ||#

1 0

[oct-cvs] Oct commit: oct .cvsignore
by rtoy 10 Oct '07

10 Oct '07

Update of /project/oct/cvsroot/oct In directory clnet:/tmp/cvs-serv22969 Modified Files: .cvsignore Log Message: Ignore *.lib and *.fas too. (Files from clisp) --- /project/oct/cvsroot/oct/.cvsignore 2007/09/12 02:03:42 1.3 +++ /project/oct/cvsroot/oct/.cvsignore 2007/10/10 15:25:41 1.4 @@ -3,3 +3,5 @@ *.x86f *.err *.fasl +*.lib +*.fas

1 0

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

10 Oct '07

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

1 0

[oct-cvs] Oct commit: oct qd-fun.lisp
by rtoy 10 Oct '07

10 Oct '07

Update of /project/oct/cvsroot/oct In directory clnet:/tmp/cvs-serv22545 Modified Files: qd-fun.lisp Log Message: o Move logb-finite before first use. o Declare arg as double-float, not float. --- /project/oct/cvsroot/oct/qd-fun.lisp 2007/09/12 02:31:14 1.80 +++ /project/oct/cvsroot/oct/qd-fun.lisp 2007/10/10 15:21:47 1.81 @@ -34,6 +34,17 @@ (in-package #:qdi) +(defun logb-finite (x) + "Same as logb but X is not infinity and non-zero and not a NaN, so +that we can always return an integer" + (declare (type cl:double-float x)) + (multiple-value-bind (signif expon sign) + (cl:decode-float x) + (declare (ignore signif sign)) + ;; decode-float is almost right, except that the exponent + ;; is off by one + (1- expon))) + #+cmu (declaim (maybe-inline sqrt-qd)) (defun sqrt-qd (a) @@ -69,17 +80,6 @@ (setf r (add-qd r (mul-qd r (sub-d-qd half (mul-qd h (sqr-qd r))))))) (scale-float-qd (mul-qd r new-a) (ash k -1))))) -(defun logb-finite (x) - "Same as logb but X is not infinity and non-zero and not a NaN, so -that we can always return an integer" - (declare (type cl:float x)) - (multiple-value-bind (signif expon sign) - (cl:decode-float x) - (declare (ignore signif sign)) - ;; decode-float is almost right, except that the exponent - ;; is off by one - (1- expon))) - (defun hypot-aux-qd (x y) (declare (type %quad-double x y)) (let ((k (- (logb-finite (max (cl:abs (qd-0 x))

1 0

[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 24 Sep '07

24 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