Re: [Bug-gsl] diff.c (numerical differentiation)

[spasmous]

Hi Brian,

Sorry that may also be foiled if the CPU has high precision registers.
I think it's necessary to add a volatile variable to force the compiler to
use the same precision throughout.


  volatile double y;
  ...
  h = sqrt (GSL_SQRT_DBL_EPSILON / (2.0 * a2));
  ...
  // ensure difference between (x+h) and h are representable in double
  y = x+h;
  h = y-x;

  *result = (GSL_FN_EVAL (f, x + h) - GSL_FN_EVAL (f, x)) / h;


Thanks!



On Tue, Apr 8, 2008 at 1:55 PM, Brian Gough <address@hidden>
wrote:

> At Sun, 6 Apr 2008 13:29:41 -0700,
> spasmous wrote:
> >   h = sqrt (GSL_SQRT_DBL_EPSILON / (2.0 * a2));
> >   ...
> >   h = (x+h)-x; // ensure difference between (x+h) and x is exactly h
> >   *result = (GSL_FN_EVAL (f, x + h) - GSL_FN_EVAL (f, x)) / h;
>
> Thanks for the suggestion.
>
> --
> Brian Gough
>