Lazarus

Programming => General => Topic started by: linuxbob on August 05, 2022, 02:49:06 am

Title: Strange error with arccos()
Post by: linuxbob on August 05, 2022, 02:49:06 am
I am working on a project that needs to convert between rectangular and polar coordinates. I am sometimes getting an error:

"Project R2P raised exception class 'External: FLT INVALID OPERATION'. If I proceed, the debugger jumps to the assembler instruction window at the instruction 'sqrtsd xmm0, xmm2'.

In the debugger, the error occurs on this FP code:

rx1 := rx * arccos(ry);

where ry is an angular value between -Pi and +Pi radians. The next code line is:

ry1 := rx * arcsin(ry);

which never produces an error.

The intention of this code is, given a vector rx at angle ry, to load rx1 with the X coordinate of the vector and load ry1 with the y coordinate of the vector. I am an electrical engineer by profession, and I use these values to determine the amount of kvars I need to mitigate a poor power factor.

Suggestions?
Title: Re: Strange error with arccos()
Post by: mtrsoft on August 05, 2022, 06:11:06 am
You are trying to take the square root of a negative number.

The algorithym that calculates arccos has a square root in it.

The "principal range" for arccos(y) is 0 <= y <= Pi and for arcsin(y) is -Pi/2 <= y <= Pi/2.

See: https://en.wikipedia.org/wiki/Inverse_trigonometric_functions

Regards,
John
Title: Re: Strange error with arccos()
Post by: Paolo on August 05, 2022, 08:22:22 am
From what you say the input range of ry is invalid for arccos (you said +-pi) Argument for arccos must be |arg|<= 1. I suspect that you are doing something conceptually wrong, argument of arccos should not be an angular value. But please provide some line of code.
Title: Re: Strange error with arccos()
Post by: wp on August 05, 2022, 08:48:51 am
The best function to get the polar angle from cartesian coordinates (x, y) is the function arctan2(y, x) because it also finds the correct quadrant for the angle.
Title: Re: Strange error with arccos()
Post by: Paolo on August 05, 2022, 08:55:08 am
Just to complete the answer. Given (R,T) polar coordinate and (X,Y) rectangular coordinate the relationship is
X=R*cos(T)
Y=R*sin(T)
And
R=SQRT(X^2+Y^2)
T=ARCTAN2(Y,X)  --> arctan2 handle correctly the sign of X and Y
I hope this can help.
Title: Re: Strange error with arccos()
Post by: wp on August 05, 2022, 09:15:23 am
Or, to complete it furthermore, here's a ready-to-use procedure for the conversion:
Code: Pascal  [Select][+][-]
  1. uses
  2.   Math;
  3. procedure CartesianToPolar(x, y: Double; var r, phi: Double);
  4. begin
  5.   if (x = 0) and (y = 0) then
  6.     raise Exception.Create('Undefined polar coordinates for x=0 and y=0');  // Or maybe use some "acceptable" result (r=0, phi=0) to avoid the exception
  7.   r := sqrt(x*x + y*y);
  8.   phi := arctan2(y, x);
  9. end;
Title: Re: Strange error with arccos()
Post by: Paolo on August 05, 2022, 09:37:17 am
@wp, that is mathematically fine.
Infact when X=Y=0 phase is undefined... but for practical purposes once X=Y=0 that is the vector is null, or other phisical parameter does not exist, who cares of phase value ?
so usually the phase is put to 0. And in this sense it is what ARCTAN2 do in its code: if X=Y=0 it says result=0.
So this should be usually enough
Code: Pascal  [Select][+][-]
  1.   r := sqrt(x*x + y*y);
  2.   phi := arctan2(y, x);
  3.  
Title: Re: Strange error with arccos()
Post by: Paolo on August 05, 2022, 09:49:32 am
on the occasion of this discussion I want to underline some aspects of the arctan2 calculus in math units which in my opinion should be modified

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.       ArcTan2:=ArcTan(y/x);
  14.     if x<0.0 then
  15.       ArcTan2:=ArcTan2+pi;
  16.     if ArcTan2>pi then
  17.       ArcTan2:=ArcTan2-2*pi;
  18.   end;  
  19.  

in the code above I see something inefficient

1- there should be a "begin..end" block after else, otherwise the final part is execute even if the result is already determinated

Code: Pascal  [Select][+][-]
  1. function arctan2(y,x : float) : float;
  2.   begin
  3.     if (x=0) then begin
  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. begin <----- begin here, otherwise the two next "if" are executed but if X=0 this is useless !!!
  14.       ArcTan2:=ArcTan(y/x);
  15.     if x<0.0 then
  16.       ArcTan2:=ArcTan2+pi;
  17.     if ArcTan2>pi then
  18.       ArcTan2:=ArcTan2-2*pi;
  19. end; <-----------------------------
  20.   end;  
  21.  

2- why the third check in the first "if" chain, if alerady checked Y=0 and Y>0 what is the meaning of testing furhter if  y<0 ?

Code: Pascal  [Select][+][-]
  1. function arctan2(y,x : float) : float;
  2.   begin
  3.     if (x=0) then begin
  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  <---- remove this
  10.           arctan2:=-pi/2;   <--- just do that
  11.       end
  12.     else
  13. begin <----- here, otherwise the two next if are perfomed but if X=0 this is useless !!!
  14.       ArcTan2:=ArcTan(y/x);
  15.     if x<0.0 then
  16.       ArcTan2:=ArcTan2+pi;
  17.     if ArcTan2>pi then
  18.       ArcTan2:=ArcTan2-2*pi;
  19. end; <----
  20.   end;  
  21.  

3- finally the last part I simplified as follow (final implementation,) why add and subtract twice some values ?

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

