Lazarus

Programming => General => Topic started by: taazz on November 15, 2017, 07:00:08 am

Title: Cos(Theta). wrong results?
Post by: taazz on November 15, 2017, 07:00:08 am
In the process of trying a matrix library (DynMatrix by 5dpo (https://paginas.fe.up.pt/~paco/pmwiki/index.php?n=DynMatrix.DynMatrix)) I'm a bit confused with the results of COS they seem, wrong.
Here is a small test app.
Code: Pascal  [Select][+][-]
  1. program project1;
  2.  
  3. {$mode delphi}{$H+}
  4. {$APPTYPE CONSOLE}
  5. uses
  6.   {$IFDEF UNIX}{$IFDEF UseCThreads}
  7.   cthreads,
  8.   {$ENDIF}{$ENDIF}
  9.   Classes, math
  10.   { you can add units after this };
  11.  
  12. {$R *.res}
  13. var
  14.   vRad, vCos, vSin, vRCos, vRSin : float;
  15.  
  16. begin
  17.   vRad := degtorad(90);
  18.   vCos:= cos(90);
  19.   vRCos:= cos(vRad);
  20.   vSin := sin(90);
  21.   vRSin := sin(vRad);
  22.   WriteLn('Deg/Rad', 90,':',vRad);
  23.   WriteLn('Cos ', vCos ,':',vRCos);
  24.   WriteLn('Sin ', vSin ,':',vRSin);
  25.   ReadLn;
  26. end.
  27.  
and the result from lazarus 1.6.4 is
Quote
Deg/Rad90 : 1.57079632679489661926E+0000
Cos -4.48073616129170152476E-0001:-2.71050543121376108502E-0020
Sin   8.93996663600557890477E-0001: 1.00000000000000000000E+0000
I was under the impression that cos of 90 degrees is 0, sin at least seems to be correct. I'm almost certain that I miss something crucial but I spend 4 hours looking and I can't see it.
Title: Re: Cos(Theta). wrong results?
Post by: Eugene Loza on November 15, 2017, 07:10:26 am
-2.71050543121376108502E-0020 seems close enough to zero, isn't it? I.e. that's all the accuracy digits the compiler actually has, and it isn't even able to determine Pi/2 with greater accuracy without special effort.
Title: Re: Cos(Theta). wrong results?
Post by: Thaddy on November 15, 2017, 08:54:25 am
That's correct. Note that the math unit has IsZero functions for the case you need to determine if a float value should be considered equal to zero:
Code: Pascal  [Select][+][-]
  1. function IsZero(const A: Single; Epsilon: Single): Boolean; overload;
  2. function IsZero(const A: Single): Boolean;inline; overload;
  3. {$ifdef FPC_HAS_TYPE_DOUBLE}
  4. function IsZero(const A: Double; Epsilon: Double): Boolean; overload;
  5. function IsZero(const A: Double): Boolean;inline; overload;
  6. {$endif FPC_HAS_TYPE_DOUBLE}
  7. {$ifdef FPC_HAS_TYPE_EXTENDED}
  8. function IsZero(const A: Extended; Epsilon: Extended): Boolean; overload;
  9. function IsZero(const A: Extended): Boolean;inline; overload;
  10. {$endif FPC_HAS_TYPE_EXTENDED}
  11.  

Which gives using your slightly modified code:
Code: [Select]
Deg/Rad90: 1.5707963267948966E+000 FALSE
Cos -4.4807361612917013E-001: 6.1232339957367660E-017  TRUE
Sin  8.9399666360055796E-001: 1.0000000000000000E+000 FALSE

using:
Code: Pascal  [Select][+][-]
  1. program test_freepascal;
  2. {$mode delphi}{$H+}
  3. {$APPTYPE CONSOLE}
  4. uses
  5.   {$IFDEF UNIX}{$IFDEF UseCThreads}
  6.   cthreads,
  7.   {$ENDIF}{$ENDIF}
  8.   Classes, math
  9.   { you can add units after this };
  10.  
  11. var
  12.   vRad, vCos, vSin, vRCos, vRSin : float;
  13.  
  14. begin
  15.   vRad := degtorad(90);
  16.   vCos:= cos(90);
  17.   vRCos:= cos(vRad);
  18.   vSin := sin(90);
  19.   vRSin := sin(vRad);
  20.   WriteLn('Deg/Rad', 90,':',vRad,Iszero(vRad):6);
  21.   WriteLn('Cos ', vCos ,':',vRCos,IsZero(vRCos):6);
  22.   WriteLn('Sin ', vSin ,':',vRSin,IsZero(vRSin):6);
  23.   ReadLn;
  24. end.
Title: Re: Cos(Theta). wrong results?
Post by: munair on November 15, 2017, 09:03:44 am
E-20? How accurate do you want it to be? It's the same old story. Binary systems cannot handle floating point numbers very well. It is always a close approximation.

A simple comparison: my good old calculator already rounds off DegToRad (which is simply x PI/180) to 1.5707963267. In Math a rule of thumb is to round off using the value with the least digits within your calculations. So first you should determine how much accuracy you want and then apply that to your calculations. Accurate Float is up to 18 digits if I'm not mistaken so E-20 would then render 0. Beyond that, results become unpredictable because the computer returns the closest approximation that can be expressed as powers of 2.
Title: Re: Cos(Theta). wrong results?
Post by: Thaddy on November 15, 2017, 09:05:55 am
@Munair: look at my answer, they crossed. There is no rule of thumb, there is a mathematical convention to test for zero.
Title: Re: Cos(Theta). wrong results?
Post by: munair on November 15, 2017, 09:10:15 am
@Munair: look at my answer, they crossed. There is no rule of thumb, there is a mathematical convention to test for zero.
Yes, but with a certain precision in mind.
In my program Solpage that computes planetary angles, I had the same problem. So I established a specific precision I wanted to work with and rounded the angles accordingly:
Code: Pascal  [Select][+][-]
  1. L1 := RoundTo(a.b2.Coord.OSC.L, -kHighPrecision);
In this case an angle between 0 and 180 degrees is returned with 8 digits precision. Anything beyond that would be rediculous.
Title: Re: Cos(Theta). wrong results?
Post by: munair on November 15, 2017, 09:17:50 am
Actually, my example computes the longitude of two planets. These can be >= 0 and < 360 degrees. Then the smaller angle (>= 0 and < 180) is returned:
Code: Pascal  [Select][+][-]
  1. L1 := RoundTo(a.b2.Coord.OSC.L, -kHighPrecision);
  2. L2 := RoundTo(a.b3.Coord.OSC.L, -kHighPrecision);
  3. a.Size := SmallAngle(Abs(L1 - L2));
Without rounding to a specific precision, the engine would either become terribly slow, or inaccurate preventing me to see when both planets reach the same longitude (a.Size = 0.00000000).  :D
Title: Re: Cos(Theta). wrong results?
Post by: Thaddy on November 15, 2017, 09:46:15 am
Then you should also look at the SameValue functions in math, etc...... You are doing double work...
Title: Re: Cos(Theta). wrong results?
Post by: munair on November 15, 2017, 10:54:42 am
Then you should also look at the SameValue functions in math, etc...... You are doing double work...
No, because any angle (a.Size) must be returned. I just need to make sure that angles are rounded to a specific number of digits so that later on IF a.Size = 0, or IF a.Size = 90 returns true.  ;)
Title: Re: Cos(Theta). wrong results?
Post by: Thaddy on November 15, 2017, 11:41:06 am
Did you check the Epsilon parameters in that set?..... <slightly amused, not yet grumpy.. 8-)>
Title: Re: Cos(Theta). wrong results?
Post by: taazz on November 15, 2017, 04:14:38 pm
-2.71050543121376108502E-0020 seems close enough to zero, isn't it? I.e. that's all the accuracy digits the compiler actually has, and it isn't even able to determine Pi/2 with greater accuracy without special effort.
Well after a good night's sleep I managed to scroll my eyes to the right far enough to see it too, thank you.

