Strange error with arccos()

[linuxbob]

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?

[mtrsoft]

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

[Paolo]

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.

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

[Paolo]

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.

[wp]

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;

[Paolo]

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

[Paolo]

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 ?

[Paolo]

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

 

[wp]

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.

[Josh]

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

[Paolo]

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

[wp]

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.

Testing Math.arctan2... done: 0.816 sec
Testing PaoloArctan2... done: 1.062 sec
Testing JoshArctan2... done: 1.079 sec
Finished. Press ENTER to close.

[winni]

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

[Paolo]

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