Strange error with arccos()

[BobDog]

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.  

[Josh]

@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.  

[Paolo]

@Josh, I am Paolo not Paola 

I'll check your code, thanks

[Paolo]

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

[Paolo]

@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).

[linuxbob]

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.

[wp]

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

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:

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.

[linuxbob]

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.

[wp]

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.

[Paolo]

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...

[Paolo]

@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..

[Paolo]

@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 !)

[winni]

Hi!

As mentioned above Linux64 is ASM.

Winni

[Paolo]

@winni, ok for linux ASM code is used.

[wp]

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.