E-20? How accurate do you want it to be?
Oh not much, its just for screen coordinates so 2 decimal digits should suffice.  I just did not look far enough for 4 hours!!

That's correct. Note that the math unit has IsZero functions for the case you need to determine if a float value should be considered equal to zero:
Code: Pascal  [Select][+][-]
  1. function IsZero(const A: Single; Epsilon: Single): Boolean; overload;
  2. function IsZero(const A: Single): Boolean;inline; overload;
  3. {$ifdef FPC_HAS_TYPE_DOUBLE}
  4. function IsZero(const A: Double; Epsilon: Double): Boolean; overload;
  5. function IsZero(const A: Double): Boolean;inline; overload;
  6. {$endif FPC_HAS_TYPE_DOUBLE}
  7. {$ifdef FPC_HAS_TYPE_EXTENDED}
  8. function IsZero(const A: Extended; Epsilon: Extended): Boolean; overload;
  9. function IsZero(const A: Extended): Boolean;inline; overload;
  10. {$endif FPC_HAS_TYPE_EXTENDED}
  11.  

Which gives using your slightly modified code:
Code: [Select]
Deg/Rad90: 1.5707963267948966E+000 FALSE
Cos -4.4807361612917013E-001: 6.1232339957367660E-017  TRUE
Sin  8.9399666360055796E-001: 1.0000000000000000E+000 FALSE
Thank you sir, that will prove very useful in the near future.

sorry to waste your time guys and thanks for heads up.
Title: Re: Cos(Theta). wrong results?
Post by: munair on November 15, 2017, 09:43:30 pm
Did you check the Epsilon parameters in that set?..... <slightly amused, not yet grumpy.. 8-)>
Yes I did. Already back then. What you're looking at is part of a general procedure spitting out angles at a given time T. So it is not a comparison right away. Rather, it's the level of accuracy that actually makes the given angles meaningful as a function of time. From there it is much simpler to test (because the difference is already computed) if a.Size is 0 or whatever depending on the request, than having to compare the angles again using SameValue getting either true or false, which would be completely meaningless. ;D
TinyPortal © 2005-2018