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

[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

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

12 Sep '07

Update of /project/oct/cvsroot/oct In directory clnet:/tmp/cvs-serv13071 Modified Files: .cvsignore Log Message: Ignore fasl. --- /project/oct/cvsroot/oct/.cvsignore 2007/08/27 13:10:50 1.2 +++ /project/oct/cvsroot/oct/.cvsignore 2007/09/12 02:03:42 1.3 @@ -2,3 +2,4 @@ *.sparcf *.x86f *.err +*.fasl

1 0

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

06 Sep '07

Update of /project/oct/cvsroot/oct In directory clnet:/tmp/cvs-serv23948 Modified Files: qd-class.lisp Log Message: Fix qd-complex printer so it prints out something that is readable. --- /project/oct/cvsroot/oct/qd-class.lisp 2007/08/31 21:13:36 1.24 +++ /project/oct/cvsroot/oct/qd-class.lisp 2007/09/06 02:58:38 1.25 @@ -67,7 +67,7 @@ (make-instance 'qd-real :value (qd-value x))) (defmethod print-object ((qd qd-complex) stream) - (format stream "#q(~<~/qdi::qd-format/ ~/qdi::qd-format/~:@>)" + (format stream "#q(~<#q~/qdi::qd-format/ #q~/qdi::qd-format/~:@>)" (list (qd-real qd) (qd-imag qd))))

1 0