Recent

Author Topic: Cos(Theta). wrong results?  (Read 4378 times)

taazz

  • Hero Member
  • *****
  • Posts: 5368
Cos(Theta). wrong results?
« on: November 15, 2017, 07:00:08 am »
In the process of trying a matrix library (DynMatrix by 5dpo) 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.
Good judgement is the result of experience … Experience is the result of bad judgement.

OS : Windows 7 64 bit
Laz: Lazarus 1.4.4 FPC 2.6.4 i386-win32-win32/win64

Eugene Loza

  • Hero Member
  • *****
  • Posts: 663
    • My games in Pascal
Re: Cos(Theta). wrong results?
« Reply #1 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.
My FOSS games in FreePascal&CastleGameEngine: https://decoherence.itch.io/ (Sources: https://gitlab.com/EugeneLoza)

Thaddy

  • Hero Member
  • *****
  • Posts: 14205
  • Probably until I exterminate Putin.
Re: Cos(Theta). wrong results?
« Reply #2 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.
« Last Edit: November 15, 2017, 09:04:23 am by Thaddy »
Specialize a type, not a var.

munair

  • Hero Member
  • *****
  • Posts: 798
  • compiler developer @SharpBASIC
    • SharpBASIC
Re: Cos(Theta). wrong results?
« Reply #3 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.
keep it simple

Thaddy

  • Hero Member
  • *****
  • Posts: 14205
  • Probably until I exterminate Putin.
Re: Cos(Theta). wrong results?
« Reply #4 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.
Specialize a type, not a var.

munair

  • Hero Member
  • *****
  • Posts: 798
  • compiler developer @SharpBASIC
    • SharpBASIC
Re: Cos(Theta). wrong results?
« Reply #5 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.
« Last Edit: November 15, 2017, 09:13:09 am by Munair »
keep it simple

munair

  • Hero Member
  • *****
  • Posts: 798
  • compiler developer @SharpBASIC
    • SharpBASIC
Re: Cos(Theta). wrong results?
« Reply #6 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
« Last Edit: November 15, 2017, 09:33:13 am by Munair »
keep it simple

Thaddy

  • Hero Member
  • *****
  • Posts: 14205
  • Probably until I exterminate Putin.
Re: Cos(Theta). wrong results?
« Reply #7 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...
Specialize a type, not a var.

munair

  • Hero Member
  • *****
  • Posts: 798
  • compiler developer @SharpBASIC
    • SharpBASIC
Re: Cos(Theta). wrong results?
« Reply #8 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.  ;)
« Last Edit: November 15, 2017, 10:56:25 am by Munair »
keep it simple

Thaddy

  • Hero Member
  • *****
  • Posts: 14205
  • Probably until I exterminate Putin.
Re: Cos(Theta). wrong results?
« Reply #9 on: November 15, 2017, 11:41:06 am »
Did you check the Epsilon parameters in that set?..... <slightly amused, not yet grumpy.. 8-)>
Specialize a type, not a var.

taazz

  • Hero Member
  • *****
  • Posts: 5368
Re: Cos(Theta). wrong results?
« Reply #10 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.
Good judgement is the result of experience … Experience is the result of bad judgement.

OS : Windows 7 64 bit
Laz: Lazarus 1.4.4 FPC 2.6.4 i386-win32-win32/win64

munair

  • Hero Member
  • *****
  • Posts: 798
  • compiler developer @SharpBASIC
    • SharpBASIC
Re: Cos(Theta). wrong results?
« Reply #11 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
« Last Edit: November 15, 2017, 10:00:07 pm by Munair »
keep it simple

 

TinyPortal © 2005-2018