[oct-cvs] Oct commit: oct qd-dd.lisp

[rtoy]

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)