Recent

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

linuxbob

  • New member
  • *
  • Posts: 8
Strange error with arccos()
« 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?

mtrsoft

  • New Member
  • *
  • Posts: 38
Re: Strange error with arccos()
« Reply #1 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

Paolo

  • Sr. Member
  • ****
  • Posts: 333
Re: Strange error with arccos()
« Reply #2 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.
« Last Edit: August 05, 2022, 08:45:58 am by Paolo »

wp

  • Hero Member
  • *****
  • Posts: 10297
Re: Strange error with arccos()
« Reply #3 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.

Paolo

  • Sr. Member
  • ****
  • Posts: 333
Re: Strange error with arccos()
« Reply #4 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.

wp

  • Hero Member
  • *****
  • Posts: 10297
Re: Strange error with arccos()
« Reply #5 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;

Paolo

  • Sr. Member
  • ****
  • Posts: 333
Re: Strange error with arccos()
« Reply #6 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.  

Paolo

  • Sr. Member
  • ****
  • Posts: 333
Re: Strange error with arccos()
« Reply #7 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 ?

Paolo

  • Sr. Member
  • ****
  • Posts: 333
Re: Strange error with arccos()
« Reply #8 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.  

 :)

wp

  • Hero Member
  • *****
  • Posts: 10297
Re: Strange error with arccos()
« Reply #9 on: August 05, 2022, 10:30:24 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.

Josh

  • Hero Member
  • *****
  • Posts: 1024
Re: Strange error with arccos()
« Reply #10 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
« Last Edit: August 05, 2022, 12:19:49 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: 333
Re: Strange error with arccos()
« Reply #11 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 !

wp

  • Hero Member
  • *****
  • Posts: 10297
Re: Strange error with arccos()
« Reply #12 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.

winni

  • Hero Member
  • *****
  • Posts: 3165
Re: Strange error with arccos()
« Reply #13 on: August 05, 2022, 01:58:30 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

Paolo

  • Sr. Member
  • ****
  • Posts: 333
Re: Strange error with arccos()
« Reply #14 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.  

 

TinyPortal © 2005-2018