Page 1 of 2
Numerical Analysis
Posted: Fri Nov 15, 2013 5:11 pm
by Robert
I have a function, say 'Foo', that returns
with u & v declared as REAL.
My application has the line:
Code: Select all
IF u * u + v * v <= 1. THEN w := Foo (u, v) ELSE ...
This has worked fine for years, and I assumed that it was 'safe'. Yesterday it crashed (Square root of a negative number).
My first guess was it was due to some operations being 64 bit, and others leaving 80-bit partial results on the floating point stack.
However the first experiment I tried was to rewrite Foo as
ie to simply change the evaluation order.
That change stopped the TRAP with yesterday's data, however I am not sure if the code is now completely reliable, or even if this change will cause some pairs of values that did work to now fail, ie it may have made the problem worse!
Does anyone have any useful insight?
Incidentally I have reasons to not to like the obvious 'hacks':
or
Code: Select all
Math.Sqrt (MAX (1. - u * u - v * v), 0.)
Regards
Robert
PS - Ivan: Please don't answer this, give someone else a chance to get involved with the forum!
Re: Numerical Analysis
Posted: Tue Nov 19, 2013 7:25 am
by Josef Templ
I would avoid the common sub-expression u * u + v * v by evaluating
it only once and storing it in an auxiliary variable.
This is (1) faster and (2) more predictable with respect to floating
point precision.
- Josef
Re: Numerical Analysis
Posted: Tue Nov 19, 2013 9:03 am
by manumart1
Instead of
you can use
Although they seem equivalent, they are not because of floating point issues.
I have made some test and found some values which show that:
Code: Select all
MODULE TestFloatPoint;
IMPORT StdLog, Math;
PROCEDURE Do*;
CONST ONE = 0.99999999999999999;
VAR u, v: REAL;
k, j, tot, hits: INTEGER;
BEGIN
StdLog.String("ONE = "); StdLog.Real(ONE); StdLog.Ln;
tot := 0;
hits := 0;
FOR k := 1 TO 99 DO
FOR j := 1 TO k DO
u := Math.Sqrt((ONE / k) * j);
v := Math.Sqrt((ONE / k) * (k - j));
INC(tot, 1);
IF u * u + v * v <= 1. THEN
IF 1. - u * u - v * v < 0 THEN
INC(hits, 1);
StdLog.String(" *************************** "); StdLog.Ln;
StdLog.String(" u = "); StdLog.Real(u); StdLog.Ln;
StdLog.String(" v = "); StdLog.Real(v); StdLog.Ln;
StdLog.String(" u * u + v * v = "); StdLog.Real(u * u + v * v); StdLog.Ln;
StdLog.String(" 1 - u * u - v * v = "); StdLog.Real(1 - u * u - v * v); StdLog.Ln;
StdLog.String(" *************************** "); StdLog.Ln
END
END
END
END;
StdLog.String("END. tot = "); StdLog.Int(tot); StdLog.String(" hits = "); StdLog.Int(hits); StdLog.Ln
END Do;
END TestFloatPoint.
***************************
u = 0.9925092578236596
v = 0.1221694443563052
u * u + v * v = 1.0
1 - u * u - v * v = -2.456395547037471E-20
***************************
Regards,
Manuel
Re: Numerical Analysis
Posted: Tue Nov 19, 2013 11:46 am
by gerard94
Robert,
I'd like to think that, whenever a <= 1, then 1 - a >= 0. So, if you replace a by u * u + v * v, then 1 - (u * u + v * v) >= 0.
I've changed 1 - u * u - v * v by 1 - (u * u + v * v) in manumart1's example, and there is no more hit.
Cheers.
Gérard
Re: Numerical Analysis
Posted: Tue Nov 19, 2013 10:27 pm
by Robert
My first guess was it was a 64-bit / 80-bit issue, but here is an example using 3-digit base 10 floating point:
The Test:
0.501 + 0.502 <= 1
becomes
1.00(3) rounded to 1.00 <= 1
becomes
TRUE.
The Expression:
Sqrt (1.00 - 0.501 - 0.502) ~> Sqrt (0.499 - 0.502) ~> Sqrt (-0.003) ~> TRAP!
It is the order of evaluation that matters.
But is Sqrt (1.00 - (0.501 + 0.502)) safe?
If I 'model' the 80-bit floating point register as a 4 decimal digit intermediate result this last expression remains:
Sqrt (1.00 - 1.003 (not yet rounded!)) ~> Sqrt (-0.003) ~> TRAP.
Does the BlackBox compiler round intermediate results to 64-bit prior to the comparison "<= 1.0" ?
Robert
Re: Numerical Analysis
Posted: Wed Nov 20, 2013 9:01 am
by manumart1
Robert wrote:But is Sqrt (1.00 - (0.501 + 0.502)) safe?
I think yes,
Sqrt (1.00 - (0.501 + 0.502)) is more safe than
Sqrt (1.00 - 0.501 - 0.502)
I say this because (as Gérard has pointed) in my previous program to find anomalous values of u and v, I see this:
Code: Select all
...
IF u * u + v * v <= 1. THEN
IF 1. - (u * u) - (v * v) < 0 THEN
INC(hits, 1);
...
--> Anomalous values found
Code: Select all
...
IF u * u + v * v <= 1. THEN
IF 1. - (u * u + v * v) < 0 THEN
INC(hits, 1);
...
--> Anomalous values NOT found
So perhaps the expression
(u * u + v * v) is rounded from 80 bits (CPU registers) to 64 bits (memory cells) prior to subtract from 1.
On the other hand, perhaps the expression
1. - (u * u) - (v * v) is is totally computed using only in CPU registers (80 bits).
I don't know.
Regards
Re: Numerical Analysis
Posted: Wed Nov 20, 2013 6:46 pm
by Rex Couture
I think manumart1 has the right idea (or at least one right idea), which is to test the value of the actual, entire parameter, rather than test part of the expression (u*u +v*v). For example, you might write
IF 1. -u*u -v*v >= 0. THEN w := Foo (u, v) ELSE ....
Even so, I'm not entirely certain that the expression will always be evaluated with the same precision in both the calling module and the procedure.
To guard against that problem, I really like Robert's hack:
RETURN Math.Sqrt (MAX (1. - u*u -v*v), 0.)
(or some equivalent). This would work if the problem is structured so that 1 -u*u -v*v is never negative (except as a numerical accident). However, it would return a very wrong answer in the case of a substantially negative number.
This brings up a related question on how to structure such a module. Do you put the negative number test in the calling module or the called procedure, or both? Each has advantages and disadvantages. If you put it in the calling module, any expression using the function could get messy. Furthermore, the procedure would lack generality. But if you put the test in the procedure, the procedure will have unavoidable overhead, which could cause a penalty if it's in an inner loop.
Re: Numerical Analysis
Posted: Wed Nov 20, 2013 8:36 pm
by Robert
The following Module:
- Does NOT TRAP on the ASSERT 30
- Does TRAP on the ASSERT 40.
Code: Select all
MODULE RdcRounding;
IMPORT Math;
PROCEDURE MyPi () : REAL;
BEGIN
RETURN Math.Pi ()
END MyPi;
PROCEDURE TestEQ (a, b : REAL);
BEGIN
ASSERT (a = b, 30)
END TestEQ;
PROCEDURE Do*;
VAR
pi : REAL;
BEGIN
pi := MyPi ();
TestEQ (pi, MyPi ());
ASSERT (pi = MyPi (), 40)
END Do;
END RdcRounding.
From this I draw three conclusions:
1)- REAL values are passed in to procedures as 64-bit numbers
2)- REAL values can be returned from procedures as 80-bit numbers
3)- Comparisons (eg a + b >= 1.) operate on 80-bit numbers.
Because the comparison is 80-bit I think the code:
IF a + b <= 1 THEN c := Sqrt (1. - (a+b)) END
is indeed safe.
But what is to stop another compiler from rounding both sides of a + b <= 1. to 64-bit, thereby breaking this code? The above behaviour is not defined in the language.
Gérard wrote
I'd like to think that, whenever a <= 1, then 1 - a >= 0.
I think this is correct if the comparison and subtraction both operate to the same precision. I am confident that the subtraction operates at 80-bits, so this implication requires that the comparison also operates at 80-bits.
Re: Numerical Analysis
Posted: Wed Nov 27, 2013 12:04 pm
by cfbsoftware
But what is to stop another compiler from rounding both sides of a + b <= 1. to 64-bit, thereby breaking this code?
Nothing as far as I know.
Intel have written an article which discusses how they tackle these sorts of issues in their compilers:
http://software.intel.com/en-us/article ... l-compiler
There is a download link to the complete article on that page or you can get it from here:
http://software.intel.com/sites/default ... 2712_1.pdf
Re: Numerical Analysis
Posted: Tue Feb 11, 2014 8:37 am
by Romiras
Rounding errors are the issue? See
The Floating-Point Guide