any comment on that ?
Title: Re: Strange error with arccos()
Post by: Paolo on August 05, 2022, 09:55:35 am
@wp, I see now the comment in the code you showed

Code: Pascal  [Select][+][-]
  1. // Or maybe use some "acceptable" result (r=0, phi=0) to avoid the exception
  2.  

 :)
Title: Re: Strange error with arccos()
Post by: wp on August 05, 2022, 10:30:24 am
Quote from: Paolo on August 05, 2022, 09:49:32 am
Code: Pascal  [Select][+][-]
  1. function MyArcTan2(Y, X : float) : float;
  2. begin
  3.   if (X = 0) then begin
  4.     if (Y = 0) then
  5.       Result:=0
  6.     else
  7.       if (Y > 0) then
  8.         Result:=0.5*pi
  9.       else
  10.         Result:=-0.5*pi
  11.   end
  12.   else begin
  13.     Result:=ArcTan(Y/X);
  14.     if (X < 0) then
  15.       if (Y < 0) then
  16.         Result:=Result-pi
  17.       else
  18.         Result:=Result+pi
  19.   end;
  20. end;
  21.  
any comment on that ?
Submit it to bugtracker so that i can be considered by the FPC devs.
Title: Re: Strange error with arccos()
Post by: Josh on August 05, 2022, 11:39:36 am
hi

not a math guru, but hopefully my logic flow is working, does this do the same

Code: Pascal  [Select][+][-]
  1. function JoshArcTan2 (y, x : double) : double;
  2. begin
  3.  result:=0.0; // set to 0.0 if it falls through with x and y =0.00 also stops hint result not initialized
  4.  if x>0.0 then result := arctan (y/x)                      // using /2 rather than *.5 as old school / was generally more accurate and for me more readable
  5.    else if x <0.0 then result := arctan (y/x) + pi
  6.      else if y<>0.0 then result := pi/2 * sign (y);   // x=0 so if y is also 0.0 result stays set at 0.00 else result gets set +-pi/2
  7.  
  8.  if result>pi then result:=result-2*pi    // keep result in range of +- pi   could use abs sign here but then would be used every iteration of function
  9.    else if result<-pi then result:=result+2*pi;    // not sure if needed but safguarded like original code
  10. end;    

edited added various comments and range check of +- pi
Title: Re: Strange error with arccos()
Post by: Paolo on August 05, 2022, 01:07:11 pm
here my consideration (see attached image comparing the code and the passes executed for the test cases).
I did some test cases, at the moment my code seems the shorter in execution  :D
In particular if X=Y=0 is quickly evaluated... instead of pass through all the code.

your code seems suffering from point-2 I highlighted in previous post, that is even if some value computation are ended you still go through the code

I attached also a mini project (LAZ 2.2.2/fpc 3.2.2) where some attempt in evaluation arctan routines is done (PS: note that ARCTAN2 of lazarus is executed because I working on Win10 64 bit, other OS have different implementation)

in the code you can put X and Y values as you like and in the code of right button you can select the routine you like to execute.

Obviously all the above if I am not wrong !
Title: Re: Strange error with arccos()
Post by: wp on August 05, 2022, 01:45:28 pm
My 2 cents: this is not worth the effort. Run attached test which populates an array with 20 million pairs of random numbers between -10..10, calculates the arctan2 for them and measures the time. The math.arctan2 function is faster by about 20% than the optimizations presented here because the math unit calls an assembler routine, at least on my system with Win11.

Quote
Testing Math.arctan2... done: 0.816 sec
Testing PaoloArctan2... done: 1.062 sec
Testing JoshArctan2... done: 1.079 sec
Finished. Press ENTER to close.
Title: Re: Strange error with arccos()
Post by: winni on August 05, 2022, 01:58:30 pm
Quote from: wp on August 05, 2022, 01:45:28 pm
The math.arctan2 function is faster by about 20% than the optimizations presented here because the math unit calls an assembler routine, at least on my system with Win11.

Hi!

Also for Linux.
That's all:

Code: Pascal  [Select][+][-]
  1. function arctan2(y,x : float) : float;assembler;
  2.   asm
  3.      fldt y
  4.      fldt x
  5.      fpatan
  6.      fwait
  7.   end;
  8.  

Winni
Title: Re: Strange error with arccos()
Post by: Paolo on August 05, 2022, 02:03:49 pm
@wp, please check on win10 64-bit, I understand that in this case the code called is this one, (no ASM)

Code: Pascal  [Select][+][-]
  1.  
  2. {$ifndef FPC_MATH_HAS_ARCTAN2}
  3. function arctan2(y,x : float) : float;
  4.   begin
  5.     if (x=0) then
  6.       begin
  7.         if y=0 then
  8.           arctan2:=0.0
  9.         else if y>0 then
  10.           arctan2:=pi/2
  11.         else if y<0 then
  12.           arctan2:=-pi/2;
  13.       end
  14.     else
  15.       ArcTan2:=ArcTan(y/x);
  16.     if x<0.0 then
  17.       ArcTan2:=ArcTan2+pi;
  18.     if ArcTan2>pi then
  19.       ArcTan2:=ArcTan2-2*pi;
  20.   end;
  21. {$endif FPC_MATH_HAS_ARCTAN2}
  22.  
