At Fri, 15 Feb 2008 00:00:55 +0100, Frank Reininghaus wrote: > Maybe these (at least for me) surprising observations could be used to > isolate the cause of the bug. Does anyone have a good idea?
It seems that in double precision the error term (delta_frac - 1) is exactly zero at some point, so the continued fraction terminates. It looks like this happens before the desired accuracy is actually reached though, so it seems to be undesirable. In extended precision that doesn't happen, the loop continues up to its maximum iteration limit and NAN is returned to signal that. If you add #include <gsl/gsl_ieee_utils.h> and gsl_ieee_env_setup(); then run with GSL_IEEE_MODE=double-precision ./a.out you can see the difference.
