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

[oct-cvs] Oct commit: oct qd-methods.lisp
by rtoy 31 Jul '08

31 Jul '08

Update of /project/oct/cvsroot/oct In directory clnet:/tmp/cvs-serv30205 Modified Files: qd-methods.lisp Log Message: For CMUCL, define compiler macros to convert two-arg-foo into the appropriate CL function or QD-REAL function so we don't have to do CLOS dispatch, if the types are known at compile-time. --- /project/oct/cvsroot/oct/qd-methods.lisp 2008/07/18 17:02:04 1.66 +++ /project/oct/cvsroot/oct/qd-methods.lisp 2008/07/31 19:13:42 1.67 @@ -990,6 +990,76 @@ (if (cdr more-numbers) form `(not (two-arg-= ,number ,(car more-numbers))))) + + +;; Define compiler macro the convert two-arg-foo into the appropriate +;; CL function or QD-REAL function so we don't have to do CLOS +;; dispatch. +#+cmu +(macrolet + ((frob (name cl-op qd-op) + `(define-compiler-macro ,name (&whole form x y &environment env) + (flet ((arg-type (arg) + (multiple-value-bind (def-type localp decl) + (ext:variable-information arg env) + (declare (ignore localp)) + (when def-type + (cdr (assoc 'type decl)))))) + (let ((x-type (arg-type x)) + (y-type (arg-type y))) + (cond ((and (subtypep x-type 'cl:number) + (subtypep y-type 'cl:number)) + `(,',cl-op ,x ,y)) + ((and (subtypep x-type 'qd-real) + (subtypep y-type 'qd-real)) + `(make-instance 'qd-real :value (,',qd-op (qd-value ,x) + (qd-value ,y)))) + (t + ;; Don't know how to handle this, so give up. + form))))))) + (frob two-arg-+ cl:+ add-qd) + (frob two-arg-- cl:- sub-qd) + (frob two-arg-* cl:* mul-qd) + (frob two-arg-/ cl:/ div-qd)) + +#+cmu +(macrolet + ((frob (name cl-op qd-op cl-qd-op qd-cl-op) + `(define-compiler-macro ,name (&whole form x y &environment env) + (flet ((arg-type (arg) + (multiple-value-bind (def-type localp decl) + (ext:variable-information arg env) + (declare (ignore localp)) + (when def-type + (cdr (assoc 'type decl)))))) + (let ((x-type (arg-type x)) + (y-type (arg-type y))) + (cond ((subtypep x-type 'cl:float) + (cond ((subtypep y-type 'cl:number) + `(,',cl-op ,x ,y)) + ((subtypep y-type 'qd-real) + (if ,cl-qd-op + `(make-instance 'qd-real :value (,',cl-qd-op (cl:float ,x 1d0) + (qd-value ,y))) + form)) + (t form))) + ((subtypep x-type 'qd-real) + (cond ((subtypep y-type 'cl:float) + (if ,qd-cl-op + `(make-instance 'qd-real :value (,',qd-cl-op (qd-value ,x) + (float ,y 1d0))) + form)) + ((subtypep y-type 'qd-real) + `(make-instance 'qd-real :value (,',qd-op (qd-value ,x) + (qd-value ,y)))) + (t form))) + (t + ;; Don't know how to handle this, so give up. + form))))))) + (frob two-arg-+ cl:+ add-qd add-d-qd add-qd-d) + (frob two-arg-- cl:- sub-qd sub-d-qd sub-qd-d) + (frob two-arg-* cl:* mul-qd mul-d-qd mul-qd-d) + (frob two-arg-/ cl:/ div-qd nil nil)) (defun read-qd-real-or-complex (stream)

1 0

[oct-cvs] Oct commit: oct qd-methods.lisp
by rtoy 18 Jul '08

18 Jul '08

