Recent

Author Topic: Strange error with arccos()  (Read 5418 times)

Paolo

  • Sr. Member
  • ****
  • Posts: 499
Re: Strange error with arccos()
« Reply #45 on: August 09, 2022, 02:47:38 pm »
Done some testing, some improvements are observed even if very little.

Josh

  • Hero Member
  • *****
  • Posts: 1271
Re: Strange error with arccos()
« Reply #46 on: August 09, 2022, 03:23:08 pm »
Hi

I think it needs to be fixed, as the adding pi and subtracting 2pi, is causing a round error.


compact code; which is equally as fast but highly doubt if would be accepted
Code: Pascal  [Select][+][-]
  1. function JoshArcTan2(Y,X:Float):Float;
  2. function SignGreaterEqualZero(Const AValue:Float):Integer;inline;
  3. begin
  4.   If AValue<0.0 Then SignGreaterEqualZero:=-1 Else SignGreaterEqualZero:=+1;
  5. end;
  6. begin
  7.   If X>0.0 Then JoshArcTan2:=Arctan(Y/X)
  8.   Else
  9.     If X<0.0 then JoshArcTan2:=Arctan(Y/X)+SignGreaterEqualZero(Y)*pi
  10.     else JoshArcTan2:=Sign(Y)*pi/2;
  11. end;      
  12.  

or just correct existing code to stop the rounding issue

Code: Pascal  [Select][+][-]
  1. function arctan2(y,x : float) : float;
  2. begin
  3.   if (x=0) then
  4.     begin
  5.       if y=0 then
  6.         arctan2:=0.0
  7.       else if y>0 then
  8.         arctan2:=pi/2
  9.       else if y<0 then
  10.         arctan2:=-pi/2;
  11.     end
  12.   else
  13.   if x<0.0 then
  14.     begin
  15.       if y<0.0 then
  16.         arctan2:=arctan(y/x)-pi
  17.       else
  18.         arctan2:= arctan(y/x)+pi;
  19.     end
  20.   else arctan2:=ArcTan(y/x);
  21. end;                    
  22.  

I started investigating FPC tanh function, on initial test it is flawed, as the range check is not correct,
try calling it with a value of -356, you get a FT error, according to the range check it is in range..
Not too sure about the formula used, if mem serves it should be using cosh,sinh for calculation, but there is probably other equations to solve it, I just have a gut feeling something is wrong; but at the moment i do not have a good time slot to check it out.
« Last Edit: August 09, 2022, 03:27:29 pm by Josh »
The best way to get accurate information on the forum is to post something wrong and wait for corrections.

Paolo

  • Sr. Member
  • ****
  • Posts: 499
Re: Strange error with arccos()
« Reply #47 on: August 09, 2022, 03:43:18 pm »
I'll submit the second version, for the same reasons you exposed. I'll check tanh !
« Last Edit: August 09, 2022, 03:47:21 pm by Paolo »

Paolo

  • Sr. Member
  • ****
  • Posts: 499
Re: Strange error with arccos()
« Reply #48 on: August 09, 2022, 04:43:20 pm »
for tanh the suspect is:

