oct-cvs

oct-cvs@common-lisp.net

202 discussions

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

05 Nov '07

Update of /project/oct/cvsroot/oct In directory clnet:/tmp/cvs-serv7545 Modified Files: Tag: THREE-ARG-BRANCH timing2.lisp Log Message: Update timings. --- /project/oct/cvsroot/oct/timing2.lisp 2007/10/16 14:21:13 1.4 +++ /project/oct/cvsroot/oct/timing2.lisp 2007/11/05 16:00:53 1.4.2.1 @@ -117,140 +117,143 @@ Test Time qd oct ---- ----------- -add 0.023 0.09 -mul 0.075 0.13 -div 0.299 0.29 -sqrt 0.105 0.11 -sin 0.115 0.14 -log 0.194 0.12 - -Times are in sec for the test. The default number of iterations were -used. Most of the timings match my expectations, including the log -test. Oct uses a different algorithm (Halley's method) which is -faster (in Lisp) than the algorithm used in qd (Newtwon iteration). +add 0.23 1.16 +mul 0.749 1.54 +div 3.00 3.11 +sqrt 10.57 12.2 +sin 57.33 64.5 +log 194 119 + +Times are in microsec/operation for the test. The default number of +iterations were used. Most of the timings match my expectations, +including the log test. Oct uses a different algorithm (Halley's +method) which is faster (in Lisp) than the algorithm used in qd +(Newtwon iteration). +Not also that these times include the 3-arg versions of the routines. ------------------------------------------------------------------------------- The raw data: The output from qd_timer -qd -v: +Timing qd_real +-------------- + Timing addition... -n = 100000 t = 0.0231462 -b = 1.428571e+04 -100000 operations in 0.0231462 s. - 0.231462 us +n = 1000000 t = 0.236154 +b = 142857.142857 +1000000 operations in 0.236154 s. + 0.236154 us Timing multiplication ... -n = 100000 t = 0.0749929 -b = 2.718268e+00 -100000 operations in 0.0749929 s. - 0.749929 us +n = 1000000 t = 0.748933 +b = 2.718280 +1000000 operations in 0.748933 s. + 0.748933 us Timing division ... -n = 100000 t = 0.298858 -b = 0.367881 -100000 operations in 0.298858 s. - 2.988580 us +n = 1000000 t = 3.004328 +b = 0.367880 +1000000 operations in 3.004328 s. + 3.004328 us Timing square root ... -n = 10000 t = 0.105049 +n = 100000 t = 1.057170 a = 2.821980 -10000 operations in 0.105049 s. - 10.504860 us +100000 operations in 1.057170 s. + 10.571696 us Timing sin ... -n = 2000 t = 0.114943 +n = 20000 t = 1.146667 a = 3.141593 -2000 operations in 0.114943 s. - 57.471350 us +20000 operations in 1.146667 s. + 57.333335 us Timing log ... -n = 1000 t = 0.193698 +n = 10000 t = 1.939869 a = -50.100000 -1000 operations in 0.193698 s. -193.697800 us -The output from CMUCL: +10000 operations in 1.939869 s. +193.986900 us + +-------------------------------------------------- +CMUCL results: -QD> (time-add) +CL-USER> (oct::time-add 1000000) ; 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 +; 1.16 seconds of real time +; 0.98 seconds of user run time +; 0.18 seconds of system run time +; 1,845,637,904 CPU cycles ; 0 page faults and -; 7,200,016 bytes consed. +; 72,000,248 bytes consed. ; -n = 100000 -b = #q14285.7142857142857142857142857142857142857142857142857142857142855q0 -NIL -QD> (time-mul) +n = 1000000 +b = #q142857.142857142857142857142857142857142857142857142857142857142854q0 + +CL-USER> (oct::time-mul 1000000) ; 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 +; 1.53 seconds of real time +; 1.27 seconds of user run time +; 0.25 seconds of system run time +; 2,430,859,732 CPU cycles ; 0 page faults and -; 7,200,824 bytes consed. +; 72,000,248 bytes consed. ; -n = 100000 -b = #q2.71826823717448966803506482442604644797444693267782286300915989397q0 -NIL -QD> (time-div) +n = 1000000 +b = #q2.71828046931937688381979970845435639275164502668250771294016782123q0 + +CL-USER> (oct::time-div 1000000) ; 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 +; 3.11 seconds of real time +; 2.94 seconds of user run time +; 0.16 seconds of system run time +; 4,957,512,968 CPU cycles ; 0 page faults and -; 7,200,016 bytes consed. +; 72,000,248 bytes consed. ; -n = 100000 -b = #q0.36788128056098406210328658773118942247132502490133718973918140856q0 -NIL -QD> (time-sqrt 10000) +n = 1000000 +b = #q0.367879625111086265804761271038216553876450599098470428879260437304q0 +CL-USER> (oct::time-sqrt 100000) ; 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 +; 1.22 seconds of real time +; 1.1 seconds of user run time +; 0.1 seconds of system run time +; 1,938,798,996 CPU cycles ; 0 page faults and -; 2,402,560 bytes consed. +; 24,000,128 bytes consed. ; -n = 10000 +n = 100000 a = #q2.82198033014704783016853125515542796898998765943212617578596649019q0 -NIL -QD> (time-sin) + +CL-USER> (oct::time-sin 20000) ; 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 +; 1.29 seconds of real time +; 1.24 seconds of user run time +; 0.05 seconds of system run time +; 2,053,157,408 CPU cycles ; 0 page faults and -; 3,105,800 bytes consed. +; 27,751,144 bytes consed. ; -n = 2000 -a = #q3.14159265358979323846264338327950288419716939937510582097494459409q0 -NIL -QD> (time-log) +n = 20000 +a = #q3.14159265358979323846264338327950288419716939937510582097494458294q0 + +CL-USER> (oct::time-log 10000) ; 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 +; 1.19 seconds of real time +; 1.13 seconds of user run time +; 0.04 seconds of system run time +; 1,890,677,952 CPU cycles ; 0 page faults and -; 1,621,792 bytes consed. +; 16,197,536 bytes consed. ; -n = 1000 -a = #q-50.100000000000000000000000000000000000000000000000000000000208796q0 -NIL -QD> +n = 10000 +a = #q-50.100000000000000000000000000000000000000000000000000000552824575q0 ---------------------------------------------- ||#

