[SciPy-user] round off errors

[Robert Kern]

Ryan Krauss wrote:

> How will SciPy handle sinh(b)+ (or -) sin(b) when sinh(b) is very
> large and sin(b) is obviously between +/-1?  As I understand it, the
> IEEE standard requires exact rounding for this operation which seems
> difficult and costly to me.

Scipy doesn't really do anything special here. We just add/subtract the values
in C. After that, it's the FPU's job.

One thing to note is that the C functions that compute sinh() and sin() are
usually defined in double precision only. If you give them a single precision
array as input, the function will upcast it to a double, do the computation in
double precision, and downcast the answer to a single precision value again.
Also, some processors may do addition and subtraction with 80-bit extended
precision in the intermediate stages.

Raymond Hettinger wrote a floating point simulator to help explore the silliness
of floating point math. You might find it useful.

  http://aspn.activestate.com/ASPN/Cookbook/Python/Recipe/265894

One thing you might want to consider is breaking apart the sums into lists of
values. You can then do various things to make the total sum behave better.
Google "Kahan summation."

You may also want to look into extended/unlimited precision libraries. The SAGE
environment has bindings to lots of good stuff for this kind of work.

  http://sage.sourceforge.net/

-- 
Robert Kern
robert.kern at gmail.com

"In the fields of hell where the grass grows high
 Are the graves of dreams allowed to die."
  -- Richard Harter