Numerical Analysis

Programming language questions
User avatar
Robert
Posts: 177
Joined: Sat Sep 28, 2013 11:04 am
Location: Edinburgh, Scotland

Numerical Analysis

Post by Robert »

I have a function, say 'Foo', that returns

Code: Select all

Math.Sqrt (1. - u * u - v * v) 
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

Code: Select all

Math.Sqrt (1 - (u * u + v * v)) 
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':

Code: Select all

IF u * u + v * v  <=  1. - 1.E-15
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!
Josef Templ
Posts: 262
Joined: Tue Sep 17, 2013 6:50 am

Re: Numerical Analysis

Post 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
manumart1
Posts: 67
Joined: Tue Sep 17, 2013 6:25 am

Re: Numerical Analysis

Post by manumart1 »

Instead of

Code: Select all

IF u * u + v * v  <=  1. THEN
you can use

Code: Select all

IF (1. - u * u - v * v) >= 0 THEN
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
gerard94
Posts: 1
Joined: Tue Sep 17, 2013 11:43 am

Re: Numerical Analysis

Post 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
User avatar
Robert
Posts: 177
Joined: Sat Sep 28, 2013 11:04 am
Location: Edinburgh, Scotland

Re: Numerical Analysis

Post 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
manumart1
Posts: 67
Joined: Tue Sep 17, 2013 6:25 am

Re: Numerical Analysis

Post 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
Last edited by manumart1 on Wed Dec 04, 2013 3:32 pm, edited 1 time in total.
Rex Couture
Posts: 2
Joined: Tue Nov 19, 2013 6:26 pm

Re: Numerical Analysis

Post 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.
User avatar
Robert
Posts: 177
Joined: Sat Sep 28, 2013 11:04 am
Location: Edinburgh, Scotland

Re: Numerical Analysis

Post 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.
cfbsoftware
Posts: 55
Joined: Wed Sep 18, 2013 10:06 pm
Contact:

Re: Numerical Analysis

Post 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
Romiras
Posts: 35
Joined: Tue Sep 17, 2013 5:55 am
Location: Tel-Aviv

Re: Numerical Analysis

Post by Romiras »

Rounding errors are the issue? See The Floating-Point Guide
Post Reply