1 0

[oct-cvs] Oct commit: oct qd.lisp
by rtoy 05 Nov '07

05 Nov '07

Update of /project/oct/cvsroot/oct In directory clnet:/tmp/cvs-serv7505 Modified Files: Tag: THREE-ARG-BRANCH qd.lisp Log Message: Cosmetic --- /project/oct/cvsroot/oct/qd.lisp 2007/11/04 16:31:45 1.60.2.4 +++ /project/oct/cvsroot/oct/qd.lisp 2007/11/05 16:00:39 1.60.2.5 @@ -870,7 +870,8 @@ (declaim (inline invert-qd)) -(defun invert-qd(v) ;; a quartic newton iteration for 1/v +(defun invert-qd (v) + ;; a quartic newton iteration for 1/v ;; to invert v, start with a good guess, x. ;; let h= 1-v*x ;; h is small ;; return x+ x*(h+h^2+h^3) . compute h3 in double-float

1 0

[oct-cvs] Oct commit: oct qd.lisp
by rtoy 04 Nov '07

04 Nov '07

Update of /project/oct/cvsroot/oct In directory clnet:/tmp/cvs-serv18778 Modified Files: Tag: THREE-ARG-BRANCH qd.lisp Log Message: Oops. Make sqr-qd-t inline. --- /project/oct/cvsroot/oct/qd.lisp 2007/11/04 03:00:56 1.60.2.3 +++ /project/oct/cvsroot/oct/qd.lisp 2007/11/04 16:31:45 1.60.2.4 @@ -167,6 +167,7 @@ mul-qd mul-qd-t sqr-qd + sqr-qd-t div-qd div-qd-t div-qd-d

1 0

[oct-cvs] Oct commit: oct qd-rep.lisp qd.lisp
by rtoy 04 Nov '07

04 Nov '07

Update of /project/oct/cvsroot/oct In directory clnet:/tmp/cvs-serv10089 Modified Files: Tag: THREE-ARG-BRANCH qd-rep.lisp qd.lisp Log Message: o Add support for SQR-QD-T o Add compiler macro for SQR-QD. --- /project/oct/cvsroot/oct/qd-rep.lisp 2007/11/04 02:45:01 1.10.2.3 +++ /project/oct/cvsroot/oct/qd-rep.lisp 2007/11/04 03:00:56 1.10.2.4 @@ -258,3 +258,14 @@ #-cmu (define-compiler-macro sub-qd (a b &optional c) `(add-qd-t ,a (neg-qd ,b) ,c)) + +#+cmu +(define-compiler-macro sqr-qd (a &optional c) + (if c + `(setf ,c (sqr-qd-t ,a nil)) + `(sqr-qd-t ,a nil))) + +#-cmu +(define-compiler-macro sqr-qd (a b &optional c) + `(sqr-qd-t ,a ,c)) + --- /project/oct/cvsroot/oct/qd.lisp 2007/11/04 02:45:01 1.60.2.2 +++ /project/oct/cvsroot/oct/qd.lisp 2007/11/04 03:00:56 1.60.2.3 @@ -626,7 +626,9 @@ (declare (type %quad-double a b) (optimize (speed 3) (space 0)) - (inline float-infinity-p)) + (inline float-infinity-p) + #+cmu + (ignore target)) (with-qd-parts (a0 a1 a2 a3) a (declare (double-float a0 a1 a2 a3)) @@ -787,11 +789,16 @@ (multiple-value-call #'%make-qd-d (renorm-5 p0 p1 s0 t0 t1)))))))))))))))))))) -(defun sqr-qd (a) +(defun sqr-qd (a &optional (target #-cmu (%make-qd-d 0d0 0d0 0d0 0d0))) + (sqr-qd-t a target)) + +(defun sqr-qd-t (a target) "Square A" (declare (type %quad-double a) (optimize (speed 3) - (space 0))) + (space 0)) + #+cmu + (ignore target)) (multiple-value-bind (p0 q0) (two-sqr (qd-0 a)) (multiple-value-bind (p1 q1) @@ -831,7 +838,7 @@ (multiple-value-bind (a0 a1 a2 a3) (renorm-5 p0 p1 p2 p3 p4) - (%make-qd-d a0 a1 a2 a3))))))))) + (%store-qd-d target a0 a1 a2 a3))))))))) (defun div-qd (a b &optional (target #-cmu (%make-qd-d 0d0 0d0 0d0 0d0)))

1 0

[oct-cvs] Oct commit: oct qd-extra.lisp
by rtoy 04 Nov '07

04 Nov '07

Update of /project/oct/cvsroot/oct In directory clnet:/tmp/cvs-serv8162 Modified Files: Tag: THREE-ARG-BRANCH qd-extra.lisp Log Message: Update timing data. --- /project/oct/cvsroot/oct/qd-extra.lisp 2007/10/26 15:48:15 1.5 +++ /project/oct/cvsroot/oct/qd-extra.lisp 2007/11/04 02:51:39 1.5.2.1 @@ -1099,14 +1099,14 @@ ;; Time Sparc PPC x86 PPC (fma) Sparc2 ;; exp-qd/reduce 2.06 3.18 10.46 2.76 6.12 ;; expm1-qd/series 8.81 12.24 18.87 3.26 29.0 -;; expm1-qd/dup 5.68 4.34 18.47 3.64 18.78 -;; exp-qd/pade 1.53 +;; expm1-qd/dup 5.68 4.34 18.47 3.64 9.77 +;; exp-qd/pade 1.53 4.51 ;; ;; Consing (MB) Sparc ;; exp-qd/reduce 45 45 638 44.4 45 ;; expm1-qd/series 519 519 1201 14.8 519 ;; expm1-qd/dup 32 32 1224 32.0 32 -;; exp-qd/pade 44 +;; exp-qd/pade 44 44 ;; ;; Speeds seem to vary quite a bit between architectures. ;;

1 0

[oct-cvs] Oct commit: oct qd-rep.lisp qd.lisp
by rtoy 04 Nov '07

04 Nov '07

Update of /project/oct/cvsroot/oct In directory clnet:/tmp/cvs-serv6695 Modified Files: Tag: THREE-ARG-BRANCH qd-rep.lisp qd.lisp Log Message: o Move compiler macros from qd.lisp to qd-rep.lisp o Declare add-qd-t, mul-qd-t and div-qd-t as inline functions so that everything is still fast on cmucl. --- /project/oct/cvsroot/oct/qd-rep.lisp 2007/11/02 20:45:32 1.10.2.2 +++ /project/oct/cvsroot/oct/qd-rep.lisp 2007/11/04 02:45:01 1.10.2.3 @@ -233,3 +233,28 @@ (declare (ignore x)) nil) ) ; end progn + + +(macrolet + ((frob (qd qd-t) + #+cmu + `(define-compiler-macro ,qd (a b &optional c) + (if c + `(setf ,c (,',qd-t ,a ,b nil)) + `(,',qd-t ,a ,b nil))) + #-cmu + `(define-compiler-macro ,qd (a b &optional (c (%make-qd-d 0d0 0d0 0d0 0d0))) + `(,',qd-t ,a ,b ,c)))) + (frob add-qd add-qd-t) + (frob mul-qd mul-qd-t) + (frob div-qd div-qd-t)) + +#+cmu +(define-compiler-macro sub-qd (a b &optional c) + (if c + `(setf ,c (add-qd-t ,a (neg-qd ,b) nil)) + `(add-qd-t ,a (neg-qd ,b) nil))) + +#-cmu +(define-compiler-macro sub-qd (a b &optional c) + `(add-qd-t ,a (neg-qd ,b) ,c)) --- /project/oct/cvsroot/oct/qd.lisp 2007/11/02 20:11:42 1.60.2.1 +++ /project/oct/cvsroot/oct/qd.lisp 2007/11/04 02:45:01 1.60.2.2 @@ -161,11 +161,14 @@ add-qd-d add-qd-dd add-qd + add-qd-t mul-qd-d mul-qd-dd mul-qd + mul-qd-t sqr-qd div-qd + div-qd-t div-qd-d div-qd-dd)) @@ -187,6 +190,7 @@ mul-qd-t sqr-qd div-qd div-qd-d div-qd-dd + div-qd-t make-qd-dd )) @@ -394,7 +398,9 @@ (defun add-qd-t (a b target) (declare (type %quad-double a b) (optimize (speed 3) - (space 0))) + (space 0)) + #+cmu + (ignore target)) ;; This is the version that is NOT IEEE. Should we use the IEEE ;; version? It's quite a bit more complicated. ;; @@ -454,16 +460,6 @@ ;; 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. -#+cmu -(define-compiler-macro add-qd (a b &optional c) - (if c - `(setf c (add-qd-t ,a ,b nil)) - `(add-qd-t ,a ,b nil))) - -#-cmu -(define-compiler-macro add-qd (a b &optional (c (%make-qd-d 0d0 0d0 0d0 0d0))) - `(add-qd-t ,a ,b ,c)) - (defun neg-qd (a) (declare (type %quad-double a)) (with-qd-parts (a0 a1 a2 a3) @@ -476,16 +472,6 @@ (add-qd-t a (neg-qd b) target)) #+cmu -(define-compiler-macro sub-qd (a b &optional c) - (if c - `(setf ,c `(add-qd-t ,a (neg-qd ,b) nil)) - `(add-qd-t ,a (neg-qd ,b) nil))) - -#-cmu -(define-compiler-macro sub-qd (a b &optional (c (%make-qd-d 0d0 0d0 0d0 0d0))) - `(add-qd-t ,a (neg-qd ,b) ,c)) - -#+cmu (defun sub-qd-dd (a b) (declare (type %quad-double a) (type double-double-float b)) @@ -699,16 +685,6 @@ (%store-qd-d target p0 0d0 0d0 0d0) (%store-qd-d target r0 r1 s0 s1)))))))))))))) -#+cmu -(define-compiler-macro mul-qd (a b &optional c) - (if c - `(setf ,c `(mul-qd-t ,a ,b nil)) - `(mul-qd-t ,a ,b nil))) - -#-cmu -(define-compiler-macro mul-qd (a b &optional (c (%make-qd-d 0d0 0d0 0d0 0d0))) - `(mul-qd-t ,a ,b ,c)) - ;; This is the non-sloppy version. I think this works just fine, but ;; since qd defaults to the sloppy multiplication version, we do the @@ -865,7 +841,9 @@ (declare (type %quad-double a b) (optimize (speed 3) (space 0)) - (inline float-infinity-p)) + (inline float-infinity-p) + #+cmu + (ignore target)) (let ((a0 (qd-0 a)) (b0 (qd-0 b))) (let* ((q0 (cl:/ a0 b0)) @@ -882,17 +860,6 @@ (renorm-4 q0 q1 q2 q3) (%store-qd-d target q0 q1 q2 q3))))))) -#+cmu -(define-compiler-macro div-qd (a b &optional c) - (if c - `(setf ,c `(div-qd-t ,a ,b nil)) - `(div-qd-t ,a ,b nil))) - -#-cmu -(define-compiler-macro div-qd (a b &optional (c (%make-qd-d 0d0 0d0 0d0 0d0))) - `(div-qd-t ,a ,b ,c)) - - (declaim (inline invert-qd)) (defun invert-qd(v) ;; a quartic newton iteration for 1/v

1 0

[oct-cvs] Oct commit: oct qd-rep.lisp
by rtoy 02 Nov '07

02 Nov '07

Update of /project/oct/cvsroot/oct In directory clnet:/tmp/cvs-serv14009 Modified Files: Tag: THREE-ARG-BRANCH qd-rep.lisp Log Message: Oops. %STORE-QD-D wasn't returning the target value. --- /project/oct/cvsroot/oct/qd-rep.lisp 2007/11/02 20:11:42 1.10.2.1 +++ /project/oct/cvsroot/oct/qd-rep.lisp 2007/11/02 20:45:32 1.10.2.2 @@ -178,7 +178,8 @@ (setf (aref ,dest 0) ,q0) (setf (aref ,dest 1) ,q1) (setf (aref ,dest 2) ,q2) - (setf (aref ,dest 3) ,q3)))) + (setf (aref ,dest 3) ,q3) + ,dest))) (defun qd-parts (qd) "Extract the four doubles comprising a quad-double and return them

1 0

[oct-cvs] Oct commit: oct qd-rep.lisp qd.lisp
by rtoy 02 Nov '07

02 Nov '07

Update of /project/oct/cvsroot/oct In directory clnet:/tmp/cvs-serv8022 Modified Files: Tag: THREE-ARG-BRANCH qd-rep.lisp qd.lisp Log Message: First cut at adding a 3-arg versions of the basic operations to reduce consing by allowing the third argument to be a place where the result can be stored. This is intended to help reduce allocation and gc costs for Lisps that use arrays to represent quad-doubles. More work is needed to make the compiler macros do the right thing for CMUCL. qd-rep.lisp: o Add %STORE-QD-D to store a quad-double into a place. For CMUCL, there place argument is ignored and a fresh quad-double is created. qd.lisp: o Modify ADD-QD, SUB-QD, MUL-QD, and DIV-QD to take an optional third argument indicating where the result can be stored. Ignored on CMUCL. o Add ADD-QD-T, SUB-QD-T, MUL-QD-T, and DIV-QD-T, which are 3-arg functions with the third arg always required which is the storage area to hold the result. Ignored on CMUCL. o Add compiler macros to convert ADD-QD and friends to ADD-QD-T if the third arg is always given. The effect is, essentially, inlining ADD-QD. --- /project/oct/cvsroot/oct/qd-rep.lisp 2007/10/16 17:09:46 1.10 +++ /project/oct/cvsroot/oct/qd-rep.lisp 2007/11/02 20:11:42 1.10.2.1 @@ -81,6 +81,9 @@ (kernel:%make-double-double-float a2 a3))) ) +(defmacro %store-qd-d (target q0 q1 q2 q3) + (declare (ignore target)) + `(%make-qd-d ,q0 ,q1, q2, q3)) (defun qd-parts (qd) "Extract the four doubles comprising a quad-double and return them @@ -169,6 +172,14 @@ (setf (aref ,a 3) ,a3) ,a))) +(defmacro %store-qd-d (target q0 q1 q2 q3) + (let ((dest (gensym "TARGET-"))) + `(let ((,dest ,target)) + (setf (aref ,dest 0) ,q0) + (setf (aref ,dest 1) ,q1) + (setf (aref ,dest 2) ,q2) + (setf (aref ,dest 3) ,q3)))) + (defun qd-parts (qd) "Extract the four doubles comprising a quad-double and return them as multiple values. The most significant double is the first value." --- /project/oct/cvsroot/oct/qd.lisp 2007/10/18 14:38:11 1.60 +++ /project/oct/cvsroot/oct/qd.lisp 2007/11/02 20:11:42 1.60.2.1 @@ -179,10 +179,12 @@ make-qd-d add-qd-d add-d-qd add-qd-dd add-dd-qd - add-qd + add-qd add-qd-t neg-qd sub-qd sub-qd-dd sub-qd-d sub-d-qd - mul-qd-d mul-qd-dd mul-qd + mul-qd-d mul-qd-dd + mul-qd + mul-qd-t sqr-qd div-qd div-qd-d div-qd-dd make-qd-dd @@ -385,7 +387,11 @@ ;; which don't do a very good job with dataflow. CMUCL is one of ;; those compilers. -(defun add-qd (a b) +(defun add-qd (a b &optional (target #-cmu (%make-qd-d 0d0 0d0 0d0 0d0))) + (add-qd-t a b target)) + + +(defun add-qd-t (a b target) (declare (type %quad-double a b) (optimize (speed 3) (space 0))) @@ -407,7 +413,7 @@ (inline float-infinity-p)) (when (float-infinity-p s0) - (return-from add-qd (%make-qd-d s0 0d0 0d0 0d0))) + (return-from add-qd-t (%store-qd-d target s0 0d0 0d0 0d0))) (let ((v0 (cl:- s0 a0)) (v1 (cl:- s1 a1)) (v2 (cl:- s2 a2)) @@ -441,8 +447,22 @@ (multiple-value-setq (s0 s1 s2 s3) (renorm-5 s0 s1 s2 s3 t0)) (if (and (zerop a0) (zerop b0)) - (%make-qd-d (+ a0 b0) 0d0 0d0 0d0) - (%make-qd-d s0 s1 s2 s3))))))))))) + (%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. +#+cmu +(define-compiler-macro add-qd (a b &optional c) + (if c + `(setf c (add-qd-t ,a ,b nil)) + `(add-qd-t ,a ,b nil))) + +#-cmu +(define-compiler-macro add-qd (a b &optional (c (%make-qd-d 0d0 0d0 0d0 0d0))) + `(add-qd-t ,a ,b ,c)) (defun neg-qd (a) (declare (type %quad-double a)) @@ -451,9 +471,19 @@ (declare (double-float a0 a1 a2 a3)) (%make-qd-d (cl:- a0) (cl:- a1) (cl:- a2) (cl:- a3)))) -(defun sub-qd (a b) +(defun sub-qd (a b &optional (target #-cmu (%make-qd-d 0d0 0d0 0d0 0d0))) (declare (type %quad-double a b)) - (add-qd a (neg-qd b))) + (add-qd-t a (neg-qd b) target)) + +#+cmu +(define-compiler-macro sub-qd (a b &optional c) + (if c + `(setf ,c `(add-qd-t ,a (neg-qd ,b) nil)) + `(add-qd-t ,a (neg-qd ,b) nil))) + +#-cmu +(define-compiler-macro sub-qd (a b &optional (c (%make-qd-d 0d0 0d0 0d0 0d0))) + `(add-qd-t ,a (neg-qd ,b) ,c)) #+cmu (defun sub-qd-dd (a b) @@ -602,7 +632,11 @@ ;; ;; Clisp says ;; 14.142135623730950488016887242096980785696718753769480731766797379908L0 -(defun mul-qd (a b) + +(defun mul-qd (a b &optional (target #-cmu (%make-qd-d 0d0 0d0 0d0 0d0))) + (mul-qd-t a b target)) + +(defun mul-qd-t (a b target) (declare (type %quad-double a b) (optimize (speed 3) (space 0)) @@ -617,7 +651,7 @@ (two-prod a0 b0) #+cmu (when (float-infinity-p p0) - (return-from mul-qd (%make-qd-d p0 0d0 0d0 0d0))) + (return-from mul-qd-t (%store-qd-d target p0 0d0 0d0 0d0))) (multiple-value-bind (p1 q1) (two-prod a0 b1) (multiple-value-bind (p2 q2) @@ -662,8 +696,19 @@ (multiple-value-bind (r0 r1 s0 s1) (renorm-5 p0 p1 s0 s1 s2) (if (zerop r0) - (%make-qd-d p0 0d0 0d0 0d0) - (%make-qd-d r0 r1 s0 s1)))))))))))))) + (%store-qd-d target p0 0d0 0d0 0d0) + (%store-qd-d target r0 r1 s0 s1)))))))))))))) + +#+cmu +(define-compiler-macro mul-qd (a b &optional c) + (if c + `(setf ,c `(mul-qd-t ,a ,b nil)) + `(mul-qd-t ,a ,b nil))) + +#-cmu +(define-compiler-macro mul-qd (a b &optional (c (%make-qd-d 0d0 0d0 0d0 0d0))) + `(mul-qd-t ,a ,b ,c)) + ;; This is the non-sloppy version. I think this works just fine, but ;; since qd defaults to the sloppy multiplication version, we do the @@ -813,7 +858,10 @@ (%make-qd-d a0 a1 a2 a3))))))))) -(defun div-qd (a b) +(defun div-qd (a b &optional (target #-cmu (%make-qd-d 0d0 0d0 0d0 0d0))) + (div-qd-t a b target)) + +(defun div-qd-t (a b target) (declare (type %quad-double a b) (optimize (speed 3) (space 0)) @@ -825,14 +873,25 @@ (q1 (cl:/ (qd-0 r) b0))) (when (float-infinity-p q0) - (return-from div-qd (%make-qd-d q0 0d0 0d0 0d0))) + (return-from div-qd-t (%store-qd-d target 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))))))) + (%store-qd-d target q0 q1 q2 q3))))))) + +#+cmu +(define-compiler-macro div-qd (a b &optional c) + (if c + `(setf ,c `(div-qd-t ,a ,b nil)) + `(div-qd-t ,a ,b nil))) + +#-cmu +(define-compiler-macro div-qd (a b &optional (c (%make-qd-d 0d0 0d0 0d0 0d0))) + `(div-qd-t ,a ,b ,c)) + (declaim (inline invert-qd))

1 0

[oct-cvs] Oct commit: oct qd-extra.lisp
by rtoy 26 Oct '07

26 Oct '07

Update of /project/oct/cvsroot/oct In directory clnet:/tmp/cvs-serv1229 Modified Files: qd-extra.lisp Log Message: Add Pade approximation for exp. From Richard Fateman. Experiments show that it is 25% faster than exp-qd/reduce. However, it loses 3 bits of accuracy compared with exp-qd/reduce. --- /project/oct/cvsroot/oct/qd-extra.lisp 2007/10/15 18:54:02 1.4 +++ /project/oct/cvsroot/oct/qd-extra.lisp 2007/10/26 15:48:15 1.5 @@ -1018,6 +1018,71 @@ (mul-qd +cordic-scale+ y)))) +;; Evaluate the polynomial poly at the point q. The polynomial is +;; a list of the coefficients arranged in descending powers. +(defun poly-eval-qd (q poly) + (let ((sum (%make-qd-d 0d0 0d0 0d0 0d0))) + (dolist (c poly) + (setf sum (add-qd (mul-qd q sum) c))) + sum)) + +;; Like poly-eval-qd, except the polynomial is a list of double-floats +(defun poly-eval-qd-d (q poly) + (let ((sum (%make-qd-d 0d0 0d0 0d0 0d0))) + (dolist (c poly) + (setf sum (add-qd-d (mul-qd q sum) c))) + sum)) + + +;; This idea is from Richard Fateman. I've added a few additional +;; notes, but these are mostly from Richard. +;; +;; We can evaluate exp(x) for |x|< 0.01, say, by a pade approximation, +;; computed using Maxima: +;; +;; taylor(exp(x), x, 0, 20)$ +;; +;; pade(%, 10, 10) -> +;; +;; [(x^10+110*x^9+5940*x^8+205920*x^7+5045040*x^6+90810720*x^5+1210809600*x^4 +;; +11762150400*x^3+79394515200*x^2+335221286400*x^+670442572800) +;; /(x^10-110*x^9+ +;; 5940*x^8-205920*x^7+5045040*x^6-90810720*x^5+1210809600*x^4-11762150400*x^3+ +;; 79394515200*x^2-335221286400*x^+670442572800)] +;; +;; The numerator and denominator have the same coefficients with +;; different signs on odd terms. so we note that num and den here can +;; be evaluated as f(x^2)+x*g(x^2) and f(x^2)-x*g(x^2) respectively using +;; half the number of multiplies. + +(defconstant f-exp-pade + '(1.0d0 5940.0d0 5045040.0d0 1.2108096d+9 + 7.93945152d+10 6.704425728d+11)) + +(defconstant g-exp-pade + '(110.0d0 205920.0d0 9.081072d+7 + 1.17621504d+10 3.352212864d+11 )) + +(defun exp-qd/pade-small (x) + ;; exp(x) = (f(x^2) + x*g(x^2))/(f(x^2) - x*g(x^2)) + (let* ((x^2 (sqr-qd x)) + (fx (poly-eval-qd-d x^2 f-exp-pade)) + (x*gx (mul-qd x (poly-eval-qd-d x^2 g-exp-pade)))) + (div-qd (add-qd fx x*gx) + (sub-qd fx x*gx)))) + +(defun exp-qd/pade (a) + ;; Same idea as in exp-qd/reduce, except we use our Pade + ;; approximation instead of a Taylor series. + (let* ((k 256) + (z (truncate (qd-0 (nint-qd (div-qd a +qd-log2+))))) + (r (div-qd-d (sub-qd a (mul-qd-d +qd-log2+ (float z 1d0))) + (float k 1d0))) + + (e (exp-qd/pade-small r))) + (scale-float-qd (npow e k) z))) + + ;; Some timing and consing tests. ;; ;; The tests are run using the following: @@ -1035,11 +1100,13 @@ ;; exp-qd/reduce 2.06 3.18 10.46 2.76 6.12 ;; expm1-qd/series 8.81 12.24 18.87 3.26 29.0 ;; expm1-qd/dup 5.68 4.34 18.47 3.64 18.78 +;; exp-qd/pade 1.53 ;; ;; Consing (MB) Sparc ;; exp-qd/reduce 45 45 638 44.4 45 ;; expm1-qd/series 519 519 1201 14.8 519 ;; expm1-qd/dup 32 32 1224 32.0 32 +;; exp-qd/pade 44 ;; ;; Speeds seem to vary quite a bit between architectures. ;; @@ -1066,17 +1133,22 @@ (fixnum n)) (let ((y +qd-zero+)) (declare (type %quad-double y)) - #+cmu (gc :full t) + #+cmu (ext:gc :full t) (format t "exp-qd/reduce~%") (time (dotimes (k n) (declare (fixnum k)) (setf y (exp-qd/reduce x)))) - #+cmu (gc :full t) + #+cmu (ext:gc :full t) + (format t "exp-qd/pade~%") + (time (dotimes (k n) + (declare (fixnum k)) + (setf y (exp-qd/pade x)))) + #+cmu (ext:gc :full t) (format t "expm1-qd/series~%") (time (dotimes (k n) (declare (fixnum k)) (setf y (expm1-qd/series x)))) - #+cmu (gc :full t) + #+cmu (ext:gc :full t) (format t "expm1-qd/duplication~%") (time (dotimes (k n) (declare (fixnum k)) @@ -1084,6 +1156,87 @@ )) +#|| + +;; This bit of code is meant to be run with Clisp, which has arbitrary +;; precision long floats. +;; +;; We use this to compare the accuracy of exp-qd/reduce and +;; exp-qd/pade against Clisp's exp implementation. +;; +;; Some results: +;; +;; k = 159000 +;; pade min = 4.733L-68 +;; pade max = 5.153L-61 +;; pade sumsq = 7.899L-119 +;; red min = 3.838L-69 +;; red max = 4.318L-62 +;; red sumsq = 6.076L-120 +;; +;; So, worst case error is: +;; +;; exp-qd/reduce: 4.318L-62 (203.8 bits) +;; exp-qd/pade: 5.153L-61 (200.3 bits) +;; +;; We lose 3 bits for Pade. exp-qd/reduce is "missing" 8 bits of +;; accuracy. Bummer. + +(in-package #:oct) + +;; 240 bit fractions for long-floats. I hope this is accurate enough +;; compared against our 212 bit fractions for quad-doubles. +(setf (ext:long-float-digits) 240) +(defun compare-exp (n &optional verbose) + (let ((pade-max -1L0) + (pade-min 1L100) + (pade-sum 0L0) + (red-max -1L0) + (red-min 1L100) + (red-sum 0L0)) + (dotimes (k n) + (let* ((f (random (make-qd 1/256))) + (pade-exp (octi::exp-qd/pade (qd-value f))) + (red-exp (octi::exp-qd/reduce (qd-value f)))) + (flet ((qd-to-lf (qd) + (octi:with-qd-parts (q0 q1 q2 q3) + qd + (+ (float q0 1L0) + (float q1 1L0) + (float q2 1L0) + (float q3 1L0))))) + (let* ((lf (qd-to-lf (qd-value f))) + (l-exp (exp lf)) + (pade-exp-lf (qd-to-lf pade-exp)) + (red-exp-lf (qd-to-lf red-exp)) + (pade-diff (abs (- l-exp pade-exp-lf))) + (red-diff (abs (- l-exp red-exp-lf)))) + (when verbose + (format t "f = ~S~%" f) + (format t "pade-exp = ~A~%" (make-instance 'qd-real :value pade-exp)) + (format t "l-exp = ~A~%" l-exp) + (format t "diff = ~A~%" pade-diff)) + (when (and (plusp k) (zerop (mod k 1000))) + (format t "k = ~A~%" k) + (format t "pade min = ~15,3g~%" pade-min) + (format t "pade max = ~15,3g~%" pade-max) + (format t "pade sumsq = ~15,3g~%" pade-sum) + (format t "red min = ~15,3g~%" red-min) + (format t "red max = ~15,3g~%" red-max) + (format t "red sumsq = ~15,3g~%" red-sum)) + (setf pade-max (max pade-max pade-diff)) + (setf pade-min (min pade-min pade-diff)) + (incf pade-sum (* pade-diff pade-diff)) + (setf red-max (max red-max red-diff)) + (setf red-min (min red-min red-diff)) + (incf red-sum (* red-diff red-diff))) + ))) + (values pade-max pade-min pade-sum + red-max red-min red-sum + ))) + +||# + ;; (time-log #c(3w0 0) 50000) ;; ;; Time (s) Sparc PPC x86 PPC (fma) Sparc2 @@ -1147,33 +1300,33 @@ (fixnum n)) (let ((y +qd-zero+)) (declare (type %quad-double y)) - #+cmu (gc :full t) + #+cmu (ext:gc :full t) (format t "log-qd/newton~%") (time (dotimes (k n) (declare (fixnum k)) (setf y (log-qd/newton x)))) - #+cmu (gc :full t) + #+cmu (ext:gc :full t) (format t "log1p-qd/duplication~%") (time (dotimes (k n) (declare (fixnum k)) (setf y (log1p-qd/duplication x)))) - #+cmu (gc :full t) + #+cmu (ext:gc :full t) (format t "log-qd/agm~%") (time (dotimes (k n) (declare (fixnum k)) (setf y (log-qd/agm x)))) - #+cmu (gc :full t) + #+cmu (ext:gc :full t) (format t "log-qd/agm2~%") (time (dotimes (k n) (declare (fixnum k)) (setf y (log-qd/agm2 x)))) - #+cmu (gc :full t) + #+cmu (ext:gc :full t) (format t "log-qd/agm3~%") (time (dotimes (k n) (declare (fixnum k)) (setf y (log-qd/agm3 x)))) - #+cmu (gc :full t) + #+cmu (ext:gc :full t) (format t "log-qd/halley~%") (time (dotimes (k n) (declare (fixnum k)) @@ -1217,17 +1370,17 @@ (fixnum n)) (let ((y +qd-zero+) (one +qd-one+)) - #+cmu (gc :full t) + #+cmu (ext:gc :full t) (format t "atan2-qd/newton~%") (time (dotimes (k n) (declare (fixnum k)) (setf y (atan2-qd/newton x one)))) - #+cmu (gc :full t) + #+cmu (ext:gc :full t) (format t "atan2-qd/cordic~%") (time (dotimes (k n) (declare (fixnum k)) (setf y (atan2-qd/cordic x one)))) - #+cmu (gc :full t) + #+cmu (ext:gc :full t) (format t "atan-qd/duplication~%") (time (dotimes (k n) (declare (fixnum k)) @@ -1260,12 +1413,12 @@ (declare (type %quad-double x) (fixnum n)) (let ((y +qd-zero+)) - #+cmu (gc :full t) + #+cmu (ext:gc :full t) (format t "tan-qd/cordic~%") (time (dotimes (k n) (declare (fixnum k)) (setf y (tan-qd/cordic x)))) - #+cmu (gc :full t) + #+cmu (ext:gc :full t) (format t "tan-qd/sincos~%") (time (dotimes (k n) (declare (fixnum k))

1 0

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

18 Oct '07

Update of /project/oct/cvsroot/oct In directory clnet:/tmp/cvs-serv27836 Modified Files: qd-fun.lisp Log Message: Compare against pi/4 quad-double, not double-float. --- /project/oct/cvsroot/oct/qd-fun.lisp 2007/10/17 03:44:12 1.89 +++ /project/oct/cvsroot/oct/qd-fun.lisp 2007/10/18 14:38:56 1.90 @@ -678,7 +678,7 @@ ;; We want the remainder part to be <= pi/4 because the trig ;; functions want that. So if the fraction is too big, adjust ;; it, and mod value. - (if (<= (abs (qd-0 f)) (/ pi 4)) + (if (qd-<= (abs-qd f) +qd-pi/4+) (values mod f) (values (mod (1+ mod) 4) (sub-qd f +qd-pi/2+)))))) @@ -723,7 +723,7 @@ +qd-pi/2+))) ;; If we did this right, (ash prod (- frac-bits)) should be 2 ;; bits long at most. - (if (<= (abs (qd-0 f)) (/ pi 4)) + (if (qd-<= (abs-qd f) +qd-pi/4+) (values mod f) (values (mod (1+ mod) 4) (sub-qd f +qd-pi/2+))))))

1 0

← Newer
1
...
9
10
11
12
13
14
15
...
21
Older →

Jump to page:

2024

2023

2022

2021

2020

2019

2018

2017

2016

2015

2014

2013

2012

2011

2010

2009

2008

2007

oct-cvs