Hello,
I was looking at some of W. Kahan's papers on round-off errors and
found the following test for estimation of the roundoff error in a computer.
In some notes on his web page entitled
"OLD Notes on Errors and Equation-Solving"
he states the following
If y fop z denotes the floating point representation of the mathematical
operation y op z where
op is one of +,-,*,/, then
y fop z = (y op z)/(1-a)
where
abs(a) < abs( ( ((4.0/3.0 rounded) - 1.0)*3.0 - 1.0 )
There is a difference between the stated estimate for standard double
precision and the output using 'bfloat'.
Kahan mentions that his estimate is an over-estimataion, so there is no
contradiction.
My question is whether the computed outcome with 'bfloat' is also a
correct upper estimate for 'a'. Note: I did not check the mathematics
involved. I thought someone (maybe RJF) would know the answer immediately.
Here are some outputs (maxima 5.18.1 with cmucl).
(%i5) fpprec;
(%o5) 16
(%i6) abs(((4.0/3.0) - 1.0)*3.0 -1.0);
(%o6) 2.220446049250313e-16
Here is the computation using 'bfloats'.
(%i7) abs(((bfloat(4.0)/bfloat(3.0)) - bfloat(1.0))*bfloat(3.0)
-bfloat(1.0));
(%o7) 2.775557561562891b-17
Is this smaller number an accurate estimate ?
More generally, does the method produce valid estimations for any
precision?
In particular, are the following upper bounds valid?
%i10) fpprec: 32;
(%o10) 32
(%i11) abs(((bfloat(4.0)/bfloat(3.0)) - bfloat(1.0))*bfloat(3.0)
-bfloat(1.0));
(%o11) 3.0814879110195773648895647081359b-33
(%i12) fpprec: 64;
(%o12) 64
(%i13) abs(((bfloat(4.0)/bfloat(3.0)) - bfloat(1.0))*bfloat(3.0)
-bfloat(1.0));
(%o13) 3.798227098303919498989296907824782861688386333447977986511911996b-65
TIA,
-sen