Title: Re: Strange error with arccos()
Post by: BobDog on August 05, 2022, 02:44:46 pm
Purely for fun..
Code: Pascal  [Select][+][-]
  1. program compare;
  2.  uses
  3.  math;    // only for built in arctan2
  4.    
  5.  
  6.  Function MyAtan2(x:Double;y:Double):Double ;
  7.     Function MyATN(x:Double):Double;
  8.        var
  9.       c:longint=1;
  10.       counter:longint=0;
  11.       sign:longint;
  12.       xx,term,temp1,temp2,p,accum:double;
  13.  
  14.     procedure TaylorSeries(x:double);
  15.     begin
  16.     xx:=x*x;
  17.     p:=x;
  18.     temp1:=x;
  19.     term:=0;
  20.     repeat
  21.         c:=c+2;
  22.         temp2:=temp1;
  23.         counter:=counter+1;
  24.         p:=p*XX;
  25.         term:=p/c;
  26.         sign:=1;
  27.         If (counter And 1) <>0 Then sign:=-1;
  28.         accum:=temp1+sign*term;
  29.         temp1:=accum
  30.        Until (temp1=temp2)
  31.     end;
  32.  
  33.    var
  34.    y,temp:double;
  35.    limit:longint=16;
  36.    f:longint=65536;
  37.    n:longint;
  38.    begin
  39.      accum:=0.0;
  40.      temp:=x;
  41.      y:=0;
  42.     For n  :=1 To limit do
  43.     begin
  44.         y:=temp/(1+Sqrt(1+temp*temp));
  45.         temp:=y;
  46.     end;
  47.     TaylorSeries(y);
  48.     exit(f*accum);
  49.     end;
  50.  
  51.     function sgn(x:double):integer;
  52.     begin
  53.     if x>0 then exit(1);
  54.     if x<0 then exit(-1);
  55.     exit(0);
  56.     end;
  57.  
  58.     var
  59.     k:longint;
  60.     pi,temp:double;
  61.     begin
  62.     k:=0;
  63.     pi:= 4*myatn(1);
  64.     If (y=0) Then exit (Sgn(x)*pi/2);
  65.     k:=1;
  66.     If ((Sgn(y)=1) Or (Sgn(y)=0)) Then k:=0 ;
  67.     k:=k*Sgn(x) ;
  68.     If ((x=0) And (y=0)) Then exit(0.0);
  69.     If ((Sgn(x)=0) And (Sgn(y)=-1)) Then exit(Pi);
  70.     temp:=MyATN(x/y);
  71.     exit ((temp)+k*Pi);
  72. End;
  73.  
  74.  
  75.     var
  76.     x,y:double;
  77.     i:integer;
  78.  
  79.     begin
  80.     writeln('First principles','                                      ','Unit math');
  81.     for i:=1 to 30 do
  82.     begin
  83.     x:=random*100-random*100;
  84.     y:=random*100-random*100;
  85.     writeln(Myatan2(x,y),'               ',arctan2(x,y));
  86.     end;
  87.  
  88.     writeln('Press return to finish');
  89.     readln;
  90.     end.      
  91.      
  92.  
  93.  
Title: Re: Strange error with arccos()
Post by: Josh on August 05, 2022, 03:13:03 pm
@wp

I have slightly moded your code, so that the range is more within the min and max values for arctan function,
also stored the results in array for comparison of result.

Also i have limited ram so had to change how many N numbers

also moded my routine so the range check is same as the check used by the math routine.
moded paola routine to use pi/2 to elinate calculation differences.

prog give different readings...... and differences in calculated results.
Code: Pascal  [Select][+][-]
  1. program arctan2_speed;
  2.  
  3. uses
  4.   sysutils, math;
  5.  
  6. function PaoloArcTan2(Y, X : Float) : Float;
  7. begin
  8.   if (X = 0) then begin
  9.     if (Y = 0) then
  10.       Result:=0
  11.     else
  12.       if (Y > 0) then
  13.         Result:=pi/2
  14.       else
  15.         Result:=-pi/2
  16.   end
  17.   else begin
  18.     Result:=ArcTan(Y/X);
  19.     if (X < 0) then
  20.       if (Y < 0) then
  21.         Result:=Result-pi
  22.       else
  23.         Result:=Result+pi
  24.   end;
  25. end;
  26.  
  27. function JoshArcTan2 (y, x : Float) : Float;
  28. begin
  29.  result:=0.0; // set to 0.0 if it falls through with x and y =0.00 also stops hint result not initialized
  30.  if x>0.0 then result := arctan (y/x)                      // using /2 rather than *.5 as old school / was generally more accurate and for me more readable
  31.    else if x <0.0 then result := arctan (y/x) + pi
  32.      else if y<>0.0 then result := pi/2 * sign (y);   // x=0 so if y is also 0.0 result stays set at 0.00 else result gets set +-pi/2
  33.  
  34.  if result>pi then result:=result-2*pi;    // keep result in range of +- pi   could use abs sign here but then would be used every iteration of function
  35.  //  else if result<-pi then result:=result+2*pi;    // not sure if needed but safguarded like original code
  36.  // same check as math unit
  37. end;
  38.  
  39. const
  40.   MILLION = 1000 * 1000;
  41.   N = 5 * MILLION;
  42.  
  43. type
  44.   TArcTan2Func = function(y, x: Float): Float;
  45.  
  46.   TData = record
  47.     x, y: Float;
  48.   end;
  49.  
  50. var
  51.   Data: array[1..n] of TData;
  52.   j,m,p:Array[1..n] of Float;
  53.  
  54. procedure PrepareData;
  55. var
  56.   i: Integer;
  57. begin
  58.   RandSeed := 0;
  59.   for i := 1 to N do
  60.   begin
  61.     Data[i].X := (Random()-0.5) * high(dword);
  62.     Data[i].Y := (Random()-0.5) * high(dword);
  63.   end;
  64. end;
  65.  
  66. procedure Test(Descr: String; ArcTan2Func: TArcTan2Func; w:string);
  67. var
  68.   t: TDateTime;
  69.   i: Integer;
  70.   ans:float;
  71. begin
  72.   Write('Testing ' + Descr + '...');
  73.   RandSeed := 0;
  74.   t := Now;
  75.   for i := 1 to N do
  76.   begin
  77.     ans:=ArcTan2Func(Data[i].Y, Data[i].X);
  78.     case w of
  79.       'j':j[i]:=ans;
  80.       'm':m[i]:=ans;
  81.       'p':p[i]:=ans;
  82.     end;
  83.   end;
  84.   t := Now - t;
  85.   WriteLn(' done: ', FormatDateTime('s.zzz', t), ' sec');
  86. end;
  87. var l,c:integer;
  88. begin
  89.   PrepareData;
  90.   Test('Math.arctan2', @arctan2,'m');
  91.   Test('PaoloArctan2', @paoloarctan2,'p');
  92.   Test('JoshArctan2', @josharctan2,'j');
  93.  
  94.   WriteLn('Finished. Calculating');
  95.   WriteLn('Check difference between josh and math Press ENTER.');
  96.   readln;
  97.   c:=0;
  98.   for l:=1 to n do
  99.   begin
  100.     if j[l]<>m[l] then
  101.       begin
  102.         writeln(l,':',j[l]-m[l]);
  103.         inc(c);
  104.       end;
  105.   end;
  106.   WriteLn('Differences '+c.ToString);
  107.   WriteLn('Check difference between josh and paola Press ENTER.');
  108.   readln;
  109.   c:=0;
  110.   for l:=1 to n do
  111.   begin
  112.     if j[l]<>p[l] then
  113.       begin
  114.         writeln(l,':',j[l]-p[l]);
  115.         inc(c);
  116.       end;
  117.   end;
  118.   WriteLn('Differences '+c.ToString);
  119.   ReadLn;
  120. end.
  121.  