Code: Pascal  [Select][+][-]
  1. x:=-359;
  2. ...
  3. y:=sinh(x)/cosh(x);   //works as expected
  4. y:=(Exp(x)-Exp(-x))/(Exp(X)+Exp(-x); //works as above
  5. y:=Tanh(x); //crash !
  6. y:=(1-exp(-2x))/(1+Exp(-2*x)); //that is could be as implemented tanh and crash again,
  7. //and you see way : -2*X is > 0 and Exp itself crash ! too large value
  8.  
  9.  

if the hypotesis above are right and the computation of tanh over extended range
should be (using the tanh(-x)=-tanh(x), x > 0)

Code: Pascal  [Select][+][-]
  1. funcyion NewTanh(x:double):double;
  2. if X < 0 then
  3.   result:=-(1-exp(2x))/(1+Exp(2*x)); //ie -tanh(-x) (here X < 0, tanh(x)=-tanh(-x)
  4. else
  5.   result:=(1-exp(-2x))/(1+Exp(-2*x)); //ie tanh(x)
  6.  

now working for x:=-359;

do you agree ?
« Last Edit: August 09, 2022, 04:49:47 pm by Paolo »

Paolo

  • Sr. Member
  • ****
  • Posts: 499
Re: Strange error with arccos()
« Reply #49 on: August 10, 2022, 12:42:41 pm »
but in any case Tanh(x) for x approaching 18 already gives 1, so tanh with X above such values is = 1, non need to really compute it.

Josh

  • Hero Member
  • *****
  • Posts: 1271
Re: Strange error with arccos()
« Reply #50 on: August 10, 2022, 02:18:58 pm »
hi

i agree, when going to 50 digits, the min/max is only 60 before it is +-1/

below code that test varius routine, note joshtanh is fast and just as accurate and going in the range +-60(so has an advantage)

Code: Pascal  [Select][+][-]
  1. program Project1;
  2. uses sysutils,math,crt;
  3. const n=20000000;
  4.       numberofroutines=4;
  5.  
  6. Type TTanHFunc = function(const x: Float): Float;
  7.  
  8. var i:longint;
  9.     v:float;
  10.     SampleData,Differences:array[1..N] of Float;
  11.     results:array[1..n,1..numberofroutines] of Float;
  12.     verify:boolean=true;
  13.     BaseTime:double;
  14.     TimeInRoutine:Array[1..numberofroutines] of double;
  15.     RoutineNames:Array[1..numberofroutines] of String;
  16.     mc:integer=0;
  17.  
  18. const  MaxTanh = 5678.22249441322; // Ln(MaxExtended)/2
  19.        TestMaxRange=350;
  20.  
  21. function mathtanh(const x : float) : float;
  22. var Temp : float;
  23. begin
  24.   if x>MaxTanh then exit(1.0)
  25.   else if x<-MaxTanh then exit (-1.0);
  26.   temp:=exp(-2*x);
  27.   mathtanh:=(1-temp)/(1+temp)
  28. end;
  29.  
  30.  
  31. function NewTanhIdea1(const x:float):float;
  32. var SignX:Integer;
  33. begin
  34.   SignX:=sign(x);
  35.   if SignX=0 then result:=0
  36.   else result:=SignX*(
  37.                      (1-exp(((SignX*-1)*2)*x))
  38.                           /
  39.                      (1+exp(((SignX*-1)*2)*x)));
  40. end;
  41.  
  42. function NewTanh2(const x:float):float;
  43. var SignX:Integer;
  44. begin
  45.   if X < 0 then
  46.   result:=-(1-exp(2*x))/(1+Exp(2*x)) //ie -tanh(-x) (here X < 0, tanh(x)=-tanh(-x)
  47. else
  48.   result:=(1-exp(-2*x))/(1+Exp(-2*x)); //ie tanh(x)
  49. end;
  50.  
  51.  
  52.  
  53.  
  54. function JoshNewTanh(const X: float): float;
  55. const xrange=60;
  56. begin
  57.   if X>xrange then Result := 1.0
  58.   else
  59.     if X<-xrange then Result := -1.0
  60.     else
  61.       begin
  62.         Result := Exp(X);
  63.         Result := Result*Result;
  64.         Result := (Result-1.0)/(Result+1.0);
  65.       end;
  66. end;
  67.  
  68. var t:tdatetime;
  69.     d:float;
  70.     tol:float=0.00000000000001;
  71.     dx:integer=48;
  72.     dy:integer=1;
  73.  
  74. procedure PrepareData;
  75. var j,k:longint;
  76. begin
  77.   WriteLn('Preparing Dataset');
  78.   for j:=1 to n do
  79.   begin
  80.     SampleData[j]:=random-0.5*(2*TestMaxRange);
  81.     for k:=1 to numberofroutines do results[j,k]:=0;
  82.   end;
  83.   for k:=1 to numberofroutines do differences[k]:=0;
  84. end;
  85.  
  86. Function Test(Descr: String; TanHFunc: TTanHFunc; r:integer):Double;
  87. var
  88.   t: TDateTime;
  89.   i: Integer;
  90.   ans:float;
  91. begin
  92.   Write('  Testing ' + Descr + '...');
  93.   t:=Now;
  94.   for i := 1 to N do
  95.   begin
  96.     ans:=TanHFunc(SampleData[i]);
  97.     results[i,r]:=ans;
  98.     if verify then if abs(results[i,r]-results[i,1])>tol then differences[r]:=differences[r]+abs(results[i,r]-results[i,1]);
  99.   end;
  100.   t := Now - t;
  101.   if r=1 then BaseTime:=t;
  102.   WriteLn(' done: ', FormatDateTime('s.zzz', t), ' sec');
  103.   if verify then writeLn('     Differences :',differences[r]:8:8);
  104.   WriteLn;
  105.   TimeInRoutine[r]:=TimeInRoutine[r]+t;
  106.   result:=t;
  107. end;
  108.  
  109.  
  110.  
  111. begin
  112.   RoutineNames[1]:='MathTanH....';
  113.   RoutineNames[2]:='NewTanh2....';
  114.   RoutineNames[3]:='NewTanhIdea1';
  115.   RoutineNames[4]:='JoshNewTanh.';
  116.   randomize;
  117.   randseed:=0;
  118.   mc:=0;
  119.   repeat
  120.     inc(mc);
  121.     PrepareData;
  122.     t:=test(RoutineNames[1],@mathtanh,1);
  123.     t:=test(RoutineNames[2],@NewTanh2,2);
  124.     t:=test(RoutineNames[3],@NewTanhIdea1,3);
  125.     t:=test(RoutineNames[4],@JoshNewTanh,4);
  126.     TextBackground(blue);
  127.     gotoxy(dx,dy);write('        Total Times   Count '+mc.ToString); clreol;
  128.     gotoxy(dx,dy+1);write('        hr mn se ms     Avg');  clreol;
  129.     if mc>0 then
  130.     begin
  131.       for i:=1 to numberofroutines do
  132.       begin;
  133.         gotoxy(dx,dy+1+i);
  134.         write(' '+RoutineNames[i]+':', FormatDateTime('hh.nn.ss.zzz', TimeInRoutine[i]), '  ',FormatDateTime('s.zzz',TimeInRoutine[i]/mc), ' sec'); clreol;
  135.       end;
  136.     end;
  137.     TextBackground(black);
  138.     gotoxy(1,1);
  139.   until false;
  140. end.
  141.  

What u think?
« Last Edit: August 10, 2022, 03:36:37 pm by Josh »
The best way to get accurate information on the forum is to post something wrong and wait for corrections.

wp

  • Hero Member
  • *****
  • Posts: 11853
Re: Strange error with arccos()
« Reply #51 on: August 10, 2022, 03:49:08 pm »
Code: Pascal  [Select][+][-]
  1. procedure PrepareData;
  2. [...]
  3.     SampleData[j]:=random-0.5*(2*TestMaxRange);
  4.  
What's this? Why don't you cancel the 2? Or do you mean this, with brackets around "random - 0.5"?
Code: Pascal  [Select][+][-]
  1.     SampleData[j]:=(random-0.5)*(2*TestMaxRange);
  2.  

Paolo

  • Sr. Member
  • ****
  • Posts: 499
Re: Strange error with arccos()
« Reply #52 on: August 10, 2022, 04:07:00 pm »
I like newtanh2, because in any case it computes exp function in the right direction (ie toward 0) so no overflow can occur and it should work for any x without check input value range
« Last Edit: August 10, 2022, 05:08:45 pm by Paolo »

Josh

  • Hero Member
  • *****
  • Posts: 1271
Re: Strange error with arccos()
« Reply #53 on: August 10, 2022, 05:43:06 pm »
@wp
yes it should be
Code: [Select]
SampleData[j]:=(random-0.5)*(2*TestMaxRange)Was not thinking straight, initially i was checking in small range of values to verify results across routines. then later expanded to larger values.

@Paolo
I like the newtan2 function because it does not burp on large numbers,

The fact the the Math routine checks for a +- range of 5678.22249441322, which itself cannot handle and causes a crash is a problem and needs to be fixed, a range check of +-26 should be adequate for a precision of 22 decimal places, which i think is the accuracy of the extended type.

The best way to get accurate information on the forum is to post something wrong and wait for corrections.

Paolo

  • Sr. Member
  • ****
  • Posts: 499
Re: Strange error with arccos()
« Reply #54 on: August 10, 2022, 06:21:56 pm »
My suspect is the fact is that the code below

Quote
Const MaxTanh = 5678.22249441322; // Ln(MaxExtended)/2

function tanh(x : float) : float;
  var Temp : float;
  begin
     if x>MaxTanh then
       exit(1.0)
     else
         if x<-MaxTanh then
           exit (-1.0);
     temp:=exp(-2*x);
     tanh:=(1-temp)/(1+temp)
  end;

is a problem for:
1- it does not consider correctly the large negaive numbers (as we aleady discussed)
2- as in the comment "//Ln(MaxExtended)/2" for win 64 there is not "Extended" (it is an alias of double), so the max imput is wrongly assumed, not checked ma it should be

const MAXTANH = somevalue;  //ln(Maxdouble)/2, something just less 365

Paolo

  • Sr. Member
  • ****
  • Posts: 499
Re: Strange error with arccos()
« Reply #55 on: August 10, 2022, 06:40:02 pm »

wp

  • Hero Member
  • *****
  • Posts: 11853
Re: Strange error with arccos()
« Reply #56 on: August 10, 2022, 06:47:27 pm »
for the arctan2

https://gitlab.com/freepascal.org/fpc/source/-/issues/39861
and fixed 20 minutes later - that's the power of Open Source.

Paolo

  • Sr. Member
  • ****
  • Posts: 499
Re: Strange error with arccos()
« Reply #57 on: August 10, 2022, 07:00:19 pm »
@wp fixed even before I copy paste here the link, unbeliveable  :o!!!!! faster than the arctan  :D

Paolo

  • Sr. Member
  • ****
  • Posts: 499
Re: Strange error with arccos()
« Reply #58 on: August 18, 2022, 12:06:46 am »
For who is interested on tanh issue topic

https://gitlab.com/freepascal.org/fpc/source/-/issues/39867

 

TinyPortal © 2005-2018