Update of /project/oct/cvsroot/oct In directory clnet:/tmp/cvs-serv16709 Modified Files: qd-methods.lisp Log Message: Add or cleanup some docstrings. --- /project/oct/cvsroot/oct/qd-methods.lisp 2007/10/16 13:46:01 1.65 +++ /project/oct/cvsroot/oct/qd-methods.lisp 2008/07/18 17:02:04 1.66 @@ -27,54 +27,61 @@ (defconstant +pi+ (make-instance 'qd-real :value octi:+qd-pi+) - "Quad-double value of pi") + "Pi represented as a QD-REAL") (defconstant +pi/2+ (make-instance 'qd-real :value octi:+qd-pi/2+) - "Quad-double value of pi/2") + "Pi/2 represented as a QD-REAL") (defconstant +pi/4+ (make-instance 'qd-real :value octi:+qd-pi/4+) - "Quad-double value of pi/4") + "Pi/4 represented as a QD-REAL") (defconstant +2pi+ (make-instance 'qd-real :value octi:+qd-2pi+) - "Quad-double value of 2*pi") + "2*pi represented as a QD-REAL") (defconstant +log2+ (make-instance 'qd-real :value octi:+qd-log2+) - "Quad-double value of log(2), natural log of 2") + "Natural log of 2 represented as a QD-REAL") +;; How do we represent infinity for a QD-REAL? For now, we just make +;; the QD-REAL whose most significant part is infinity. Currently +;; only supported on CMUCL. #+cmu (defconstant +quad-double-float-positive-infinity+ (make-instance 'qd-real :value (make-qd-d ext:double-float-positive-infinity)) - "Positive infinity for qd-real") + "One representation of positive infinity for QD-REAL") #+cmu (defconstant +quad-double-float-negative-infinity+ (make-instance 'qd-real :value (make-qd-d ext:double-float-negative-infinity)) - "Negative infinity for qd-real") + "One representation of negative infinity for QD-REAL") (defconstant +most-positive-quad-double-float+ (make-instance 'qd-real :value (octi::%make-qd-d most-positive-double-float (cl:scale-float most-positive-double-float (cl:* 1 -53)) (cl:scale-float most-positive-double-float (cl:* 2 -53)) - (cl:scale-float most-positive-double-float (cl:* 3 -53))))) + (cl:scale-float most-positive-double-float (cl:* 3 -53)))) + "Most positive representable QD-REAL") (defconstant +least-positive-quad-double-float+ (make-instance 'qd-real - :value (make-qd-d least-positive-double-float))) + :value (make-qd-d least-positive-double-float)) + "Least positive QD-REAL") ;; Not sure this is 100% correct, but I think if the first component ;; is any smaller than this, the last component would no longer be a ;; normalized double-float. (defconstant +least-positive-normalized-quad-double-float+ (make-instance 'qd-real - :value (make-qd-d (cl:scale-float least-positive-normalized-double-float (cl:* 3 53))))) + :value (make-qd-d (cl:scale-float least-positive-normalized-double-float (cl:* 3 53)))) + "Least positive normalized QD-REAL") (defconstant +qd-real-one+ - (make-instance 'qd-real :value (make-qd-d 1d0))) + (make-instance 'qd-real :value (make-qd-d 1d0)) + "QD-REAL representation of 1") (defmethod make-qd ((x cl:rational))

1 0

[oct-cvs] Oct commit: oct qd-class.lisp
by rtoy 18 Jul '08

18 Jul '08

Update of /project/oct/cvsroot/oct In directory clnet:/tmp/cvs-serv16645 Modified Files: qd-class.lisp Log Message: Add or cleanup some docstrings. --- /project/oct/cvsroot/oct/qd-class.lisp 2007/10/16 02:39:21 1.29 +++ /project/oct/cvsroot/oct/qd-class.lisp 2008/07/18 17:01:51 1.30 @@ -33,7 +33,8 @@ ((qd :initform +qd-zero+ :reader qd-value :initarg :value - :type %quad-double))) + :type %quad-double)) + (:documentation "QUAD-DOUBLE real number")) (defclass qd-complex () ((real :initform +qd-zero+ @@ -43,7 +44,8 @@ (imag :initform +qd-zero+ :reader qd-imag :initarg :imag - :type %quad-double))) + :type %quad-double)) + (:documentation "Complex number consisting of QUAD-DOUBLE components")) #-cmu (defmethod print-object ((qd qd-real) stream)

1 0

[oct-cvs] Oct commit: oct qd-rep.lisp qd.lisp
by rtoy 17 Jul '08

17 Jul '08