Title: Re: Strange error with arccos()
Post by: Paolo on August 05, 2022, 03:44:11 pm
@Josh, I am Paolo not Paola  :D

I'll check your code, thanks
Title: Re: Strange error with arccos()
Post by: Paolo on August 05, 2022, 04:19:59 pm
@Josh,

you are doing exactly the computation the code in math library did in the cases of (X<0) and (Y<0), you first add pi then subtract -2*pi, my code just do -pi.
Now I can argue that I do the minimal operation so the rounding should be better, but the difference is on the last digit, so who is right ?
Title: Re: Strange error with arccos()
Post by: Paolo on August 05, 2022, 04:28:55 pm
@wp I did not see any appreciable difference (See picture), at least on the run I saw.

I see that doing click in my code on math.Arctan2 I jump in the pascal code showed previously (no ASM).
Title: Re: Strange error with arccos()
Post by: linuxbob on August 05, 2022, 06:29:00 pm
Quote from: wp on August 05, 2022, 08:48:51 am
The best function to get the polar angle from cartesian coordinates (x, y) is the function arctan2(y, x) because it also finds the correct quadrant for the angle.

I'm trying to get from the angle and magnitude to the coordinates, not the other way round.
Title: Re: Strange error with arccos()
Post by: wp on August 05, 2022, 06:32:21 pm
Quote from: BobDog on August 05, 2022, 02:44:46 pm
Purely for fun..
Added this to my test application - see attachment. Of course, it is much slower than the others based on the built-in arctan, but it is an interesting alternative on systems where unit math is not applicable. I added one optimization: your code had calculated the value of pi with every call of the Arctan2 function, and this doubles the entire execution time. Now pi (called MyPI here) is calculated only once per program run. Very probably, further optimizations are possible.

These are my result for a 32-bit application on Win 11
Quote
Testing Math.arctan2... done: 0.818 sec
Testing PaoloArctan2... done: 1.057 sec
Testing JoshArctan2... done: 1.073 sec
Testing BobDogArctan2... done: 6.723 sec
Finished. Press ENTER to close.
and for a 64-bit application on Win 11:
Quote
Testing Math.arctan2... done: 0.446 sec
Testing PaoloArctan2... done: 0.471 sec
Testing JoshArctan2... done: 0.498 sec
Testing BobDogArctan2... done: 4.378 sec
Finished. Press ENTER to close.
Title: Re: Strange error with arccos()
Post by: linuxbob on August 05, 2022, 06:35:09 pm
Quote from: mtrsoft on August 05, 2022, 06:11:06 am
You are trying to take the square root of a negative number.

The algorithym that calculates arccos has a square root in it.

The "principal range" for arccos(y) is 0 <= y <= Pi and for arcsin(y) is -Pi/2 <= y <= Pi/2.

See: https://en.wikipedia.org/wiki/Inverse_trigonometric_functions

Regards,
John

John, you are correct mathematically. But, in the world of AC circuits we have both leading and lagging phase angles (from Pi/2 to -Pi/2) for the capacitive and inductive impedances respectively. That is the root of the calculations I need to do. Since there is a square root in the algorithm, I'll need to use a different approach.
Title: Re: Strange error with arccos()
Post by: wp on August 05, 2022, 06:36:29 pm
Quote from: linuxbob on August 05, 2022, 06:29:00 pm
Quote from: wp on August 05, 2022, 08:48:51 am
The best function to get the polar angle from cartesian coordinates (x, y) is the function arctan2(y, x) because it also finds the correct quadrant for the angle.

I'm trying to get from the angle and magnitude to the coordinates, not the other way round.
Then your initial post was very misleading because you don't need arccos and arcsin for this, just the ordinary cos and sin:
Code: Pascal  [Select][+][-]
  1. uses
  2.   Math;
  3. procedure PolarToCartesian(const r, phi: Double; var x, y: Double);
  4. var
  5.   sinphi, cosphi: Double;
  6. begin
  7.   SinCos(phi, sinphi, cosphi);
  8.   x := r * cosphi;
  9.   y := r * sinphi;
  10.   { faster than the straight-forward
  11.     x := r * cos(phi);
  12.     y := r * sin(phi); }
  13. end;
