X-Authentication-Warning: delorie.com: mail set sender to djgpp-bounces using -f
From: Hans-Bernhard Broeker <broeker AT physik DOT rwth-aachen DOT de>
Newsgroups: comp.os.msdos.djgpp
Subject: Re: Cross Platform Incompatabilites? - code fragments
Date: 3 Mar 2004 10:40:36 GMT
Organization: Aachen University of Technology (RWTH)
Lines: 67
Message-ID: <c24cn4$lk$1@nets3.rz.RWTH-Aachen.DE>
References: <3 DOT 0 DOT 1 DOT 16 DOT 20040216231142 DOT 38e7a7cc AT earthlink DOT net> <3 DOT 0 DOT 1 DOT 16 DOT 20040302202549 DOT 3957c428 AT earthlink DOT net>
NNTP-Posting-Host: ac3b07.physik.rwth-aachen.de
X-Trace: nets3.rz.RWTH-Aachen.DE 1078310436 692 137.226.33.205 (3 Mar 2004 10:40:36 GMT)
X-Complaints-To: abuse AT rwth-aachen DOT de
NNTP-Posting-Date: 3 Mar 2004 10:40:36 GMT
To: djgpp AT delorie DOT com
DJ-Gateway: from newsgroup comp.os.msdos.djgpp
Reply-To: djgpp AT delorie DOT com

Ethan Rosenberg <ethros AT earthlink DOT net> wrote:

[...]
>     HBR/long double
>               min           delta          etot
>               0.45     135.183805  650175.488059
>               3.15      -3.106527    1270.050546

>     SARAH/long double
>                 min             delta          etot
>               0.45     135.183805  650175.488059
>               3.15      -3.106527    1270.050546

Note that 'long double' is the only type in which both machines return
the same result?  I'm not completely convinced yet, but this seems
like it supports my assumptions: FPU register precision is part of
your problem.  Long double has the same precision as the FPU
registers, so changes to the FPU default precision shouldn't affect
them as they do the other floating-point types.

> The problem with the long double calculations, I think, is that the sqrt
> function assumes double for input.

That shouldn't too bad a problem.  Doesn't DJGPP's full-scale libm
come with a sqrtl() for long double arguments?  Or was that only the
case for the full CEPHES lib?

> The FFT routine which I have not modified at all is:

> Recognizing the the problem may be from overflow of the float, certain
> changes were made in the program:

Overflow would have manifested itself in different ways.  After you
turned on SIGFPE generation, you would have got right-out crashes of
the program.  Otherwise, you would be seeing Inf and NaN among your
outputs.

This really looks more like catastrophic loss of precision.  

Did you try -mfloat-store yet?


>     qsort(temp3, jj, sizeof(FLOAT), mycomp);
>     for(i = 0; i < jj; i++)
>         sum += 2*temp3[i];

That's not quite how the "sum smallest items first" trick I outlined is
supposed to be used.

You're supposed to have a loop like this:

     while(number of terms > 1)
       find two smallest elements (by magnitude) of the list, i and k:
       remove i and k from the list
       put their sum, i+k into it instead

The idea is to avoid summing any two items of vastly different
magnitudes, if that's at all possible.  

And there's still the odd chance that the problem you're having isn't
actually with the floating point data at all, but with some random
garbage written somewhere by some part of the source you never showed
us.

-- 
Hans-Bernhard Broeker (broeker AT physik DOT rwth-aachen DOT de)
Even if all the snow were burnt, ashes would remain.