Update of /project/oct/cvsroot/oct In directory clnet:/tmp/cvs-serv29306 Modified Files: qd-rep.lisp qd.lisp Log Message: Add or cleanup more docstrings. --- /project/oct/cvsroot/oct/qd-rep.lisp 2008/07/16 21:02:07 1.14 +++ /project/oct/cvsroot/oct/qd-rep.lisp 2008/07/17 17:26:43 1.15 @@ -78,7 +78,7 @@ (eval-when (:compile-toplevel :load-toplevel :execute) (defun %make-qd-d (a0 a1 a2 a3) - "Make a %quad-double from 4 double-floats, exactly using the given + "Make a %QUAD-DOUBLE from four double-floats, exactly using the given values. No check is made to see if the values make sense. A0 is the most significant part and A3 is the least. " @@ -243,6 +243,12 @@ ) ; end progn +;; Define some compiler macros to transform add-qd to add-qd-t +;; directly. For CMU without :oct-array, we always replace the +;; parameter C with NIL because we don't use it. For other Lisps, we +;; create the necessary object and call add-qd-t. +;; +;; Do the same mul-qd and other similar functions. (macrolet ((frob (qd qd-t) #-oct-array --- /project/oct/cvsroot/oct/qd.lisp 2008/07/16 21:02:07 1.66 +++ /project/oct/cvsroot/oct/qd.lisp 2008/07/17 17:26:43 1.67 @@ -290,8 +290,10 @@ (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 - normalizing the result from the four double-floats. + "Create a %QUAD-DOUBLE from four double-floats, appropriately +normalizing the result from the four double-floats. A0 is the most +significant part of the %QUAD-DOUBLE, and A3 is the least. The optional +parameters default to 0. " (declare (double-float a0 a1 a2 a3) (optimize (speed 3) @@ -307,7 +309,7 @@ ;; Quad-double + double (defun add-qd-d (a b &optional (target #+oct-array (%make-qd-d 0d0 0d0 0d0 0d0))) - "Add the %QUAD-DOUBLE A and the DOUBLE-FLOAT B, returning a %QUAD-DOUBLE. + "Return the sum of the %QUAD-DOUBLE A and the DOUBLE-FLOAT B. If TARGET is given, TARGET is destructively modified to contain the result." (add-qd-d-t a b target)) @@ -339,7 +341,8 @@ (%store-qd-d target r0 r1 r2 r3))))) (defun add-d-qd (a b &optional (target #+oct-array (%make-qd-d 0d0 0d0 0d0 0d0))) - "Add the DOUBLE-FLOAT A and the %QUAD-DOUBLE B, returning a %QUAD-DOUBLE" + "Return the sum of the DOUBLE-FLOAT A and the %QUAD-DOUBLE B. +If TARGET is given, TARGET is destructively modified to contain the result." (declare (double-float a) (type %quad-double b) (optimize (speed 3)) @@ -407,6 +410,8 @@ ;; those compilers. (defun add-qd (a b &optional (target #+oct-array (%make-qd-d 0d0 0d0 0d0 0d0))) + "Return the sum of the %QUAD-DOUBLE numbers A and B. +If TARGET is given, TARGET is destructively modified to contain the result." (add-qd-t a b target)) @@ -471,11 +476,9 @@ (%store-qd-d target (+ a0 b0) 0d0 0d0 0d0) (%store-qd-d target s0 s1 s2 s3))))))))))) -;; Define some compiler macros to transform add-qd to add-qd-t -;; directly. For CMU, we always replace the parameter C with NIL -;; because we don't use it. For other Lisps, we create the necessary -;; object and call add-qd-t. (defun neg-qd (a &optional (target #+oct-array (%make-qd-d 0d0 0d0 0d0 0d0))) + "Return the negative of the %QUAD-DOUBLE number A. +If TARGET is given, TARGET is destructively modified to contain the result." (neg-qd-t a target)) (defun neg-qd-t (a target) @@ -487,6 +490,8 @@ (%store-qd-d target (cl:- a0) (cl:- a1) (cl:- a2) (cl:- a3)))) (defun sub-qd (a b &optional (target #+oct-array (%make-qd-d 0d0 0d0 0d0 0d0))) + "Return the difference between the %QUAD-DOUBLE numbers A and B. +If TARGET is given, TARGET is destructively modified to contain the result." (declare (type %quad-double a b)) (add-qd-t a (neg-qd b) target)) @@ -500,12 +505,20 @@ (add-qd-dd a (cl:- b))) (defun sub-qd-d (a b &optional (target #+oct-array (%make-qd-d 0d0 0d0 0d0 0d0))) + "Return the difference between the %QUAD-DOUBLE number A and +the DOUBLE-FLOAT number B. + +If TARGET is given, TARGET is destructively modified to contain the result." (declare (type %quad-double a) (type double-float b) #+(and cmu (not oct-array)) (ignore target)) (add-qd-d a (cl:- b) #+oct-array target)) (defun sub-d-qd (a b &optional (target #+oct-array (%make-qd-d 0d0 0d0 0d0 0d0))) + "Return the difference between the DOUBLE-FLOAT number A and +the %QUAD-DOUBLE number B. + +If TARGET is given, TARGET is destructively modified to contain the result." (declare (type double-float a) (type %quad-double b) #+(and cmu (not oct-array)) (ignore target)) @@ -521,6 +534,10 @@ ;; 14.142135623730950488016887242096980785696718753769480731766797379908L0 ;; (defun mul-qd-d (a b &optional (target #+oct-array (%make-qd-d 0d0 0d0 0d0 0d0))) + "Return the product of the %QUAD-DOUBLE number A and +the DOUBLE-FLOAT number B. + +If TARGET is given, TARGET is destructively modified to contain the result." (mul-qd-d-t a b target)) (defun mul-qd-d-t (a b target) @@ -648,6 +665,8 @@ ;; 14.142135623730950488016887242096980785696718753769480731766797379908L0 (defun mul-qd (a b &optional (target #+oct-array (%make-qd-d 0d0 0d0 0d0 0d0))) + "Returns the product of the %QUAD-DOUBLE numbers A and B. +If TARGET is given, TARGET is destructively modified to contain the result." (mul-qd-t a b target)) (defun mul-qd-t (a b target) @@ -975,6 +994,7 @@ ;; quad-double / double (defun div-qd-d (a b) + "Divides the %QUAD-DOUBLE number A by the DOUBLE-FLOAT number B." (declare (type %quad-double a) (double-float b) (optimize (speed 3) @@ -1030,7 +1050,8 @@ #+cmu (defun make-qd-dd (a0 a1) - "Create a %quad-double from two double-double-floats" + "Create a %QUAD-DOUBLE from two double-double-floats. This is for CMUCL, +which supports double-double floats." (declare (double-double-float a0 a1) (optimize (speed 3) (space 0))) (make-qd-d (kernel:double-double-hi a0) @@ -1043,7 +1064,7 @@ (declaim (ext:end-block)) (defun abs-qd (a) - "Absolute value of the %QUAD-DOUBLE A" + "Returns the absolute value of the %QUAD-DOUBLE A" (declare (type %quad-double a)) (if (minusp (float-sign (qd-0 a))) (neg-qd a) @@ -1051,6 +1072,8 @@ ;; a^n for an integer n (defun npow (a n) + "Return N'th power of A, where A is a %QUAD-DOUBLE number and N +is a fixnum." (declare (type %quad-double a) (fixnum n) (optimize (speed 3) @@ -1108,7 +1131,7 @@ (div-qd (make-qd-d 1d0) r)))) (defun qd-< (a b) - "A < B" + "Returns T if A < B, where A and B are %QUAD-DOUBLE numbers." (or (< (qd-0 a) (qd-0 b)) (and (= (qd-0 a) (qd-0 b)) (or (< (qd-1 a) (qd-1 b)) @@ -1118,7 +1141,7 @@ (< (qd-3 a) (qd-3 b))))))))) (defun qd-> (a b) - "A > B" + "Returns T if A > B, where A and B are %QUAD-DOUBLE numbers." (or (> (qd-0 a) (qd-0 b)) (and (= (qd-0 a) (qd-0 b)) (or (> (qd-1 a) (qd-1 b)) @@ -1128,20 +1151,20 @@ (> (qd-3 a) (qd-3 b))))))))) (defun qd-<= (a b) - "A > B" + "Returns T if A <= B, where A and B are %QUAD-DOUBLE numbers." (not (qd-> a b))) (defun qd->= (a b) - "A > B" + "Returns T if A >= B, where A and B are %QUAD-DOUBLE numbers." (not (qd-< a b))) (defun zerop-qd (a) - "Is A zero?" + "Returns T if the %QUAD-DOUBLE number A is numerically equal to 0." (declare (type %quad-double a)) (zerop (qd-0 a))) (defun onep-qd (a) - "Is A equal to 1?" + "Returns T if the %QUAD-DOUBLE number A is numerically equal to 1." (declare (type %quad-double a)) (and (= (qd-0 a) 1d0) (zerop (qd-1 a)) @@ -1149,16 +1172,17 @@ (zerop (qd-3 a)))) (defun plusp-qd (a) - "Is A positive?" + "Returns T if the %QUAD-DOUBLE number A is strictly positive." (declare (type %quad-double a)) (plusp (qd-0 a))) (defun minusp-qd (a) - "Is A negative?" + "Returns T if the %QUAD-DOUBLE number A is strictly negative." (declare (type %quad-double a)) (minusp (qd-0 a))) (defun qd-= (a b) + "Returns T if the %QUAD-DOUBLE numbers A and B are numerically equal." (and (= (qd-0 a) (qd-0 b)) (= (qd-1 a) (qd-1 b)) (= (qd-2 a) (qd-2 b)) @@ -1239,6 +1263,13 @@ q0-sign))))))))) (defun integer-decode-qd (q) + "Like INTEGER-DECODE-FLOAT, but for %QUAD-DOUBLE numbers. + Returns three values: + 1) an integer representation of the significand. + 2) the exponent for the power of 2 that the significand must be multiplied + by to get the actual value. + 3) -1 or 1 (i.e. the sign of the argument.)" + (declare (type %quad-double q)) ;; Integer decode each component and then create the appropriate ;; integer by shifting and adding all the parts together. If any @@ -1321,6 +1352,11 @@ (defun decode-float-qd (q) + "Like DECODE-FLOAT, but for %QUAD-DOUBLE numbers. Returns three values: + 1) a %QUAD-DOUBLE number representing the significand. This is always + between 0.5 (inclusive) and 1.0 (exclusive). + 2) an integer representing the exponent. + 3) -1.0 or 1.0 (i.e. the sign of the argument.)" (declare (type %quad-double q)) (multiple-value-bind (frac exp sign) (decode-float (qd-0 q))

1 0

[oct-cvs] Oct commit: oct oct.asd
by rtoy 17 Jul '08

17 Jul '08

Update of /project/oct/cvsroot/oct In directory clnet:/tmp/cvs-serv29243 Modified Files: oct.asd Log Message: Fix minor typo. --- /project/oct/cvsroot/oct/oct.asd 2007/09/16 05:00:00 1.2 +++ /project/oct/cvsroot/oct/oct.asd 2008/07/17 17:26:21 1.3 @@ -55,7 +55,7 @@ (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" + :maintainer "See <http://www.common-lisp.net/project/oct>" :licence "MIT" :version "0.0" ; No real version yet :components

1 0

[oct-cvs] Oct commit: oct qd-const.lisp qd-fun.lisp qd-io.lisp qd-rep.lisp qd.lisp
by rtoy 16 Jul '08

16 Jul '08

Update of /project/oct/cvsroot/oct In directory clnet:/tmp/cvs-serv21392 Modified Files: qd-const.lisp qd-fun.lisp qd-io.lisp qd-rep.lisp qd.lisp Log Message: Add or cleanup some docstrings. --- /project/oct/cvsroot/oct/qd-const.lisp 2007/10/16 13:44:00 1.20 +++ /project/oct/cvsroot/oct/qd-const.lisp 2008/07/16 21:02:07 1.21 @@ -26,10 +26,12 @@ (in-package #:octi) (defconstant +qd-zero+ - (make-qd-d 0d0)) + (make-qd-d 0d0) + "%QUAD-DOUBLE representation of 0") (defconstant +qd-one+ - (make-qd-d 1d0)) + (make-qd-d 1d0) + "%QUAD-DOUBLE representation of 1") ;; The bits of 2/pi. Scale these bits by 2^(-1584) and you'll get ;; 2/pi. These are used for accurate argument reduction for the trig @@ -46,7 +48,8 @@ (scale-float (float -8753721960665020 1.0d0) -161) (scale-float (float 5857755168774013 1.0d0) -215) (scale-float (float 5380502254059520 1.0d0) -269)) - (%make-qd-d q0 q1 q2 q3))) + (%make-qd-d q0 q1 q2 q3)) + "%QUAD-DOUBLE representation of pi") ;; 6.2831853071795864769252867665590057683943387987502116419498891846156328125724L0 ;; #q6.2831853071795864769252867665590057683943387987502116419498891846q0 @@ -57,7 +60,8 @@ (scale-float (float -8753721960665020 1.0d0) -160) (scale-float (float 5857755168774013 1.0d0) -214) (scale-float (float 5380502254059520 1.0d0) -268)) - (%make-qd-d q0 q1 q2 q3))) + (%make-qd-d q0 q1 q2 q3)) + "%QUAD-DOUBLE representation of 2*pi") ;; 1.5707963267948966192313216916397514420985846996875529104874722961539082031431L0 ;; #q1.57079632679489661923132169163975144209858469968755291048747229615q0 @@ -68,7 +72,8 @@ (scale-float (float -8753721960665020 1.0d0) -162) (scale-float (float 5857755168774013 1.0d0) -216) (scale-float (float 5380502254059520 1.0d0) -270)) - (%make-qd-d q0 q1 q2 q3))) + (%make-qd-d q0 q1 q2 q3)) + "%QUAD-DOUBLE representation of pi/2") ;; 0.78539816339744830961566084581987572104929234984377645524373614807695410157155L0 ;; #q0.785398163397448309615660845819875721049292349843776455243736148076q0 @@ -79,7 +84,8 @@ (scale-float (float -8753721960665020 1.0d0) -163) (scale-float (float 5857755168774013 1.0d0) -217) (scale-float (float 5380502254059520 1.0d0) -271)) - (%make-qd-d q0 q1 q2 q3))) + (%make-qd-d q0 q1 q2 q3)) + "%QUAD-DOUBLE representation of pi/4") ;; 2.35619449019234492884698253745962716314787704953132936573120844423086230471467L0 ;; #q2.35619449019234492884698253745962716314787704953132936573120844423q0 @@ -90,7 +96,8 @@ (scale-float (float 5724553519491610 1.0d0) -160) (scale-float (float -6810541066450737 1.0d0) -214) (scale-float (float -7491566988951552 1.0d0) -273)) - (%make-qd-d q0 q1 q2 q3))) + (%make-qd-d q0 q1 q2 q3)) + "%QUAD-DOUBLE representation of 3*pi/4") ;; 0.00306796157577128245943617517898388953534879824157725177829584432842560195926387L0 ;; #q0.00306796157577128245943617517898388953534879824157725177829584432842q0 @@ -125,7 +132,8 @@ (%make-qd-d (scale-float (float 6243314768165359 1.0d0) -53) (scale-float (float 7525737178955839 1.0d0) -108) (scale-float (float 6673460182522164 1.0d0) -163) - (scale-float (float -7545482916914641 1.0d0) -217))) + (scale-float (float -7545482916914641 1.0d0) -217)) + "%QUAD-DOUBLE representation of log(2) (natural log)") ;; The rest of log(2) such that (+ +qd-log2+ +qd-log2-extra+) is ;; log(2) to twice the precision of a quad-double. --- /project/oct/cvsroot/oct/qd-fun.lisp 2007/11/28 20:00:28 1.92 +++ /project/oct/cvsroot/oct/qd-fun.lisp 2008/07/16 21:02:07 1.93 @@ -81,7 +81,7 @@ (scale-float-qd (mul-qd r new-a) (ash k -1))))) (defun sqrt-qd (a) - "Square root of the (non-negative) quad-float" + "Return the square root of the (non-negative) %QUAD-DOUBLE number A" (declare (type %quad-double a) (optimize (speed 3) (space 0))) ;; Perform the following Newton iteration: @@ -130,7 +130,8 @@ (- k)))) (defun hypot-qd (x y) - "sqrt(x^2+y^2) computed carefully without unnecessary overflow" + "sqrt(x^2+y^2) computed carefully without unnecessary overflow for +the %QUAD-DOUBLE numbers X and Y" (multiple-value-bind (abs^2 rho) (hypot-aux-qd x y) (scale-float-qd (sqrt-qd abs^2) rho))) @@ -168,7 +169,7 @@ (make-qd-d s0 s1 s2 s3)))) (defun ffloor-qd (a) - "The floor of A, returned as a quad-float" + "The floor of the %QUAD-DOUBLE A, returned as a %QUAD-DOUBLE number" (let ((x0 (ffloor (qd-0 a))) (x1 0d0) (x2 0d0) @@ -254,7 +255,8 @@ (mul-qd d (add-qd-d d 2d0)))))) (defun expm1-qd (a) - "exp(a) - 1, done accurately" + "Return exp(a) - 1 for the %QUAD-DOUBLE number A. This is more + accurate than just computing exp(a) - 1 directly." (declare (type %quad-double a)) (when (float-infinity-p (qd-0 a)) @@ -265,7 +267,7 @@ (expm1-qd/duplication a)) (defun exp-qd (a) - "exp(a)" + "Return the expnential of the %QUAD-DOUBLE number A" (declare (type %quad-double a)) ;; Should we try to be more accurate than just 709? (when (< (qd-0 a) (log least-positive-normalized-double-float)) @@ -348,8 +350,8 @@ (mul-qd-d sum 2d0))))) (defun log1p-qd (x) - "log1p(x) = log(1+x), done more accurately than just evaluating - log(1+x)" + "Return log1p(x) = log(1+x), done more accurately than just evaluating + log(1+x). X is a non-negative %QUAD-DOUBLE number" (declare (type %quad-double x)) (when (float-infinity-p (qd-0 x)) @@ -357,7 +359,7 @@ (log1p-qd/duplication x)) (defun log-qd (a) - "Log(a)" + "Return the (natural) log of the non-negative %QUAD-DOUBLE number A" (declare (type %quad-double a)) (cond ((onep-qd a) +qd-zero+) @@ -499,7 +501,7 @@ (neg-qd s))))))))))) (defun cos-qd (a) - "Cos(a)" + "Return the cosine of the %QUAD-DOUBLE number A" ;; Just like sin-qd, but for cos. (declare (type %quad-double a)) ;; To compute sin(x), choose integers a, b so that @@ -586,6 +588,8 @@ ;; Compute sin and cos of a (defun sincos-qd (a) + "Return the sine of the %QUAD-DOUBLE number A. The second returned value +is the cosine of A" (declare (type %quad-double a)) (when (zerop-qd a) (return-from sincos-qd @@ -784,7 +788,7 @@ (defun sin-qd (a) - "Sin(a)" + "Return the sine of the %QUAD-DOUBLE number A" (declare (type %quad-double a)) ;; To compute sin(x), choose integers a, b so that ;; @@ -861,7 +865,7 @@ (neg-qd c))))))))) (defun cos-qd (a) - "Cos(a)" + "Return the cosine of the %QUAD-DOUBLE number A" ;; Just like sin-qd, but for cos. (declare (type %quad-double a)) ;; To compute sin(x), choose integers a, b so that @@ -1100,18 +1104,18 @@ (atan2-qd/newton y x)) (defun atan-qd (y) - "Atan4b*(y)" + "Return the arc tangent of the %QUAD-DOUBLE number Y" (declare (type %quad-double y)) (atan-qd/newton y)) (defun asin-qd (a) - "Asin(a)" + "Return the arc sine of the %QUAD-DOUBLE number A" (declare (type %quad-double a)) (atan2-qd a (sqrt-qd (sub-d-qd 1d0 (sqr-qd a))))) (defun acos-qd (a) - "Acos(a)" + "Return the arc cosine of the %QUAD-DOUBLE number A" (declare (type %quad-double a)) (atan2-qd (sqrt-qd (sub-d-qd 1d0 (sqr-qd a))) @@ -1128,7 +1132,7 @@ (div-qd s c))) (defun tan-qd (r) - "Tan(r)" + "Return the tangent of the %QUAD-DOUBLE number A" (declare (type %quad-double r)) (if (zerop-qd r) r @@ -1136,7 +1140,7 @@ (defun sinh-qd (a) - "Sinh(a)" + "Return the hyperbolic sine of the %QUAD-DOUBLE number A" (declare (type %quad-double a)) ;; Hart et al. suggests sinh(x) = 1/2*(D(x) + D(x)/(D(x)+1)) ;; where D(x) = exp(x) - 1. This helps for x near 0. @@ -1153,7 +1157,7 @@ -1))))) (defun cosh-qd (a) - "Cosh(a)" + "Return the hyperbolic cosine of the %QUAD-DOUBLE number A" (declare (type %quad-double a)) ;; cosh(x) = 1/2*(exp(x)+exp(-x)) (let ((e (exp-qd a))) @@ -1163,7 +1167,7 @@ -1))) (defun tanh-qd (a) - "Tanh(a)" + "Return the hyperbolic tangent of the %QUAD-DOUBLE number A" (declare (type %quad-double a)) ;; Hart et al. suggests tanh(x) = D(2*x)/(2+D(2*x)) (cond ((zerop (qd-0 a)) @@ -1231,7 +1235,7 @@ (log1p-qd (sqrt-qd (add-qd-d (sqr-qd 1/a) 1d0)))))))) (defun asinh-qd (a) - "Asinh(a)" + "Return the inverse hyperbolic sine of the %QUAD-DOUBLE number A" (declare (type %quad-double a)) ;; asinh(x) = log(x + sqrt(1+x^2)) ;; @@ -1266,7 +1270,7 @@ (log1p-qd (sqrt-qd (add-qd-d (sqr-qd 1/a) 1d0))))))))) (defun acosh-qd (a) - "Acosh(a)" + "Return the inverse hyperbolic cosine of the %QUAD-DOUBLE number A" (declare (type %quad-double a)) ;; acosh(x) = log(x + sqrt(x^2-1)) #+nil @@ -1302,7 +1306,7 @@ (sqrt-qd (add-d-qd 1d0 1/a))))))))) (defun atanh-qd (a) - "Atanh(a)" + "Return the inverse hyperbolic tangent of the %QUAD-DOUBLE number A" (declare (type %quad-double a)) ;; atanh(x) = 1/2*log((1+x)/(1-x)) ;; = 1/2*log(1+(2*x)/(1-x)) @@ -1322,7 +1326,7 @@ (defun random-qd (&optional (state *random-state*)) - "Generate a quad-double random number in the range [0,1)" + "Generate a %QUAD-DOUBLE random number in the range [0,1)" (declare (optimize (speed 3))) ;; Strategy: Generate 31 bits at a time, shift the bits and repeat 7 times. (let* ((r +qd-zero+) --- /project/oct/cvsroot/oct/qd-io.lisp 2007/10/16 17:05:13 1.21 +++ /project/oct/cvsroot/oct/qd-io.lisp 2008/07/16 21:02:07 1.22 @@ -385,6 +385,7 @@ ;; convert bignums to qd. This supports converting rationals to qd ;; too. (defun rational-to-qd (rat) + "Convert a rational number RAT to a %QUAD-DOUBLE number" (declare (rational rat)) (let* ((p (coerce rat 'double-float)) (ans (make-qd-d p)) @@ -408,6 +409,9 @@ ;; This seems to work, but really needs to be rewritten! (defun read-qd (stream) + "Read a %QUAD-DOUBLE number from the stream STREAM. The format of the number +should be like a float, but with extra significant digits allowed. An exponent +marker of Q is allowed." (labels ((read-digits (s) ;; Read a sequence of digits and return the decimal ;; value, the character that terminated the sequence, and --- /project/oct/cvsroot/oct/qd-rep.lisp 2007/11/28 20:00:28 1.13 +++ /project/oct/cvsroot/oct/qd-rep.lisp 2008/07/16 21:02:07 1.14 @@ -55,18 +55,23 @@ ;; quad-double. QD-0 is the most significant part and QD-3 is the ;; least. (defun qd-0 (q) + "Return the most significant double-float in the %QUAD-DOUBLE number Q" (declare (type %quad-double q) (optimize (speed 3))) (kernel:double-double-hi (realpart q))) (defun qd-1 (q) + "Return the second most significant double-float in the %QUAD-DOUBLE number Q" (declare (type %quad-double q) (optimize (speed 3))) (kernel:double-double-lo (realpart q))) (defun qd-2 (q) + "Return the third most significant double-float in the %QUAD-DOUBLE number Q" (declare (type %quad-double q) (optimize (speed 3))) (kernel:double-double-hi (imagpart q))) (defun qd-3 (q) + "Return the fourth most significant (least significant) double-float in the + %QUAD-DOUBLE number Q" (declare (type %quad-double q) (optimize (speed 3))) (kernel:double-double-lo (imagpart q))) --- /project/oct/cvsroot/oct/qd.lisp 2008/02/11 17:04:13 1.65 +++ /project/oct/cvsroot/oct/qd.lisp 2008/07/16 21:02:07 1.66 @@ -307,6 +307,8 @@ ;; Quad-double + double (defun add-qd-d (a b &optional (target #+oct-array (%make-qd-d 0d0 0d0 0d0 0d0))) + "Add the %QUAD-DOUBLE A and the DOUBLE-FLOAT B, returning a %QUAD-DOUBLE. +If TARGET is given, TARGET is destructively modified to contain the result." (add-qd-d-t a b target)) (defun add-qd-d-t (a b target) @@ -337,6 +339,7 @@ (%store-qd-d target r0 r1 r2 r3))))) (defun add-d-qd (a b &optional (target #+oct-array (%make-qd-d 0d0 0d0 0d0 0d0))) + "Add the DOUBLE-FLOAT A and the %QUAD-DOUBLE B, returning a %QUAD-DOUBLE" (declare (double-float a) (type %quad-double b) (optimize (speed 3)) @@ -815,6 +818,8 @@ (renorm-5 p0 p1 s0 t0 t1)))))))))))))))))))) (defun sqr-qd (a &optional (target #+oct-array (%make-qd-d 0d0 0d0 0d0 0d0))) + "Return the square of the %QUAD-DOUBLE number A. If TARGET is also given, +it is destructively modified with the result." (sqr-qd-t a target)) (defun sqr-qd-t (a target) @@ -867,6 +872,8 @@ (defun div-qd (a b &optional (target #+oct-array (%make-qd-d 0d0 0d0 0d0 0d0))) + "Return the quotient of the two %QUAD-DOUBLE numbers A and B. +If TARGET is given, it destrutively modified with the result." (div-qd-t a b target)) #+nil @@ -1036,6 +1043,7 @@ (declaim (ext:end-block)) (defun abs-qd (a) + "Absolute value of the %QUAD-DOUBLE A" (declare (type %quad-double a)) (if (minusp (float-sign (qd-0 a))) (neg-qd a) @@ -1288,6 +1296,7 @@ (cl:* (qd-3 qd) scale)))) (defun scale-float-qd (qd k) + "Scale the %QUAD-DOUBLE number QD by 2^K. Just like SCALE-FLOAT" (declare (type %quad-double qd) (type (integer -2000 2000) k) (optimize (speed 3) (space 0)))

1 0