But probably I still don't understand the question correctly.
Title: Re: Strange error with arccos()
Post by: Paolo on August 05, 2022, 07:21:38 pm
Quote from: linuxbob on August 05, 2022, 06:35:09 pm
Quote from: mtrsoft on August 05, 2022, 06:11:06 am
You are trying to take the square root of a negative number.

The algorithym that calculates arccos has a square root in it.

The "principal range" for arccos(y) is 0 <= y <= Pi and for arcsin(y) is -Pi/2 <= y <= Pi/2.

See: https://en.wikipedia.org/wiki/Inverse_trigonometric_functions

Regards,
John

John, you are correct mathematically. But, in the world of AC circuits we have both leading and lagging phase angles (from Pi/2 to -Pi/2) for the capacitive and inductive impedances respectively. That is the root of the calculations I need to do. Since there is a square root in the algorithm, I'll need to use a different approach.

which algorithm are you using, show some code...
Title: Re: Strange error with arccos()
Post by: Paolo on August 05, 2022, 07:32:33 pm
@wp couriously me too, I have renamed Bobdog pi in Mypi (to differentiate wrt standard pi).

@BobDog,wp: Here attached the bobdog code revised, note that there are some semplification/corrections in the code (i am not went in deep in the code), one of them is the fact that local var counter is not initializated in the code (but strangely enough the complier does not report the usual warning).

