Mail Archives: djgpp/2004/02/17/08:01:21
Ethan Rosenberg <ethros AT earthlink DOT net> wrote:
> The following are code FRAGMENTS from the program. My feelings are that
> the error is in the integration routine. etot is defined as FLOAT.
The problem at hand is rather probably one of numerical accuracy of
your code, i.e. rounding errors. It may seem surprising that these
should depend on the operating system, while the CPU is kept the same,
but that can, indeed, happen.
But such problems are almost impossible to debug unless both actual,
complete source code that still shows the problematic behaviour and a
person with experience in investigating such matters meet each other
in the same time and place for a while. As long as you don't have the
expertise, and we don't have a reproducible example, I suspect this is
leading nowhere.
Please try to cut down your program to the smallest one you can come
up with that still demonstrates the problem:
1) Put in *lots* of printouts, to find the first calculation that
does report a different result on different operating systems.
2) throw away everything after that point, and all function only
called after by what you threw away (gcc -Wall -W will report unused
functions, if it's a single source file and all file-scope function
other than main() marked "static").
3) Get rid of all statements that don't affect values used to compute
the number that shows up as different. Careful though: removing such
lines *may* change the outcome ...
Try to repeat these steps until you don't find any more pieces that
can be discarded. Post the result here, if it's short, or on a web
site if it's not. Or, if all else fails, mail it to interested
parties directly.
> /*
> * If you want double-precision floats used everywhere,
> * define this to be "double"
> */
> #define FLOAT float
I agree with Eli: this choice rather likely is an important part of
the problem. I strongly suspect all the differences will vanish if
you change this from "float" to "long double". Going to "double" may
work, too, but "long double" is the real test.
I think it's actually even worse than Eli realized: A float has only
about 7 significant digits, but you're summing up 16384 of them,
eventually. I.e. even if all input numbers were of the same order of
magnitude, you're in danger of killing at least 4 significant decimals
just by the sheer number of terms you sum up. I.e. your sums may
easily have only three significant digits left to them, *if*
intermediate results are actually restricted to be single-precision
floats.
An important test for this kind of problem is to compile your program
using GCC flag -ffloat-store and see if that effects the result. I'm
almost willing to bet it will. The effect of this is similar, but not
quite the same as the switch 'long double' mentioned before.
> /*
> * Calculate the hypotenuse of a right triangle. The function in math.h
> * assumes the the input and return values are double. This uses FLOAT
> */
> FLOAT hypot2(FLOAT a, FLOAT b)
> {
> FLOAT c,asq,bsq;
> double aa,bb,cc;
> asq = a*a;
> bsq = b*b;
> aa = (double)asq;
> bb = (double)bsq;
> cc = sqrt(aa+bb);
> c = (FLOAT)sqrt((double)(a*a) + (double)(b*b));
> return c;
> }
This seems mangled. Why compute all those asq, aa and whatever, if
you're not going to use them to comput the actual output?
--
Hans-Bernhard Broeker (broeker AT physik DOT rwth-aachen DOT de)
Even if all the snow were burnt, ashes would remain.
- Raw text -