[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Floating point problems
Hello!
I have a few problems with floating points!
I include three test modules which shall show the
problems I have. The occur in similar ways in
different places, too.
I can't find any rules in why the calculations
have errornous results. (E.g. in a different module
with slightly different MatrixMultPoint I made it
work by introducing output (Err) at the beginning of the
prcedure.)
Testprogramme 1:
----------------
MODULE MiniTestReal1;
IMPORT
L := LowReal,
Err;
TYPE
Matrix * = ARRAY 3,3 OF REAL;
Point * = RECORD
x *, y *, w * : REAL;
END;
VAR
m : Matrix;
p : Point;
PROCEDURE PrintMatrix * ( m : Matrix );
CONST v=8; n =12;
VAR i, j : SHORTINT;
BEGIN
Err.String ("M = ");
FOR i:=0 TO 2 DO
FOR j:=0 TO 2 DO
Err.Real ( m[i,j],v,n );
Err.String (", ");
END;
IF i#2 THEN
Err.Ln;
Err.String (" ");
ELSE
Err.Ln;
END;
END;
END PrintMatrix;
PROCEDURE (VAR p : Point) Print* ;
CONST v=8; n =12;
BEGIN
Err.String ("p.x: "); Err.Real (p.x,v,n); Err.Ln;
Err.String ("p.y: "); Err.Real (p.y,v,n); Err.Ln;
Err.String ("p.w: "); Err.Real (p.w,v,n); Err.Ln;
END Print;
PROCEDURE MatrixMultPoint * ( VAR res : Point; m : Matrix; p : Point );
BEGIN
res.x := m[0,0]*p.x;
res.x := res.x + m[0,1]*p.y;
res.x := res.x + m[0,2]*p.w;
res.y := m[1,0]*p.x;
res.y := res.y + m[1,1]*p.y;
res.y := res.y + m[1,2]*p.w;
res.w := m[2,0]*p.x;
res.w := res.w + m[2,1]*p.y;
res.w := res.w + m[2,2]*p.w;
IF L.IsNaN(res.x) THEN
Err.String( "x is NaN" ); Err.Ln;
PrintMatrix(m);p.Print; Err.Ln;
Err.String("m[0,0]*p.x "); Err.Real(m[0,0]*p.x,4,8); Err.Ln;
IF m[0,0]*p.x=0.0 THEN
Err.String ("m[0,0]*p.x=0"); Err.Ln;
END;
Err.String("+m[0,1]*p.y "); Err.Real(m[0,1]*p.y,4,8); Err.Ln;
IF m[0,1]*p.y=0.0 THEN
Err.String ("m[0,1]*p.y=0"); Err.Ln;
END;
Err.String("+m[0,2]*p.w "); Err.Real(m[0,2]*p.w,4,8); Err.Ln;
IF m[0,2]*p.w=0.0 THEN
Err.String ("m[0,1]*p.w=0"); Err.Ln;
END;
END;
IF L.IsNaN(res.y) THEN
Err.String( "y is NaN" ); Err.Ln;
PrintMatrix(m);p.Print; Err.Ln;
Err.String("m[1,0]*p.x "); Err.Real(m[1,0]*p.x,4,8); Err.Ln;
Err.String("+m[1,1]*p.y "); Err.Real(m[1,1]*p.y,4,8); Err.Ln;
Err.String("+m[1,2]*p.w "); Err.Real(m[1,2]*p.w,4,8); Err.Ln;
END;
IF L.IsNaN(res.w) THEN
Err.String( "w is NaN" ); Err.Ln;
PrintMatrix(m);p.Print; Err.Ln;
Err.String("m[2,0]*p.x "); Err.Real(m[2,0]*p.x,4,8); Err.Ln;
Err.String("+m[2,1]*p.y "); Err.Real(m[2,1]*p.y,4,8); Err.Ln;
Err.String("+m[2,2]*p.w "); Err.Real(m[2,2]*p.w,4,8); Err.Ln;
END;
END MatrixMultPoint;
BEGIN
m [0,0]:=3; m [0,1]:=0; m [0,2]:=0;
m [1,0]:=0; m [1,1]:=3; m [1,2]:=0;
m [2,0]:=0; m [2,1]:=0; m [2,2]:=1;
p.x := 35;
p.y := 5;
p.w := 1;
(* Err.String ("multiply"); Err.Ln;*)
MatrixMultPoint ( p, m, p );
(* p.Print;*)
END MiniTestReal1.
------------
Execution of programme 1:
> MiniTestReal1
>
(No output!)
Testprogramme 2:
----------------
MODULE MiniTestReal2;
IMPORT
L := LowReal,
Err;
TYPE
Matrix * = ARRAY 3,3 OF REAL;
Point * = RECORD
x *, y *, w * : REAL;
END;
VAR
m : Matrix;
p : Point;
PROCEDURE PrintMatrix * ( m : Matrix );
CONST v=8; n =12;
VAR i, j : SHORTINT;
BEGIN
Err.String ("M = ");
FOR i:=0 TO 2 DO
FOR j:=0 TO 2 DO
Err.Real ( m[i,j],v,n );
Err.String (", ");
END;
IF i#2 THEN
Err.Ln;
Err.String (" ");
ELSE
Err.Ln;
END;
END;
END PrintMatrix;
PROCEDURE (VAR p : Point) Print* ;
CONST v=8; n =12;
BEGIN
Err.String ("p.x: "); Err.Real (p.x,v,n); Err.Ln;
Err.String ("p.y: "); Err.Real (p.y,v,n); Err.Ln;
Err.String ("p.w: "); Err.Real (p.w,v,n); Err.Ln;
END Print;
PROCEDURE MatrixMultPoint * ( VAR res : Point; m : Matrix; p : Point );
BEGIN
res.x := m[0,0]*p.x;
res.x := res.x + m[0,1]*p.y;
res.x := res.x + m[0,2]*p.w;
res.y := m[1,0]*p.x;
res.y := res.y + m[1,1]*p.y;
res.y := res.y + m[1,2]*p.w;
res.w := m[2,0]*p.x;
res.w := res.w + m[2,1]*p.y;
res.w := res.w + m[2,2]*p.w;
IF L.IsNaN(res.x) THEN
Err.String( "x is NaN" ); Err.Ln;
PrintMatrix(m);p.Print; Err.Ln;
Err.String("m[0,0]*p.x "); Err.Real(m[0,0]*p.x,4,8); Err.Ln;
IF m[0,0]*p.x=0.0 THEN
Err.String ("m[0,0]*p.x=0"); Err.Ln;
END;
Err.String("+m[0,1]*p.y "); Err.Real(m[0,1]*p.y,4,8); Err.Ln;
IF m[0,1]*p.y=0.0 THEN
Err.String ("m[0,1]*p.y=0"); Err.Ln;
END;
Err.String("+m[0,2]*p.w "); Err.Real(m[0,2]*p.w,4,8); Err.Ln;
IF m[0,2]*p.w=0.0 THEN
Err.String ("m[0,1]*p.w=0"); Err.Ln;
END;
END;
IF L.IsNaN(res.y) THEN
Err.String( "y is NaN" ); Err.Ln;
PrintMatrix(m);p.Print; Err.Ln;
Err.String("m[1,0]*p.x "); Err.Real(m[1,0]*p.x,4,8); Err.Ln;
Err.String("+m[1,1]*p.y "); Err.Real(m[1,1]*p.y,4,8); Err.Ln;
Err.String("+m[1,2]*p.w "); Err.Real(m[1,2]*p.w,4,8); Err.Ln;
END;
IF L.IsNaN(res.w) THEN
Err.String( "w is NaN" ); Err.Ln;
PrintMatrix(m);p.Print; Err.Ln;
Err.String("m[2,0]*p.x "); Err.Real(m[2,0]*p.x,4,8); Err.Ln;
Err.String("+m[2,1]*p.y "); Err.Real(m[2,1]*p.y,4,8); Err.Ln;
Err.String("+m[2,2]*p.w "); Err.Real(m[2,2]*p.w,4,8); Err.Ln;
END;
END MatrixMultPoint;
BEGIN
m [0,0]:=3; m [0,1]:=0; m [0,2]:=0;
m [1,0]:=0; m [1,1]:=3; m [1,2]:=0;
m [2,0]:=0; m [2,1]:=0; m [2,2]:=1;
p.x := 35;
p.y := 5;
p.w := 1;
(* Err.String ("multiply"); Err.Ln;*)
MatrixMultPoint ( p, m, p );
p.Print;
END MiniTestReal2.
-------------
Execution of programme 2:
> MiniTestReal2
p.x: 1.05000000000E+2
p.y: 9.85537172704E-13
p.w: 3.28647590912E-13
>
(Shows the results of matrixmultiplication.)
Testprogramme 3:
----------------
MODULE MiniTestReal3;
IMPORT
L := LowReal,
Err;
TYPE
Matrix * = ARRAY 3,3 OF REAL;
Point * = RECORD
x *, y *, w * : REAL;
END;
VAR
m : Matrix;
p : Point;
PROCEDURE PrintMatrix * ( m : Matrix );
CONST v=8; n =12;
VAR i, j : SHORTINT;
BEGIN
Err.String ("M = ");
FOR i:=0 TO 2 DO
FOR j:=0 TO 2 DO
Err.Real ( m[i,j],v,n );
Err.String (", ");
END;
IF i#2 THEN
Err.Ln;
Err.String (" ");
ELSE
Err.Ln;
END;
END;
END PrintMatrix;
PROCEDURE (VAR p : Point) Print* ;
CONST v=8; n =12;
BEGIN
Err.String ("p.x: "); Err.Real (p.x,v,n); Err.Ln;
Err.String ("p.y: "); Err.Real (p.y,v,n); Err.Ln;
Err.String ("p.w: "); Err.Real (p.w,v,n); Err.Ln;
END Print;
PROCEDURE MatrixMultPoint * ( VAR res : Point; m : Matrix; p : Point );
BEGIN
res.x := m[0,0]*p.x;
res.x := res.x + m[0,1]*p.y;
res.x := res.x + m[0,2]*p.w;
res.y := m[1,0]*p.x;
res.y := res.y + m[1,1]*p.y;
res.y := res.y + m[1,2]*p.w;
res.w := m[2,0]*p.x;
res.w := res.w + m[2,1]*p.y;
res.w := res.w + m[2,2]*p.w;
IF L.IsNaN(res.x) THEN
Err.String( "x is NaN" ); Err.Ln;
PrintMatrix(m);p.Print; Err.Ln;
Err.String("m[0,0]*p.x "); Err.Real(m[0,0]*p.x,4,8); Err.Ln;
IF m[0,0]*p.x=0.0 THEN
Err.String ("m[0,0]*p.x=0"); Err.Ln;
END;
Err.String("+m[0,1]*p.y "); Err.Real(m[0,1]*p.y,4,8); Err.Ln;
IF m[0,1]*p.y=0.0 THEN
Err.String ("m[0,1]*p.y=0"); Err.Ln;
END;
Err.String("+m[0,2]*p.w "); Err.Real(m[0,2]*p.w,4,8); Err.Ln;
IF m[0,2]*p.w=0.0 THEN
Err.String ("m[0,1]*p.w=0"); Err.Ln;
END;
END;
IF L.IsNaN(res.y) THEN
Err.String( "y is NaN" ); Err.Ln;
PrintMatrix(m);p.Print; Err.Ln;
Err.String("m[1,0]*p.x "); Err.Real(m[1,0]*p.x,4,8); Err.Ln;
Err.String("+m[1,1]*p.y "); Err.Real(m[1,1]*p.y,4,8); Err.Ln;
Err.String("+m[1,2]*p.w "); Err.Real(m[1,2]*p.w,4,8); Err.Ln;
END;
IF L.IsNaN(res.w) THEN
Err.String( "w is NaN" ); Err.Ln;
PrintMatrix(m);p.Print; Err.Ln;
Err.String("m[2,0]*p.x "); Err.Real(m[2,0]*p.x,4,8); Err.Ln;
Err.String("+m[2,1]*p.y "); Err.Real(m[2,1]*p.y,4,8); Err.Ln;
Err.String("+m[2,2]*p.w "); Err.Real(m[2,2]*p.w,4,8); Err.Ln;
END;
END MatrixMultPoint;
BEGIN
m [0,0]:=3; m [0,1]:=0; m [0,2]:=0;
m [1,0]:=0; m [1,1]:=3; m [1,2]:=0;
m [2,0]:=0; m [2,1]:=0; m [2,2]:=1;
p.x := 35;
p.y := 5;
p.w := 1;
Err.String ("multiply"); Err.Ln;
MatrixMultPoint ( p, m, p );
p.Print;
END MiniTestReal3.
------------
Execution of programme 3:
> MiniTestReal3
multiply
x is NaN
M = 3.00000000000, 0.00000000000, 0.00000000000,
0.00000000000, 3.00000000000, 0.00000000000,
0.00000000000, 0.00000000000, 1.00000000000,
p.x: 3.50000000000E+1
p.y: 5.00000000000
p.w: 1.00000000000
m[0,0]*p.x 1.0500000E+2
+m[0,1]*p.y 0.0000000
m[0,1]*p.y=0
+m[0,2]*p.w 0.0000000
m[0,1]*p.w=0
y is NaN
M = 3.00000000000, 0.00000000000, 0.00000000000,
0.00000000000, 3.00000000000, 0.00000000000,
0.00000000000, 0.00000000000, 1.00000000000,
p.x: 3.50000000000E+1
p.y: 5.00000000000
p.w: 1.00000000000
m[1,0]*p.x 0.0000000
+m[1,1]*p.y 1.5000000E+1
+m[1,2]*p.w 0.0000000
w is NaN
M = 3.00000000000, 0.00000000000, 0.00000000000,
0.00000000000, 3.00000000000, 0.00000000000,
0.00000000000, 0.00000000000, 1.00000000000,
p.x: 3.50000000000E+1
p.y: 5.00000000000
p.w: 1.00000000000
m[2,0]*p.x 0.0000000
+m[2,1]*p.y 0.0000000
+m[2,2]*p.w 1.0000000
p.x: NaN
p.y: NaN
p.w: NaN
>
Now the matrixmultiplication fails.
The differences in the testprogrammes are just in the last
lines of the module body. The first calls MatrixMultPoint
after initializing the matrix and point, the second prints out
the result after doing the same and the last should only print
"multiply" before doing the multiplication.
I use "OOC/ANSI-C 1.4.1 for Unix" and
"egcs-2.90.25 980302 (egcs-1.0.2 prerelease)".
Can anyone reproduce the results?