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

[oct-cvs] Oct commit: oct oct.asd oct.system
by rtoy 16 Sep '07

16 Sep '07

Update of /project/oct/cvsroot/oct In directory clnet:/tmp/cvs-serv5887 Modified Files: oct.asd oct.system Log Message: Update dependencies. --- /project/oct/cvsroot/oct/oct.asd 2007/09/13 01:07:04 1.1 +++ /project/oct/cvsroot/oct/oct.asd 2007/09/16 05:00:00 1.2 @@ -62,9 +62,9 @@ ((:file "qd-package") (:file "qd-rep" :depends-on ("qd-package")) #-cmu - (:file "qd-dd" :depends-on ("qd-package")) + (:file "qd-dd" :depends-on ("qd-package" "qd-rep")) (:file "qd" - :depends-on ("qd-rep")) + :depends-on ("qd-rep" #-cmu "qd-dd")) (:file "qd-io" :depends-on ("qd")) (:file "qd-const" --- /project/oct/cvsroot/oct/oct.system 2007/08/27 17:51:02 1.21 +++ /project/oct/cvsroot/oct/oct.system 2007/09/16 05:00:00 1.22 @@ -50,9 +50,9 @@ ((:file "qd-package") (:file "qd-rep" :depends-on ("qd-package")) #-cmu - (:file "qd-dd" :depends-on ("qd-package")) + (:file "qd-dd" :depends-on ("qd-package" "qd-rep")) (:file "qd" - :depends-on ("qd-rep")) + :depends-on ("qd-rep" #-cmu "qd-dd")) (:file "qd-io" :depends-on ("qd")) (:file "qd-const"

1 0

[oct-cvs] Oct commit: oct qd-dd.lisp qd-package.lisp qd.lisp
by rtoy 16 Sep '07

16 Sep '07

Update of /project/oct/cvsroot/oct In directory clnet:/tmp/cvs-serv12645 Modified Files: qd-dd.lisp qd-package.lisp qd.lisp Log Message: To speed up Allegro (and other Lisp's that don't support inline functions), change QUICK-TWO-SUM from a function to a macro. Note that macro has different calling convention than the function. This is needed because Allegro apparently doesn't handle VALUES without boxing. All rt tests pass. qd-package.lisp: o For CMUCL, don't import C::QUICK-TWO-SUM into the QDI package anymore. qd-dd.lisp: o New QUICK-TWO-SUM macro. qd.lisp: o Add CMUCL version of QUICK-TWO-SUM macro, which just calls C::QUICK-TWO-SUM. o Update all users of QUICK-TWO-SUM appropriately. --- /project/oct/cvsroot/oct/qd-dd.lisp 2007/09/13 17:28:30 1.8 +++ /project/oct/cvsroot/oct/qd-dd.lisp 2007/09/16 02:46:24 1.9 @@ -31,6 +31,7 @@ ;;; ;;; These routines were taken directly from CMUCL. +#|| (declaim (inline quick-two-sum)) (defun quick-two-sum (a b) "Computes fl(a+b) and err(a+b), assuming |a| >= |b|" @@ -40,6 +41,16 @@ (e (- b (- s a)))) (declare (double-float s e)) (values s e))) +||# + +(defmacro quick-two-sum (s e x y) + (let ((a (gensym)) + (b (gensym))) + `(let* ((,a ,x) + (,b ,y)) + (declare (double-float ,a ,b ,s ,e)) + (setf ,s (+ ,a ,b)) + (setf ,e (- ,b (- ,s ,a)))))) (declaim (inline two-sum)) (defun two-sum (a b) @@ -53,6 +64,21 @@ (declare (double-float s v e)) (values s e))) +#+nil +(defmacro two-sum (s e x y) + "Computes fl(a+b) and err(a+b)" + (let ((a (gensym)) + (b (gensym)) + (v (gensym)) + `(let ((,a ,x) + (,b ,y)) + (declare (double-float ,a ,b)) + (setf ,s (+ ,a ,b)) + (let ((,v (- ,s ,a))) + (declare (double-float v)) + (setf e (+ (- ,a (- ,s ,v)) + (- ,b ,v)))))))) + (declaim (inline two-prod)) (declaim (inline split)) ;; This algorithm is the version given by Yozo Hida. It has problems --- /project/oct/cvsroot/oct/qd-package.lisp 2007/09/12 21:01:13 1.37 +++ /project/oct/cvsroot/oct/qd-package.lisp 2007/09/16 02:46:24 1.38 @@ -98,7 +98,6 @@ #+cmu (:import-from #:c #:two-sum - #:quick-two-sum #:two-prod #:two-sqr)) --- /project/oct/cvsroot/oct/qd.lisp 2007/09/12 02:31:14 1.46 +++ /project/oct/cvsroot/oct/qd.lisp 2007/09/16 02:46:24 1.47 @@ -132,6 +132,16 @@ div-qd-d div-qd-dd)) +#+cmu +(defmacro quick-two-sum (s e x y) + `(multiple-value-setq (,s ,e) + (c::quick-two-sum ,x ,y))) + +#+(and nil cmu) +(defmacro two-sum (s e x y) + `(multiple-value-setq (s e) + (c::two-sum x y))) + #-(or qd-inline (not cmu)) (declaim (ext:start-block renorm-4 renorm-5 make-qd-d @@ -172,91 +182,67 @@ (defun renorm-4 (c0 c1 c2 c3) (declare (double-float c0 c1 c2 c3) - (optimize (speed 3) (safety 0))) - (let ((s2 0d0) + (optimize (speed 3) (safety 0) (debug 0))) + (let ((s0 0d0) + (s1 0d0) + (s2 0d0) (s3 0d0)) - (multiple-value-bind (s0 c3) - (quick-two-sum c2 c3) - (multiple-value-bind (s0 c2) - (quick-two-sum c1 s0) - (multiple-value-bind (c0 c1) - (quick-two-sum c0 s0) - (let ((s0 c0) - (s1 c1)) - (cond ((/= s1 0) - (multiple-value-setq (s1 s2) - (quick-two-sum s1 c2)) - (if (/= s2 0) - (multiple-value-setq (s2 s3) - (quick-two-sum s2 c3)) - (multiple-value-setq (s1 s2) - (quick-two-sum s1 c3)))) - (t - (multiple-value-setq (s0 s1) - (quick-two-sum s0 c2)) - (if (/= s1 0) - (multiple-value-setq (s1 s2) - (quick-two-sum s1 c3)) - (multiple-value-setq (s0 s1) - (quick-two-sum s0 c3))))) - (values s0 s1 s2 s3))))))) - + (declare (double-float s0 s1 s2 s3)) + (quick-two-sum s0 c3 c2 c3) + (quick-two-sum s0 c2 c1 s0) + (quick-two-sum c0 c1 c0 s0) + (setf s0 c0) + (setf s1 c1) + (cond ((/= s1 0) + (quick-two-sum s1 s2 s1 c2) + (if (/= s2 0) + (quick-two-sum s2 s3 s2 c3) + (quick-two-sum s1 s2 s1 c3))) + (t + (quick-two-sum s0 s1 s0 c2) + (if (/= s1 0) + (quick-two-sum s1 s2 s1 c3) + (quick-two-sum s0 s1 s0 c3)))) + (values s0 s1 s2 s3))) + (defun renorm-5 (c0 c1 c2 c3 c4) - (declare (double-float c0 c1 c2 c3) + (declare (double-float c0 c1 c2 c3 c4) (optimize (speed 3) (safety 0))) - (let ((s2 0d0) + (let ((s0 0d0) + (s1 0d0) + (s2 0d0) (s3 0d0)) - (declare (double-float s2 s3)) - (multiple-value-bind (s0 c4) - (quick-two-sum c3 c4) - (multiple-value-bind (s0 c3) - (quick-two-sum c2 s0) - (multiple-value-bind (s0 c2) - (quick-two-sum c1 s0) - (multiple-value-bind (c0 c1) - (quick-two-sum c0 s0) - (let ((s0 c0) - (s1 c1)) - (declare (double-float s0 s1)) - (multiple-value-setq (s0 s1) - (quick-two-sum c0 c1)) - (cond ((/= s1 0) - (multiple-value-setq (s1 s2) - (quick-two-sum s1 c2)) - (cond ((/= s2 0) - (multiple-value-setq (s2 s3) - (quick-two-sum s2 c3)) - (if (/= s3 0) - (incf s3 c4) - (incf s2 c4))) - (t - (multiple-value-setq (s1 s2) - (quick-two-sum s1 c3)) - (if (/= s2 0) - (multiple-value-setq (s2 s3) - (quick-two-sum s2 c4)) - (multiple-value-setq (s1 s2) - (quick-two-sum s1 c4)))))) - (t - (multiple-value-setq (s0 s1) - (quick-two-sum s0 c2)) - (cond ((/= s1 0) - (multiple-value-setq (s1 s2) - (quick-two-sum s1 c3)) - (if (/= s2 0) - (multiple-value-setq (s2 s3) - (quick-two-sum s2 c4)) - (multiple-value-setq (s1 s2) - (quick-two-sum s1 c4)))) - (t - (multiple-value-setq (s0 s1) - (quick-two-sum s0 c3)) - (if (/= s1 0) - (multiple-value-setq (s1 s2) - (quick-two-sum s1 c4)) - (multiple-value-setq (s0 s1) - (quick-two-sum s0 c4))))))) - (values s0 s1 s2 s3)))))))) + (declare (double-float s0 s1 s2 s3)) + (quick-two-sum s0 c4 c3 c4) + (quick-two-sum s0 c3 c2 s0) + (quick-two-sum s0 c2 c1 s0) + (quick-two-sum c0 c1 c0 s0) + (quick-two-sum s0 s1 c0 c1) + (cond ((/= s1 0) + (quick-two-sum s1 s2 s1 c2) + (cond ((/= s2 0) + (quick-two-sum s2 s3 s2 c3) + (if (/= s3 0) + (incf s3 c4) + (incf s2 c4))) + (t + (quick-two-sum s1 s2 s1 c3) + (if (/= s2 0) + (quick-two-sum s2 s3 s2 c4) + (quick-two-sum s1 s2 s1 c4))))) + (t + (quick-two-sum s0 s1 s0 c2) + (cond ((/= s1 0) + (quick-two-sum s1 s2 s1 c3) + (if (/= s2 0) + (quick-two-sum s2 s3 s2 c4) + (quick-two-sum s1 s2 s1 c4))) + (t + (quick-two-sum s0 s1 s0 c3) + (if (/= s1 0) + (quick-two-sum s1 s2 s1 c4) + (quick-two-sum s0 s1 s0 c4)))))) + (values s0 s1 s2 s3))) (defun make-qd-d (a0 &optional (a1 0d0 a1-p) (a2 0d0) (a3 0d0)) "Create a %quad-double from four double-floats, appropriately @@ -786,17 +772,13 @@ (declare (double-float t0)) (multiple-value-bind (s1 t1) (two-sum q1 p3) - (declare (double-float t1)) + (declare (double-float s1 t1)) (multiple-value-setq (s1 t0) (two-sum s1 t0)) (incf t0 t1) - - (multiple-value-setq (s1 t0) - (quick-two-sum s1 t0)) - (multiple-value-setq (p2 t1) - (quick-two-sum s0 s1)) - (multiple-value-setq (p3 q0) - (quick-two-sum t1 t0)) + (quick-two-sum s1 t0 s1 t0) + (quick-two-sum p2 t1 s0 s1) + (quick-two-sum p3 q0 t1 t0) (let ((p4 (cl:* 2 (qd-0 a) (qd-3 a))) (p5 (cl:* 2 (qd-1 a) (qd-2 a))))

1 0

[oct-cvs] Oct commit: oct qd-methods.lisp
by rtoy 16 Sep '07

16 Sep '07

Update of /project/oct/cvsroot/oct In directory clnet:/tmp/cvs-serv10199 Modified Files: qd-methods.lisp Log Message: T should be t (for Allegro case-sensitive modern-mode). --- /project/oct/cvsroot/oct/qd-methods.lisp 2007/09/12 21:01:13 1.56 +++ /project/oct/cvsroot/oct/qd-methods.lisp 2007/09/16 02:39:29 1.57 @@ -540,7 +540,7 @@ (declare (optimize (safety 2)) (dynamic-extent more-numbers)) (do ((nlist more-numbers (cdr nlist))) - ((atom nlist) T) + ((atom nlist) t) (declare (list nlist)) (if (not (two-arg-= (car nlist) number)) (return nil)))) @@ -554,7 +554,7 @@ ((atom nlist) t) (declare (list nlist)) (unless (do* ((nl nlist (cdr nl))) - ((atom nl) T) + ((atom nl) t) (declare (list nl)) (if (two-arg-= head (car nl)) (return nil)))

1 0

[oct-cvs] Oct commit: oct qd-dd.lisp
by rtoy 13 Sep '07

13 Sep '07

Update of /project/oct/cvsroot/oct In directory clnet:/tmp/cvs-serv10722 Modified Files: qd-dd.lisp Log Message: Add some comments on how split is supposed to work. Add a note that if you have a fused multiply-subtract instruction, you can replace two-prod with a much simpler and faster version. --- /project/oct/cvsroot/oct/qd-dd.lisp 2007/09/13 16:48:48 1.7 +++ /project/oct/cvsroot/oct/qd-dd.lisp 2007/09/13 17:28:30 1.8 @@ -82,6 +82,50 @@ (values a-hi a-lo))) +;; Here is my limited understanding of what SPLIT is really supposed +;; to do. +;; +;; Let A be a double-float number. We want to split the fraction bits +;; into 2 parts of 26 bits each. This is best explained by example. +;; Let use use A = 1d50. Use INTEGER-DECODE-FLOAT to display the bits +;; of A: +;; +;; (write (integer-decode-float 1d50) :base 2) -> +;; 10001000110110000111011000101011111100110010010011010 +;; +;; Break this into 2 parts with the lower part having 27 bits and the +;; upper having 26 (because of the extra "hidden" bit): +;; +;; 10001000110110000111011000 101011111100110010010011010 +;; +;; But this is not enough. Note that the bottom half has a leading 1. +;; We want to round up the upper part. Then we need to account for +;; this in the lower part: +;; +;; 10001000110110000111011001 -10100000011001101101100110 +;; +;; This is the answer we want. Convert these back to floats with the +;; appropriate exponents, and we get: +;; +;; 1.0000000087331024d50 and -8.733102285997912d41 +;; +;; While this example worked out nicely, we should note that the +;; rounding operation above should be done in an IEEE round-to-even +;; fashion. So if the lower part of the bits is exactly "half", we +;; round the upper part to even. Thus, +;; +;; (float #b10001000110110000111011000100000000000000000000000000 1d0) +;; should be split into two parts: +;; +;; 10001000110110000111011000 100000000000000000000000000 +;; +;; but +;; +;; (float #b10001000110110000111011001100000000000000000000000000 1d0) +;; is +;; +;; 10001000110110000111011010 -100000000000000000000000000 + (defun split (a) "Split the double-float number a into a-hi and a-lo such that a = a-hi + a-lo and a-hi contains the upper 26 significant bits of a and @@ -117,6 +161,21 @@ (values a-hi a-lo)))) +;; Note that if you have an architecture that has a fused +;; multiply-subtract instruction that computes a*b-c with exactly one +;; rounding operation, you can use that instead of the complicated +;; routine below. Power PC chips have such an instruction. +;; +;; Here is the code to do that, where (fused-multiply-subtract a b p) +;; computes a*b-p. +;; +;; (defun two-prod (a b) +;; "Compute fl(a*b) and err(a*b)" +;; (declare (double-float a b)) +;; (let* ((p (* a b)) +;; (err (fused-multiply-subtract a b p))) +;; (values p err))) + (defun two-prod (a b) "Compute fl(a*b) and err(a*b)" (declare (double-float a b)

1 0

[oct-cvs] Oct commit: oct qd-dd.lisp
by rtoy 13 Sep '07

13 Sep '07

Update of /project/oct/cvsroot/oct In directory clnet:/tmp/cvs-serv4360 Modified Files: qd-dd.lisp Log Message: In SPLIT, compute the constants at compile time (#.) instead of run-time, in case the Lisp doesn't do it itself. --- /project/oct/cvsroot/oct/qd-dd.lisp 2007/09/13 16:45:42 1.6 +++ /project/oct/cvsroot/oct/qd-dd.lisp 2007/09/13 16:48:48 1.7 @@ -95,7 +95,7 @@ ;; For numbers that are very large, we use a different algorithm. ;; For smaller numbers, we can use the original algorithm of Yozo ;; Hida. - (if (> (abs a) (scale-float 1d0 (- 1023 27))) + (if (> (abs a) #.(scale-float 1d0 (- 1023 27))) ;; I've tested this algorithm against Yozo's method for 1 ;; billion randomly generated double-floats between 2^(-995) and ;; 2^996, and identical results are obtained. For numbers that @@ -103,14 +103,14 @@ ;; because of underflow. For very large numbers, we, of course ;; produce different results because Yozo's method causes ;; overflow. - (let* ((tmp (* a (+ 1 (scale-float 1d0 -27)))) - (as (* a (scale-float 1d0 -27))) - (a-hi (* (- tmp (- tmp as)) (expt 2 27))) + (let* ((tmp (* a #.(+ 1 (scale-float 1d0 -27)))) + (as (* a #.(scale-float 1d0 -27))) + (a-hi (* (- tmp (- tmp as)) #.(scale-float 1d0 27))) (a-lo (- a a-hi))) (declare (double-float tmp as a-hi a-lo)) (values a-hi a-lo)) ;; Yozo's algorithm. - (let* ((tmp (* a (+ 1 (expt 2 27)))) + (let* ((tmp (* a #.(float (+ 1 (expt 2 27)) 1d0))) (a-hi (- tmp (- tmp a))) (a-lo (- a a-hi))) (declare (double-float tmp a-hi a-lo))

1 0

[oct-cvs] Oct commit: oct qd-dd.lisp
by rtoy 13 Sep '07

13 Sep '07

Update of /project/oct/cvsroot/oct In directory clnet:/tmp/cvs-serv4233 Modified Files: qd-dd.lisp Log Message: Some updates from Richard Fateman to make these routines run faster in Allegro. --- /project/oct/cvsroot/oct/qd-dd.lisp 2007/09/13 01:06:02 1.5 +++ /project/oct/cvsroot/oct/qd-dd.lisp 2007/09/13 16:45:42 1.6 @@ -35,7 +35,7 @@ (defun quick-two-sum (a b) "Computes fl(a+b) and err(a+b), assuming |a| >= |b|" (declare (double-float a b) - (optimize (speed 3) (safety 0))) + (optimize (speed 3) (safety 0) (debug 0))) (let* ((s (+ a b)) (e (- b (- s a)))) (declare (double-float s e)) @@ -45,7 +45,7 @@ (defun two-sum (a b) "Computes fl(a+b) and err(a+b)" (declare (double-float a b) - (optimize (speed 3) (safety 0))) + (optimize (speed 3) (safety 0) (debug 0))) (let* ((s (+ a b)) (v (- s a)) (e (+ (- a (- s v)) @@ -120,7 +120,7 @@ (defun two-prod (a b) "Compute fl(a*b) and err(a*b)" (declare (double-float a b) - (optimize (speed 3) (safety 0))) + (optimize (speed 3) (safety 0) (debug 0))) (let ((p (* a b))) (declare (double-float p)) (multiple-value-bind (a-hi a-lo) @@ -141,7 +141,7 @@ "Compute fl(a*a) and err(a*b). This is a more efficient implementation of two-prod" (declare (double-float a) - (optimize (speed 3) (safety 0))) + (optimize (speed 3) (safety 0) (debug 0))) (let ((q (* a a))) (declare (double-float q)) (multiple-value-bind (a-hi a-lo)

1 0

[oct-cvs] Oct commit: oct oct.asd
by rtoy 13 Sep '07

13 Sep '07

Update of /project/oct/cvsroot/oct In directory clnet:/tmp/cvs-serv22840 Added Files: oct.asd Log Message: Initial version. --- /project/oct/cvsroot/oct/oct.asd 2007/09/13 01:07:04 NONE +++ /project/oct/cvsroot/oct/oct.asd 2007/09/13 01:07:04 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. ;;; This is the asdf definition for oct. I don't normally use this, ;;; so it might be out of date. Use at your own risk. ;; If you want all core functions to be inline (like the C++ code ;; does), add :qd-inline to *features* by enabling the following line. ;; This makes compilation much, much slower, but the resulting code ;; conses much less and is significantly faster. #+(not (and cmu x86)) (eval-when (:load-toplevel :compile-toplevel :execute) (pushnew :qd-inline *features*)) ;; To be able to inline all the functions, we need to make ;; *inline-expansion-limit* much larger. ;; ;; Not sure we really want to inline everything, but the QD C++ code ;; inlines all of the functions so we do the same. This makes CMUCL ;; take a very long time to compile the code, and the resulting ;; functions are huge. (I think div-qd is 8 KB, and sqrt-qd is a ;; whopping 30 KB!) ;; #+(and cmu qd-inline) (eval-when (:load-toplevel :compile-toplevel :execute) (setf ext:*inline-expansion-limit* 1600)) (defpackage #:oct-system (:use #:cl)) (in-package #:oct-system) (asdf:defsystem oct :description "A portable implementation of quad-double arithmetic. See <http://www.common-lisp.net/project/oct>." :author "Raymond Toy" :maintainer "See <http://www.common-lisp.net/project/oct" :licence "MIT" :version "0.0" ; No real version yet :components ((:file "qd-package") (:file "qd-rep" :depends-on ("qd-package")) #-cmu (:file "qd-dd" :depends-on ("qd-package")) (:file "qd" :depends-on ("qd-rep")) (:file "qd-io" :depends-on ("qd")) (:file "qd-const" :depends-on ("qd-io")) (:file "qd-fun" :depends-on ("qd" "qd-const")) (:file "qd-class" :depends-on ("qd-fun")) (:file "qd-methods" :depends-on ("qd-class")) (:file "qd-format" :depends-on ("qd-methods")) (:file "qd-complex" :depends-on ("qd-methods")) ))

1 0

[oct-cvs] Oct commit: oct qd-dd.lisp
by rtoy 13 Sep '07

13 Sep '07

Update of /project/oct/cvsroot/oct In directory clnet:/tmp/cvs-serv22721 Modified Files: qd-dd.lisp Log Message: Add declarations. Mostly to help Allegro generate much better code. --- /project/oct/cvsroot/oct/qd-dd.lisp 2007/08/25 17:08:48 1.4 +++ /project/oct/cvsroot/oct/qd-dd.lisp 2007/09/13 01:06:02 1.5 @@ -34,19 +34,23 @@ (declaim (inline quick-two-sum)) (defun quick-two-sum (a b) "Computes fl(a+b) and err(a+b), assuming |a| >= |b|" - (declare (double-float a b)) + (declare (double-float a b) + (optimize (speed 3) (safety 0))) (let* ((s (+ a b)) (e (- b (- s a)))) + (declare (double-float s e)) (values s e))) (declaim (inline two-sum)) (defun two-sum (a b) "Computes fl(a+b) and err(a+b)" - (declare (double-float a b)) + (declare (double-float a b) + (optimize (speed 3) (safety 0))) (let* ((s (+ a b)) (v (- s a)) (e (+ (- a (- s v)) (- b v)))) + (declare (double-float s v e)) (values s e))) (declaim (inline two-prod)) @@ -103,36 +107,46 @@ (as (* a (scale-float 1d0 -27))) (a-hi (* (- tmp (- tmp as)) (expt 2 27))) (a-lo (- a a-hi))) + (declare (double-float tmp as a-hi a-lo)) (values a-hi a-lo)) ;; Yozo's algorithm. (let* ((tmp (* a (+ 1 (expt 2 27)))) (a-hi (- tmp (- tmp a))) (a-lo (- a a-hi))) + (declare (double-float tmp a-hi a-lo)) (values a-hi a-lo)))) (defun two-prod (a b) "Compute fl(a*b) and err(a*b)" - (declare (double-float a b)) + (declare (double-float a b) + (optimize (speed 3) (safety 0))) (let ((p (* a b))) + (declare (double-float p)) (multiple-value-bind (a-hi a-lo) (split a) + (declare (double-float a-hi a-lo)) (multiple-value-bind (b-hi b-lo) (split b) + (declare (double-float b-hi b-lo)) (let ((e (+ (+ (- (* a-hi b-hi) p) (* a-hi b-lo) (* a-lo b-hi)) (* a-lo b-lo)))) + (declare (double-float e)) (values p e)))))) (declaim (inline two-sqr)) (defun two-sqr (a) "Compute fl(a*a) and err(a*b). This is a more efficient implementation of two-prod" - (declare (double-float a)) + (declare (double-float a) + (optimize (speed 3) (safety 0))) (let ((q (* a a))) + (declare (double-float q)) (multiple-value-bind (a-hi a-lo) (split a) + (declare (double-float a-hi a-lo)) (values q (+ (+ (- (* a-hi a-hi) q) (* 2 a-hi a-lo)) (* a-lo a-lo))))))

1 0

[oct-cvs] Oct commit: oct qd-complex.lisp qd-methods.lisp qd-package.lisp
by rtoy 12 Sep '07

12 Sep '07

Update of /project/oct/cvsroot/oct In directory clnet:/tmp/cvs-serv9448 Modified Files: qd-complex.lisp qd-methods.lisp qd-package.lisp Log Message: qd-package.lisp: o Rearrange some exports so the CMU ones are all grouped together. o Export new constants pi/2, pi/4, 2pi, and log2. o Export the qd-real and qd-complex types. qd-methods.lisp: o Define new constants for pi/2, pi/4, 2pi, and log2. o Update some of the macrolets to work with a modern-mode lisp, like Allegro. qd-complex.lisp: o Use the new constants as needed. --- /project/oct/cvsroot/oct/qd-complex.lisp 2007/08/31 21:13:36 1.35 +++ /project/oct/cvsroot/oct/qd-complex.lisp 2007/09/12 21:01:13 1.36 @@ -350,7 +350,7 @@ (let ((t0 (/ 1 (sqrt #q2.0q0))) (t1 #q1.2q0) (t2 #q3q0) - (ln2 #.(log #q2.0)) + (ln2 +log2+) (x (realpart z)) (y (imagpart z))) (multiple-value-bind (rho k) @@ -407,7 +407,7 @@ (let* ( ;; Constants (theta (/ (sqrt most-positive-double-float) 4.0d0)) (rho (/ 4.0d0 (sqrt most-positive-double-float))) - (half-pi #.(/ +pi+ 2d0)) + (half-pi +pi/2+) (rp (realpart z)) (beta (float-sign rp)) (x (* beta rp)) --- /project/oct/cvsroot/oct/qd-methods.lisp 2007/08/31 21:13:36 1.55 +++ /project/oct/cvsroot/oct/qd-methods.lisp 2007/09/12 21:01:13 1.56 @@ -26,7 +26,24 @@ (in-package #:qd) (defconstant +pi+ - (make-instance 'qd-real :value qdi:+qd-pi+)) + (make-instance 'qd-real :value qdi:+qd-pi+) + "Quad-double value of pi") + +(defconstant +pi/2+ + (make-instance 'qd-real :value qdi:+qd-pi/2+) + "Quad-double value of pi/2") + +(defconstant +pi/4+ + (make-instance 'qd-real :value qdi:+qd-pi/4+) + "Quad-double value of pi/4") + +(defconstant +2pi+ + (make-instance 'qd-real :value qdi:+qd-2pi+) + "Quad-double value of 2*pi") + +(defconstant +log2+ + (make-instance 'qd-real :value qdi:+qd-log2+) + "Quad-double value of log(2), natural log of 2") #+cmu (defconstant +quad-double-float-positive-infinity+ @@ -200,9 +217,13 @@ (unary-divide number))) (macrolet ((frob (name &optional (type 'real)) - (let ((method-name (intern (concatenate 'string "Q" (symbol-name name)))) + (let ((method-name (intern (concatenate 'string + (string '#:q) + (symbol-name name)))) (cl-name (intern (symbol-name name) :cl)) - (qd-name (intern (concatenate 'string (symbol-name name) "-QD")))) + (qd-name (intern (concatenate 'string + (symbol-name name) + (string '#:-qd))))) `(progn (defmethod ,method-name ((x ,type)) (,cl-name x)) @@ -318,9 +339,11 @@ (macrolet ((frob (op) - (let ((method (intern (concatenate 'string "TWO-ARG-" (symbol-name op)))) + (let ((method (intern (concatenate 'string + (string '#:two-arg-) + (symbol-name op)))) (cl-fun (find-symbol (symbol-name op) :cl)) - (qd-fun (intern (concatenate 'string "QD-" (symbol-name op)) + (qd-fun (intern (concatenate 'string (string '#:qd-) (symbol-name op)) (find-package :qdi)))) `(progn (defmethod ,method ((a real) (b real)) @@ -352,9 +375,11 @@ (macrolet ((frob (name) (let ((method-name - (intern (concatenate 'string "Q" (symbol-name name)))) + (intern (concatenate 'string (string '#:q) + (symbol-name name)))) (cl-name (intern (symbol-name name) :cl)) - (qd-name (intern (concatenate 'string (symbol-name name) "-QD")))) + (qd-name (intern (concatenate 'string (symbol-name name) + (string '#:-qd))))) `(progn (defmethod ,name ((x number)) (,cl-name x)) @@ -828,7 +853,9 @@ ;; the corresponding two-arg-<foo> function. (macrolet ((frob (op) - (let ((method (intern (concatenate 'string "TWO-ARG-" (symbol-name op))))) + (let ((method (intern (concatenate 'string + (string '#:two-arg-) + (symbol-name op))))) `(define-compiler-macro ,op (number &rest more-numbers) (do* ((n number (car nlist)) (nlist more-numbers (cdr nlist)) --- /project/oct/cvsroot/oct/qd-package.lisp 2007/08/29 01:22:03 1.36 +++ /project/oct/cvsroot/oct/qd-package.lisp 2007/09/12 21:01:13 1.37 @@ -30,11 +30,9 @@ #:read-qd #:add-qd #:add-qd-d - #:cmu #:add-qd-dd #:add-d-qd #:sub-qd #:sub-qd-d - #:cmu #:sub-qd-dd #:sub-d-qd #:neg-qd #:mul-qd @@ -42,9 +40,7 @@ #:sqr-qd #:div-qd #:div-qd-d - #+cmu #:div-qd-dd #:make-qd-d - #+cmu #:make-qd-dd #:integer-decode-qd #:npow #:qd-0 @@ -55,6 +51,10 @@ #:+qd-one+ #:+qd-zero+ #:+qd-pi+ + #:+qd-pi/2+ + #:+qd-pi/4+ + #:+qd-2pi+ + #:+qd-log2+ ;; Functions #:hypot-qd #:abs-qd @@ -91,6 +91,11 @@ #:random-qd ) #+cmu + (:export #:add-qd-dd + #:sub-qd-dd + #:div-qd-dd + #:make-qd-dd) + #+cmu (:import-from #:c #:two-sum #:quick-two-sum @@ -164,6 +169,10 @@ #:decf #:float-digits ) + ;; Export types + (:export #:qd-real + #:qd-complex) + ;; Export functions (:export #:+ #:- #:* @@ -229,7 +238,11 @@ #:float-digits ) ;; Constants - (:export #:+pi+) + (:export #:+pi+ + #:+pi/2+ + #:+pi/4+ + #:+2pi+ + #:+log2+) ;; CMUCL supports infinities. #+cmu (:export #:+quad-double-float-positive-infinity+

1 0

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

12 Sep '07

Update of /project/oct/cvsroot/oct In directory clnet:/tmp/cvs-serv18006 Modified Files: qd-fun.lisp qd-io.lisp qd.lisp Log Message: qd-fun.lisp: o Remove unused var R1 in EXP-QD/REDUCE. o TAN-QD was calling ZEROP instead of ZEROP-QD. o Comment out extra copy of ASINH-QD. qd-io.lisp: o Ignore unused var in QD-PRINT-EXPONENT and QD-READER. qd.lisp: o Remove extra version of DIV-QD. --- /project/oct/cvsroot/oct/qd-fun.lisp 2007/08/25 17:08:48 1.79 +++ /project/oct/cvsroot/oct/qd-fun.lisp 2007/09/12 02:31:14 1.80 @@ -160,8 +160,6 @@ (let* ((k 256) (z (truncate (qd-0 (nint-qd (div-qd a +qd-log2+))))) - (r1 (sub-qd a (mul-qd-d +qd-log2+ (float z 1d0)))) - ;; r as above (r (div-qd-d (sub-qd a (mul-qd-d +qd-log2+ (float z 1d0))) (float k 1d0))) ;; For Taylor series. p = r^2/2, the first term @@ -740,7 +738,7 @@ (defun tan-qd (r) "Tan(r)" (declare (type %quad-double r)) - (if (zerop r) + (if (zerop-qd r) r (tan-qd/sincos r))) @@ -810,6 +808,7 @@ (d (expm1-qd a2))) (div-qd d (add-qd-d d 2d0)))))) +#+(or) (defun asinh-qd (a) "Asinh(a)" (declare (type %quad-double a)) --- /project/oct/cvsroot/oct/qd-io.lisp 2007/08/27 17:49:19 1.14 +++ /project/oct/cvsroot/oct/qd-io.lisp 2007/09/12 02:31:14 1.15 @@ -151,6 +151,7 @@ (scale r s m+ m-))))))) (defun qd-print-exponent (x exp stream) + (declare (ignore x)) (let ((*print-radix* nil)) (format stream "q~D" exp))) @@ -461,6 +462,7 @@ (make-float sign int-part frac-part scale exp))))) (defun qd-reader (stream subchar arg) + (declare (ignore subchar arg)) (read-qd stream)) (defun qd-from-string (string) --- /project/oct/cvsroot/oct/qd.lisp 2007/08/25 17:08:48 1.45 +++ /project/oct/cvsroot/oct/qd.lisp 2007/09/12 02:31:14 1.46 @@ -819,7 +819,6 @@ (renorm-5 p0 p1 p2 p3 p4)))))))))) -#-cmu (defun div-qd (a b) (declare (type %quad-double a b) (optimize (speed 3) @@ -840,25 +839,6 @@ (renorm-4 q0 q1 q2 q3) (%make-qd-d q0 q1 q2 q3))))))) -(defun div-qd (a b) - (declare (type %quad-double a b) - (optimize (speed 3) - (space 0))) - (let ((a0 (qd-0 a)) - (b0 (qd-0 b))) - (let* ((q0 (cl:/ a0 b0)) - (r (sub-qd a (mul-qd-d b q0))) - (q1 (cl:/ (qd-0 r) b0))) - (when (float-infinity-p q0) - (return-from div-qd (%make-qd-d q0 0d0 0d0 0d0))) - (setf r (sub-qd r (mul-qd-d b q1))) - (let ((q2 (cl:/ (qd-0 r) b0))) - (setf r (sub-qd r (mul-qd-d b q2))) - (let ((q3 (cl:/ (qd-0 r) b0))) - (multiple-value-bind (q0 q1 q2 q3) - (renorm-4 q0 q1 q2 q3) - (%make-qd-d q0 q1 q2 q3))))))) - ;; Non-sloppy divide #+(or) (defun div-qd (a b)

1 0