@wp can you show what arctan2 code you see when jump to the routine, you said having win11 (64 bit ? I don't know) I am expecting the same behaviour for win10 64bit..
Title: Re: Strange error with arccos()
Post by: Paolo on August 05, 2022, 07:38:48 pm
@wp I see now you tests, sorry. Your test seems confirming that win32 is doing for ASM and win64 for "Pascal implementation", am I right ? in the case of 64bit all routines are close each other in speed (the mode of compliation may affect some performances (?)).

Anway in my view the pascal ARCTAN2 should be updated as suggested... (maybe it'll go faster !)

Title: Re: Strange error with arccos()
Post by: winni on August 05, 2022, 07:42:25 pm
Hi!

As mentioned above Linux64 is ASM.

Winni
Title: Re: Strange error with arccos()
Post by: Paolo on August 05, 2022, 07:51:12 pm
@winni, ok for linux ASM code is used.
Title: Re: Strange error with arccos()
Post by: wp on August 05, 2022, 07:59:18 pm
Quote from: Paolo on August 05, 2022, 07:38:48 pm
Anway in my view the pascal ARCTAN2 should be updated as suggested... (maybe it'll go faster !)
I don't see an essential speed gain, but the current code really looks "old-fashioned" Fortran-like. As I said: Write a bugreport, refer to this discussion and let the dev decide. The report should also contain a test program comparing the results of the methods (similar to what was shown above), also with corner cases x=y=0, x=0, y=0.
Title: Re: Strange error with arccos()
Post by: Josh on August 05, 2022, 08:09:07 pm
Difference in speed I think is too dependent on CPU / architecture /OS but results I get on my HP Laptop
Intel(R) Pentium(R) CPU 5405U @ 2.30GHz  win 11 64, 8GB RAM. Laz/FPC Trunk 64 bit.

Quote
Calculating and testing 20000000 samples
Testing Math.arctan2... done: 1.796 sec
Testing PaoloArctan2... done: 2.161 sec
Testing JoshArctan2... done: 1.690 sec
Finished. Calculating
Check difference between math and josh Press ENTER.

Differences 0
Total of Differences 0.0000000000000000

Check difference between math and Paolo Press ENTER.

Differences 1143105
Total of Differences 0.0000000003991967

Finished Press ENTER
Title: Re: Strange error with arccos()
Post by: Paolo on August 05, 2022, 08:19:08 pm
@Josh, In my testing I don't see too much difference (on both intel and AMD), but can you repeat the test with "constant" input X=Y=0 ?
Title: Re: Strange error with arccos()
Post by: Josh on August 05, 2022, 08:34:58 pm
Quote
Calculating and testing 20000000 random x & y samples
Testing Math.arctan2... done: 1.613 sec
Testing PaoloArctan2... done: 2.134 sec
Testing JoshArctan2... done: 1.463 sec
Finished. Calculating
Check difference between math and josh Press ENTER.

Differences 0
Total of Differences 0.0000000000000000

Check difference between math and Paolo Press ENTER.

Differences 1143105
Total of Differences 0.0000000003991967

Finished Press ENTER
Quote
Calculating and testing 20000000 samples with all x=0
Testing Math.arctan2... done: 1.036 sec
Testing PaoloArctan2... done: 1.215 sec
Testing JoshArctan2... done: 0.428 sec
Finished. Calculating
Check difference between math and josh Press ENTER.

Differences 0
Total of Differences 0.0000000000000000

Check difference between math and Paolo Press ENTER.

Differences 0
Total of Differences 0.0000000000000000

Finished Press ENTER
Quote
Calculating and testing 20000000 samples with all y=0
Testing Math.arctan2... done: 1.107 sec
Testing PaoloArctan2... done: 1.245 sec
Testing JoshArctan2... done: 0.709 sec
Finished. Calculating
Check difference between math and josh Press ENTER.

Differences 0
Total of Differences 0.0000000000000000

Check difference between math and Paolo Press ENTER.

Differences 0
Total of Differences 0.0000000000000000

Finished Press ENTER
Quote
Calculating and testing 20000000 samples with all y=0 and x=0
Testing Math.arctan2... done: 0.880 sec
Testing PaoloArctan2... done: 0.890 sec
Testing JoshArctan2... done: 0.466 sec
Finished. Calculating
Check difference between math and josh Press ENTER.

Differences 0
Total of Differences 0.0000000000000000

Check difference between math and Paolo Press ENTER.

Differences 0
Total of Differences 0.0000000000000000

Finished Press ENTER
Title: Re: Strange error with arccos()
Post by: Paolo on August 05, 2022, 08:59:20 pm
How do you explain the result for input x=y=0, my code just test x , y and exit, how do you explain your code is twice faster?
Title: Re: Strange error with arccos()
Post by: Josh on August 05, 2022, 09:44:40 pm
code and results below on my laptop
what are others result.
Code: Pascal  [Select][+][-]
  1. program arctan2_speed;
  2.  
  3. uses
  4.   sysutils, math;
  5.  
  6.  
  7.  
  8. function PaoloArcTan2(Y, X : Float) : Float;
  9. begin
  10.   if (X = 0) then begin
  11.     if (Y = 0) then
  12.       Result:=0
  13.     else
  14.       if (Y > 0) then
  15.         Result:=0.5*pi
  16.       else
  17.         Result:=-0.5*pi
  18.   end
  19.   else begin
  20.     Result:=ArcTan(Y/X);
  21.     if (X < 0) then
  22.       if (Y < 0) then
  23.         Result:=Result-pi
  24.       else
  25.         Result:=Result+pi
  26.   end;
  27. end;
  28.  
  29. function JoshArcTan2 (y, x : Float) : Float;
  30. begin
  31.  if x>0.0 then result := arctan (y/x)                      // using /2 rather than *.5 as old school / was generally more accurate and for me more readable
  32.    else if x <0.0 then result := arctan (y/x) + pi
  33.      else if y<>0.0 then result := pi/2 * sign (y)
  34.             else result:=0.0;// x=0 so if y is also 0.0 result stays set at 0.00 else result gets set +-pi/2
  35.  
  36.  if result>pi then result:=result-2*pi;    // keep result in range of +- pi   could use abs sign here but then would be used every iteration of function
  37. end;
  38.  
  39. const
  40.   MILLION = 1000 * 1000;
  41.   N = 20 * MILLION;
  42.  
  43. type
  44.   TArcTan2Func = function(y, x: Float): Float;
  45.  
  46.   TData = record
  47.     x, y: Float;
  48.   end;
  49.  
  50. var
  51.   Data: array[1..n] of TData;
  52.   j,m,p:Array[1..n] of Float;
  53.  
  54. procedure PrepareData(fixx,fixy:boolean;v:float);
  55. var
  56.   i: Integer;
  57. begin
  58.  // randomize;
  59.   RandSeed := 0;
  60.   for i := 1 to N do
  61.   begin
  62.     if fixx then Data[i].X:=v else Data[i].X := (Random()-0.5) * high(dword);
  63.     if fixy then Data[i].Y:=v else Data[i].Y := (Random()-0.5) * high(dword);
  64.   end;
  65. end;
  66.  
  67. procedure Test(Descr: String; ArcTan2Func: TArcTan2Func; w:string);
  68. var
  69.   t: TDateTime;
  70.   i: Integer;
  71.   ans:float;
  72. begin
  73.   Write('Testing ' + Descr + '...');
  74.   RandSeed := 0;
  75.   t := Now;
  76.   for i := 1 to N do
  77.   begin
  78.     ans:=ArcTan2Func(Data[i].Y, Data[i].X);
  79.     case w of
  80.       'j':j[i]:=ans;
  81.       'm':m[i]:=ans;
  82.       'p':p[i]:=ans;
  83.     end;
  84.   end;
  85.   t := Now - t;
  86.   WriteLn(' done: ', FormatDateTime('s.zzz', t), ' sec');
  87. end;
  88. var l,c:integer;
  89.   td:extended;
  90.  
  91. procedure reportresult(heading:string);
  92. begin
  93.   WriteLn(heading);
  94.   Test('Math.arctan2', @arctan2,'m');
  95.   Test('PaoloArctan2', @paoloarctan2,'p');
  96.   Test('JoshArctan2', @josharctan2,'j');
  97.   WriteLn('Finished. Calculating');
  98.   Write('Check difference between math and Josh');
  99.   c:=0;
  100.   td:=0;
  101.   for l:=1 to n do
  102.   begin
  103.     if j[l]<>m[l] then
  104.       begin
  105.         inc(c);
  106.         td:=td+abs(j[l]-m[l]);
  107.       end;
  108.   end;
  109.   if c=0 then WriteLn(' No Differences')
  110.   else
  111.   begin
  112.     WriteLn;
  113.     WriteLn('     Differences '+c.ToString+'Total of Differences ',td:4:16);
  114.   end;
  115.   Write('Check difference between math and Paolo');
  116.   c:=0;
  117.   td:=0;
  118.   for l:=1 to n do
  119.   begin
  120.     if m[l]<>p[l] then
  121.       begin
  122.         inc(c);
  123.         td:=td+abs(p[l]-m[l]);
  124.       end;
  125.   end;
  126.   if c=0 then WriteLn(' No Differences')
  127.   else
  128.   begin
  129.     WriteLn;
  130.     WriteLn('     Differences '+c.ToString+'Total of Differences ',td:4:16);
  131.   end;
  132.   Writeln('** End of Test ** '+heading);
  133.   Writeln;Writeln;
  134. end;
  135.  
  136. begin
  137.   randomize;
  138.   PrepareData(false,false,0);
  139.   reportresult('Calculating and testing '+n.ToString+' random x & y samples');
  140.   PrepareData(true,false,0);
  141.   reportresult('Calculating and testing '+n.ToString+' samples with all x=0');
  142.   PrepareData(false,true,0);
  143.   reportresult('Calculating and testing '+n.ToString+' samples with all y=0');
  144.   PrepareData(true,true,0);
  145.   reportresult('Calculating and testing '+n.ToString+' samples with all y=0 and x=0');
  146.   ReadLn;
  147. end.
  148.  

result from above code
Quote
Calculating and testing 20000000 random x & y samples
Testing Math.arctan2... done: 1.527 sec
Testing PaoloArctan2... done: 1.858 sec
Testing JoshArctan2... done: 1.165 sec
Finished. Calculating
Check difference between math and Josh No Differences
Check difference between math and Paolo
     Differences 1143105Total of Differences 0.0000000003991967
** End of Test ** Calculating and testing 20000000 random x & y samples


Calculating and testing 20000000 samples with all x=0
Testing Math.arctan2... done: 0.694 sec
Testing PaoloArctan2... done: 0.880 sec
Testing JoshArctan2... done: 0.275 sec
Finished. Calculating
Check difference between math and Josh No Differences
Check difference between math and Paolo No Differences
** End of Test ** Calculating and testing 20000000 samples with all x=0


Calculating and testing 20000000 samples with all y=0
Testing Math.arctan2... done: 0.764 sec
Testing PaoloArctan2... done: 0.976 sec
Testing JoshArctan2... done: 0.465 sec
Finished. Calculating
Check difference between math and Josh No Differences
Check difference between math and Paolo No Differences
** End of Test ** Calculating and testing 20000000 samples with all y=0


Calculating and testing 20000000 samples with all y=0 and x=0
Testing Math.arctan2... done: 0.507 sec
Testing PaoloArctan2... done: 0.727 sec
Testing JoshArctan2... done: 0.191 sec
Finished. Calculating
Check difference between math and Josh No Differences
Check difference between math and Paolo No Differences
** End of Test ** Calculating and testing 20000000 samples with all y=0 and x=0

I assume to check for >0.00 the compiler would just check that the MSB is set to on; to evaluate the expression, to compare with a specific value it would need to check the whole value for comparison 8 bytes.
I always where possible use this technique <0 or >0 as it is generally much faster,even with integers.
Title: Re: Strange error with arccos()
Post by: d4eva on August 05, 2022, 10:18:41 pm
Results on Mac M1 mini compiled as native aarch64

Calculating and testing 20000000 random x & y samples
Testing Math.arctan2... done: 0.761 sec
Testing PaoloArctan2... done: 0.924 sec
Testing JoshArctan2... done: 0.596 sec
Finished. Calculating
Check difference between math and Josh No Differences
Check difference between math and Paolo
     Differences 1141138Total of Differences 0.0000000003992473
** End of Test ** Calculating and testing 20000000 random x & y samples


Calculating and testing 20000000 samples with all x=0
Testing Math.arctan2... done: 0.339 sec
Testing PaoloArctan2... done: 0.493 sec
Testing JoshArctan2... done: 0.098 sec
Finished. Calculating
Check difference between math and Josh No Differences
Check difference between math and Paolo No Differences
** End of Test ** Calculating and testing 20000000 samples with all x=0


Calculating and testing 20000000 samples with all y=0
Testing Math.arctan2... done: 0.376 sec
Testing PaoloArctan2... done: 0.523 sec
Testing JoshArctan2... done: 0.239 sec
Finished. Calculating
Check difference between math and Josh No Differences
Check difference between math and Paolo No Differences
** End of Test ** Calculating and testing 20000000 samples with all y=0


Calculating and testing 20000000 samples with all y=0 and x=0
Testing Math.arctan2... done: 0.238 sec
Testing PaoloArctan2... done: 0.375 sec
Testing JoshArctan2... done: 0.095 sec
Finished. Calculating
Check difference between math and Josh No Differences
Check difference between math and Paolo No Differences
** End of Test ** Calculating and testing 20000000 samples with all y=0 and x=0
Title: Re: Strange error with arccos()
Post by: BobDog on August 05, 2022, 10:41:36 pm
Quote from: wp on August 05, 2022, 06:32:21 pm
Quote from: BobDog on August 05, 2022, 02:44:46 pm
Purely for fun..
Added this to my test application - see attachment. Of course, it is much slower than the others based on the built-in arctan, but it is an interesting alternative on systems where unit math is not applicable. I added one optimization: your code had calculated the value of pi with every call of the Arctan2 function, and this doubles the entire execution time. Now pi (called MyPI here) is calculated only once per program run. Very probably, further optimizations are possible.

These are my result for a 32-bit application on Win 11
Quote
Testing Math.arctan2... done: 0.818 sec
Testing PaoloArctan2... done: 1.057 sec
Testing JoshArctan2... done: 1.073 sec
Testing BobDogArctan2... done: 6.723 sec
Finished. Press ENTER to close.
and for a 64-bit application on Win 11:
Quote
Testing Math.arctan2... done: 0.446 sec
Testing PaoloArctan2... done: 0.471 sec
Testing JoshArctan2... done: 0.498 sec
Testing BobDogArctan2... done: 4.378 sec
Finished. Press ENTER to close.
Hi wp
I have set the counter at zero in my code.
Thanks for noticing that.
Couldn't you use (for speed)
 const MyPi=3.141592653589793;
just above the var near the function end.

Title: Re: Strange error with arccos()
Post by: Paolo on August 05, 2022, 10:45:51 pm
@Josh and 4deva,

try to modify the case test sequence

from
Code: Pascal  [Select][+][-]
  1.     case w of
  2.       'j':j[i]:=ans;
  3.       'm':m[i]:=ans;
  4.       'p':p[i]:=ans;
  5.     end;
  6.  

to
Code: Pascal  [Select][+][-]
  1.     case w of
  2.       'p':p[i]:=ans;
  3.       'j':j[i]:=ans;
  4.       'm':m[i]:=ans;
  5.     end;
  6.  

and tell me the results :)
Title: Re: Strange error with arccos()
Post by: Paolo on August 05, 2022, 11:06:54 pm
@Josh, trying to compile in "debug" mode I see error in the picture ("default" and "Release" are fine)

I did a working copy...
Title: Re: Strange error with arccos()
Post by: jamie on August 06, 2022, 03:42:00 am
For a good time in bug hunting..

 ArcCos(YourValue Mod 1.0);

 :o
Title: Re: Strange error with arccos()
Post by: d4eva on August 06, 2022, 08:33:01 am
Results after changing the switch statement:

Calculating and testing 20000000 random x & y samples
Testing Math.arctan2... done: 0.937 sec
Testing PaoloArctan2... done: 0.627 sec
Testing JoshArctan2... done: 0.731 sec
Finished. Calculating
Check difference between math and Josh No Differences
Check difference between math and Paolo
     Differences 1141138Total of Differences 0.0000000003992473
** End of Test ** Calculating and testing 20000000 random x & y samples


Calculating and testing 20000000 samples with all x=0
Testing Math.arctan2... done: 0.499 sec
Testing PaoloArctan2... done: 0.188 sec
Testing JoshArctan2... done: 0.238 sec
Finished. Calculating
Check difference between math and Josh No Differences
Check difference between math and Paolo No Differences
** End of Test ** Calculating and testing 20000000 samples with all x=0


Calculating and testing 20000000 samples with all y=0
Testing Math.arctan2... done: 0.532 sec
Testing PaoloArctan2... done: 0.209 sec
Testing JoshArctan2... done: 0.393 sec
Finished. Calculating
Check difference between math and Josh No Differences
Check difference between math and Paolo No Differences
** End of Test ** Calculating and testing 20000000 samples with all y=0


Calculating and testing 20000000 samples with all y=0 and x=0
Testing Math.arctan2... done: 0.395 sec
Testing PaoloArctan2... done: 0.095 sec
Testing JoshArctan2... done: 0.232 sec
Finished. Calculating
Check difference between math and Josh No Differences
Check difference between math and Paolo No Differences
** End of Test ** Calculating and testing 20000000 samples with all y=0 and x=0


Title: Re: Strange error with arccos()
Post by: Paolo on August 06, 2022, 09:02:30 am
@d4eva, thanks to check.
This explain the huge difference in the results.
Any speed test should avoid things that interfere with the measurement. I have also the suspect on the array dimension playing some role.
Title: Re: Strange error with arccos()
Post by: Josh on August 06, 2022, 07:29:23 pm
hi

dont normally use internet on weekends..  ;)

yep case is expensive (more than i thought), so replaced with if then statements, so everytime the same comparisons are done, no favour to routines.
.

Also moded test app, so that it repeats the samples unitil stops, but keeps a running total and average time.
The app also nows has a copy of the function from the math unit, so that any optimizations are kept the same.
also moded the matharctan2 function so that it does not -2pi the result; this has the effect of correcting the rounding error and all function generate equal data.

on line 26 if you comment out the define it will compile as factory. ie add // before the define

on my laptop the longer it runs the more they even out...

images are runs of 200+ so each routine visited over 4 Billion Times.



Title: Re: Strange error with arccos()
Post by: Paolo on August 06, 2022, 08:16:14 pm
Quote from: Josh on August 06, 2022, 07:29:23 pm
dont normally use internet on weekends..  ;)
..and it is summer time..  :D, .. and arctan can wait  :D

