Mail Archives: djgpp/2004/03/03/21:16:09
On 3 Mar 2004 10:40:36 GMT in comp.os.msdos.djgpp, Hans-Bernhard
Broeker <broeker AT physik DOT rwth-aachen DOT de> wrote:
>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.
Treating the array as a heap would allow this to be done efficiently.
See Knuth, Sedgewick, or any algorithms book, for details on heaps.
Heap creation replaces the qsort. Keeping the smallest items at the
end of the array/heap simplifies reheaping and reduces the work
required after removing the two smallest items and inserting their
sum. Ensure the comparison operations are the right way round to keep
the array/heap in largest first order.
--
Thanks. Take care, Brian Inglis Calgary, Alberta, Canada
Brian DOT Inglis AT CSi DOT com (Brian dot Inglis at SystematicSw dot ab dot ca)
fake address use address above to reply
- Raw text -