specfunc error analysis - generated and propagated errors.

[Jonathan Underwood <j dot underwood at open dot ac dot uk>]

Dear All,

In an effort to fix up the error analysis in the wigner coefficient code I 
posted recently, I started to look at how the error analysis is 
implemented for GSL in the special functions (specfunc directory) and came 
to the following conclusions and questions. The reason I'm posting these 
is in the hope that my conclusions will be confirmed or corrected, and the 
questions answered :). I'll split the discussion up into generated errors, 
and propagated errors.

(1) Generated error.
---------------------
By which i mean the error from carrying out a function on a variable. For 
most of the gsl_sf functions these are calculated and returned, however, 
I'm a bit unsure  as to how things have been standardized for GSL for 
addition, subtraction, multiplication and division i.e. for exact A and B, 
what is the generated error for the operations A {+,-,*,/} B.

(a) Multiplication: Looking at coupling.c, it looks like the _relative_ 
generated error for multplication (A*B) is (2 * GSL_DBL_EPSILON). [See 
lines 73 and 74 of coupling.c]. However, I then looked at legendre_poly.c, 
as this contains some very simple algebra. Looking at line 84, I see that 
my assertion must be wrong. But then, comparing line 122, which is the 
same algebraic expression as line 84, gives an entirely different 
generated error analysis. What gives?

(b) Addition: Again, the relative generated error seems to be (2 * 
GSL_DBL_EPSILON) - see eg. line 177 of coupling.c

(c) Subtraction: Absolute generated error seems to be (2*GSL_DBL_EPSILON) 
* (|A|+|B|) [this seems very odd]

(d) Division: Didn't consider this, but it would be nice to know.

It would be really helpful if someone could clarify on the issue of 
generated errors, with examples. This should eventually go in the design 
document (I'm happy to collate the text and add it to the documentation).


(2) Propagated Error
---------------------
This is simpler
(a) Addition and subtraction: Algebraically add the absolute errors in all 
variables being added together to get the result's absolute error. eg. for 
A with absolute error a, and B with absolute error b, the error in (A +/- 
B) will be a+b. As Brian said, we're not adding errors in quadrature here 
(a simplification).

(b) Multiplication and division: the resulting relative error is the sum 
of the relative errors of the values, i.e. for A*B, the relative error in 
the answer is (a/A + b/B), and the absolute error in the answer is 
|AB|*(|a/A| + |b/B|). Ignoring terms of (ab), can then write the resultant 
absolute error as a|A| +b|B|. This is extended to A*B*C, where the 
absolute error in the answer will be a|BC| + b|AC| + c|AB|. etc.

Further points
---------------
Brian in a previous email alluded to keeping the relative errors signed 
such as to allow for error cancellation during propagation, however I 
don't see any implementation of that.

Related to this, I wonder what the error analysis is actually trying to 
achieve. It seems to be giving an estimation of the worst case error 
resulting from roundoff error, but I don't think it is in anyway checking 
or estimating truncation error (is it even possible to have a general 
pescription for estimating truncation error numerically?). When I say 
"worst case", I mean it doesn't seem to allow for error cancellation, and 
assumes all steps contribute maximum error (there will of course be a 
distribution of error).

As a final design point - heretically, I wonder about the decision to be 
always calculating the error - this puts a big overhead on each 
calculation as it roughly doubles the number of operations. I realize that 
by always giving the user the option to examine the error (s)he can then 
write new functions and propagate the error accordingly. I'd reckon that 
> 98 percent of the time, noone uses the _e functions however. Would it not 
be better to have error calculation only done when the library is compiled 
with a certain flag set?

Finally, I apologize for all ignorance displayed by myself herein.

Thanks in advance,
Jonathan