Re: Same .m file: different results with different versions of Octave

[Jaroslav Hajek]

On Mon, Apr 19, 2010 at 5:44 PM, Judd Storrs <address@hidden> wrote:
> I'm not sure that's really a good fix for the problem. GSL's approach
> is fundamentally superior because it avoids cancellation entirely by
> factoring the results into a product. All of the additions in GSL's
> implementation are between probably positive values that cannot have
> cancellation problems. glibc's algorithim is comparatively naive and
> should just be replaced but I haven't had a chance to verify the math
> behind the GSL terms yet.
>

I verified both the GSL and libc formulas, they're indeed correct. I
tried the GSL implementation; it seems to be about 40% slower on
doubles, somewhat more on singles. I hope one can do better. Also,
note that the libc implementation handles the NaN + 0*i case better.
The libc formula suffers a possible cancellation problem when re is
tiny and im is close to (k+1/2)*pi, but I didn't yet succeed in
figuring out an instance. Maybe there is an approach that takes the
best of both.

-- 
RNDr. Jaroslav Hajek, PhD
computing expert & GNU Octave developer
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz