[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?