[cmucl-commit] [git] CMU Common Lisp branch master updated. snapshot-2014-08-34-g1345f3b
Raymond Toy
rtoy at common-lisp.net
Wed Aug 27 03:30:14 UTC 2014
This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "CMU Common Lisp".
The branch, master has been updated
via 1345f3b146c9840447fa7a943daf85f111011a7f (commit)
from df548f6e69bc4fa0aec3825d10074135c9e4ac03 (commit)
Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.
- Log -----------------------------------------------------------------
commit 1345f3b146c9840447fa7a943daf85f111011a7f
Author: Raymond Toy <toy.raymond at gmail.com>
Date: Tue Aug 26 20:29:56 2014 -0700
Use the two-prod algorithm from crlibm documentation.
This handles overflows better and
(c::two-prod 1.7976931214684583d308 (1+ (scale-float 1d0 -28)))
doesn't signal an overflow like the old algorithm.
diff --git a/src/compiler/float-tran.lisp b/src/compiler/float-tran.lisp
index 7da6b95..ccb6e06 100644
--- a/src/compiler/float-tran.lisp
+++ b/src/compiler/float-tran.lisp
@@ -2289,22 +2289,9 @@
(kernel:%make-double-double-float hi lo))))
(declaim (maybe-inline split))
-;; This algorithm is the version given by Yozo Hida. It has problems
-;; with overflow because we multiply by 1+2^27.
-;;
-;; But be very careful about replacing this with a new algorithm. The
-;; values computed here are very important to get the rounding right.
-;; If you change this, the rounding may be different, which will
-;; affect other parts of the algorithm.
-;;
-;; I (rtoy) tried a different algorithm that split the number in two
-;; as described, but without overflow. However, that caused
-;; -9.4294948327242751340284975915175w0/1w14 to return a value that
-;; wasn't really close to -9.4294948327242751340284975915175w-14.
-;;
-;; This also means we can't print numbers like 1w308 with the current
-;; printing algorithm, or even divide 1w308 by 10.
-#+nil
+;; See Listing 2.6: Mul12 in "CR-LIBM: A library of correctly rounded
+;; elementary functions in double-precision". Also known as Dekker's
+;; algorithm.
(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
@@ -2315,68 +2302,61 @@
(a-lo (- a a-hi)))
(values a-hi a-lo)))
-;; +split-limit+ is the largest number for which Yozo's algorithm
-;; still works. Basically we want a*(1+2^27) <=
-;; most-positive-double-float < 2^1024. Therefore, a < 2^1024/(1+2^27)
-;; If we calculate that, we get a = 1.3393857490036326d300. A quick
-;; test shows that this would cause overflow, but previous float would
-;; not. This is the value we want.
-(defconstant +split-limit+
- (scale-float (/ (float (1+ (expt 2 27)) 1d0)) 1024))
+;; Values used for scaling in two-prod. These are used to determine
+;; if SPLIT might overflow so the value (and result) can be scaled to
+;; prevent overflow.
+(defconstant +two970+
+ (scale-float 1d0 970))
-(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
- a-lo contains the lower 26 bits."
- (declare (double-float a)
- (optimize (speed 3)
- (inhibit-warnings 3)))
- ;; This splits the number a into 2 parts of 27 and 26 bits each, but
- ;; the parts are, I think, supposed to be properly rounded in an
- ;; IEEE fashion.
- ;;
- ;; 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) +split-limit+)
- ;; 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
- ;; are very small, this algorithm produces different numbers
- ;; 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)))
- (a-lo (- a a-hi)))
- (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)))
- (values a-hi a-lo))))
+(defconstant +two53+
+ (scale-float 1d0 53))
+(defconstant +two-53+
+ (scale-float 1d0 -53))
(declaim (inline two-prod))
+
+;; This is essentially the algorithm given by Listing 2.7 Mul12Cond
+;; given in "CR-LIBM: A library of correctly rounded elementary
+;; functions in double-precision".
#-ppc
(defun two-prod (a b)
_N"Compute fl(a*b) and err(a*b)"
- (declare (double-float a b))
- (let ((p (* a b)))
- (multiple-value-bind (a-hi a-lo)
- (split a)
- ;;(format t "a-hi, a-lo = ~S ~S~%" a-hi a-lo)
- (multiple-value-bind (b-hi b-lo)
- (split b)
- ;;(format t "b-hi, b-lo = ~S ~S~%" b-hi b-lo)
- (let ((e (+ (+ (- (* a-hi b-hi) p)
- (* a-hi b-lo)
- (* a-lo b-hi))
- (* a-lo b-lo))))
- (locally
- (declare (optimize (inhibit-warnings 3)))
- (values p e)))))))
+ (declare (double-float a b)
+ (optimize (speed 3)))
+ ;; If the numbers are too big, scale them done so SPLIT doesn't overflow.
+ (multiple-value-bind (aa bb)
+ (values (if (> a +two970+)
+ (* a +two-53+)
+ a)
+ (if (> b +two970+)
+ (* b +two-53+)
+ b))
+ (let ((p (* aa bb)))
+ (declare (double-float p)
+ (inline split))
+ (multiple-value-bind (aa-hi aa-lo)
+ (split aa)
+ ;;(format t "aa-hi, aa-lo = ~S ~S~%" aa-hi aa-lo)
+ (multiple-value-bind (bb-hi bb-lo)
+ (split bb)
+ ;;(format t "bb-hi, bb-lo = ~S ~S~%" bb-hi bb-lo)
+ (let ((e (+ (+ (- (* aa-hi bb-hi) p)
+ (* aa-hi bb-lo)
+ (* aa-lo bb-hi))
+ (* aa-lo bb-lo))))
+ (declare (double-float e))
+ (locally
+ (declare (optimize (inhibit-warnings 3)))
+ ;; If the numbers was scaled down, we need to scale the
+ ;; result back up.
+ (when (> a +two970+)
+ (setf p (* p +two53+)
+ e (* e +two53+)))
+ (when (> b +two970+)
+ (setf p (* p +two53+)
+ e (* e +two53+)))
+ (values p e))))))))
#+ppc
(defun two-prod (a b)
@@ -2402,6 +2382,27 @@
(values q (+ (+ (- (* a-hi a-hi) q)
(* 2 a-hi a-lo))
(* a-lo a-lo)))))))
+(defun two-sqr (a)
+ _N"Compute fl(a*a) and err(a*b). This is a more efficient
+ implementation of two-prod"
+ (declare (double-float a))
+ (let ((aa (if (> a +two970+)
+ (* a +two-53+)
+ a)))
+ (let ((q (* aa aa)))
+ (declare (double-float q)
+ (inline split))
+ (multiple-value-bind (a-hi a-lo)
+ (split aa)
+ (locally
+ (declare (optimize (inhibit-warnings 3)))
+ (let ((e (+ (+ (- (* a-hi a-hi) q)
+ (* 2 a-hi a-lo))
+ (* a-lo a-lo))))
+ (if (> a +two970+)
+ (values (* q +two53+)
+ (* e +two53+))
+ (values q e))))))))
#+ppc
(defun two-sqr (a)
-----------------------------------------------------------------------
Summary of changes:
src/compiler/float-tran.lisp | 143 +++++++++++++++++++++---------------------
1 file changed, 72 insertions(+), 71 deletions(-)
hooks/post-receive
--
CMU Common Lisp
More information about the cmucl-commit
mailing list