Can you post your latest version of code ? From my testing I see something still not clear.

Your conclusion on the basis of such results ?
Title: Re: Strange error with arccos()
Post by: Josh on August 06, 2022, 08:29:02 pm
test code is added in previous post
Title: Re: Strange error with arccos()
Post by: Paolo on August 09, 2022, 02:47:38 pm
Done some testing, some improvements are observed even if very little.
Title: Re: Strange error with arccos()
Post by: Josh 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.
Title: Re: Strange error with arccos()
Post by: Paolo on August 09, 2022, 03:43:18 pm
I'll submit the second version, for the same reasons you exposed. I'll check tanh !
Title: Re: Strange error with arccos()
Post by: Paolo 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 ?
Title: Re: Strange error with arccos()
Post by: Paolo 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.
Title: Re: Strange error with arccos()
Post by: Josh 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?
Title: Re: Strange error with arccos()
Post by: wp on August 10, 2022, 03:49:08 pm
Quote from: Josh on August 10, 2022, 02:18:58 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.  
Title: Re: Strange error with arccos()
Post by: Paolo 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
Title: Re: Strange error with arccos()
Post by: Josh 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.

Title: Re: Strange error with arccos()
Post by: Paolo 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
Title: Re: Strange error with arccos()
Post by: Paolo on August 10, 2022, 06:40:02 pm
for the arctan2

https://gitlab.com/freepascal.org/fpc/source/-/issues/39861
Title: Re: Strange error with arccos()
Post by: wp on August 10, 2022, 06:47:27 pm
Quote from: Paolo on August 10, 2022, 06:40:02 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.
Title: Re: Strange error with arccos()
Post by: Paolo 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
Title: Re: Strange error with arccos()
Post by: Paolo 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
 (https://gitlab.com/freepascal.org/fpc/source/-/issues/39867)
TinyPortal © 2005-2018