Recent

Author Topic: Arctan2  (Read 12393 times)

Paolo

  • Sr. Member
  • ****
  • Posts: 499
Arctan2
« on: August 02, 2021, 06:49:09 pm »
Hello,

following the discussion in https://forum.lazarus.freepascal.org/index.php/topic,54902.msg408121.html#msg408121, I understood that in win-64 arctan2 is computed as follow.

Code: Pascal  [Select][+][-]
  1. {$ifndef FPC_MATH_HAS_ARCTAN2}
  2. function arctan2(y,x : float) : float;
  3.   begin
  4.     if (x=0) then
  5.       begin
  6.         if y=0 then
  7.           arctan2:=0.0
  8.         else if y>0 then
  9.           arctan2:=pi/2
  10.         else if y<0 then
  11.           arctan2:=-pi/2;
  12.       end
  13.     else
  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. {$endif FPC_MATH_HAS_ARCTAN2}
  21.  

However looking on internet I found this

Code: Pascal  [Select][+][-]
  1. function PartArcTan(X: double): double;
  2. begin
  3. asm
  4. {$IFDEF CPUX64}
  5.         movq [rsp-8], xmm0
  6.         fld    qword ptr [rsp-8]
  7. {$ELSE}
  8.         fld    qword ptr X
  9. {$ENDIF}
  10.         fld1
  11.         fpatan
  12.         fwait
  13. {$IFDEF CPUX64}
  14.         fstp   qword ptr [rsp-8]
  15.         movq   xmm0, [rsp-8]
  16. {$ENDIF}
  17. end;
  18. end;
  19.  

And I adapt it as follow to work on all four quadrants

Code: Pascal  [Select][+][-]
  1. function PartArcTan2(Y, X: double): double;
  2. begin
  3. {$ASMMODE intel}
  4. asm
  5.         movq [rsp-8], xmm0
  6.         movq [rsp-16], xmm1
  7.         fld    qword ptr [rsp-8]
  8.         fld    qword ptr [rsp-16]
  9.         fpatan
  10.         fwait
  11.         fstp   qword ptr [rsp-8]
  12.         movq   xmm0, [rsp-8]
  13. end;
  14. end;  
  15.  

Being not so skilled in ASM I would like to ask (I am interested on Win64 solution):
- My code is correct ?
- Is it worth to do that ? I mean for example in term of speed or it does not make sense to do that ?
- I have noted that for X=Y=0 arctan=0 and if X=0 the arctan is +-90 according to the sign of Y (so X=0 is not raising any exception, in some way this make sense considering divion by zero at the limit +-inf that happens for +-90)

Thank you in advance.

win64, laz 2.0.12, fpc 3.2.0

winni

  • Hero Member
  • *****
  • Posts: 3197
Re: Arctan2
« Reply #1 on: August 02, 2021, 07:39:05 pm »
Hi!

The unit math for 64 bit (win/lin) shows just this

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

Kays

  • Hero Member
  • *****
  • Posts: 569
  • Whasup!?
    • KaiBurghardt.de
Re: Arctan2
« Reply #2 on: August 02, 2021, 08:35:44 pm »
Being not so skilled in ASM I would like to ask (I am interested on Win64 solution):
- My code is correct ?
Well, does it produce the right result? (Rhetorical question.)Considering all of this, and I see winni just posted that, you arrive at
Code: Pascal  [Select][+][-]
  1. {$asmMode Intel}
  2. function arcTangent(const x, y: valReal): valReal; assembler;
  3. asm
  4.         fld x  // x
  5.         fld y  // y x
  6.         fpatan // arctan(x) / y
  7. end;

- Is it worth to do that ? I mean for example in term of speed or it does not make sense to do that ?
Frankly? No. It’s too trivial. If you have some complicated calculations requiring like all 8 FPU stack registers, then maaaaayyyybe.

- I have noted that for X=Y=0 arctan=0 and if X=0 the arctan is +-90 according to the sign of Y (so X=0 is not raising any exception, in some way this make sense considering divion by zero at the limit +-inf that happens for +-90)
It’s all documented. Take a look at the instruction reference (the comment below the table).
Yours Sincerely
Kai Burghardt

Paolo

  • Sr. Member
  • ****
  • Posts: 499
Re: Arctan2
« Reply #3 on: August 02, 2021, 08:46:26 pm »
Thanks to all,

here my notes:

- the routine shown by Winni is grayed and I assume that the code executed is the first one shown by me being not gayed
- I copy and paste the code as winni/Kays shown in the main code: here the worning: "unrecognized opcode FLDT (or FLD)"

so, if I missed something, how can I execute the code shown by winni/Kays ?



« Last Edit: August 02, 2021, 08:50:25 pm by Paolo »

winni

  • Hero Member
  • *****
  • Posts: 3197
Re: Arctan2
« Reply #4 on: August 02, 2021, 09:04:51 pm »
Thanks to all,

here my notes:

- the routine shown by Winni is grayed and I assume that the code executed is the first one shown by me being not gayed

If your code is grayed you are on a 32 bit installation



Paolo

  • Sr. Member
  • ****
  • Posts: 499
Re: Arctan2
« Reply #5 on: August 02, 2021, 09:10:33 pm »
I am on win64, and installed lazarus for 64bit code.

PascalDragon

  • Hero Member
  • *****
  • Posts: 5444
  • Compiler Developer
Re: Arctan2
« Reply #6 on: August 03, 2021, 09:00:14 am »
Thanks to all,

here my notes:

- the routine shown by Winni is grayed and I assume that the code executed is the first one shown by me being not gayed

If your code is grayed you are on a 32 bit installation

Or on Win64, cause Win64 does not use the FPU.


It is needed if one wants a potential FPU exception raised at the correct location. Otherwise it will be raised at the next FPU instruction which might be far away.

  • If you haven’t altered the calling convention, the default ABI on Linux (and maybe on Windoze, I don’t know) allows you to pass the result of an real function on top of the FPU stack, so you don’t need to move anything anywhere. Just place the result in st(0).

On Win64 the FPU stack is not used, thus no non-assembly code calling that function will check there.

winni

  • Hero Member
  • *****
  • Posts: 3197
Re: Arctan2
« Reply #7 on: August 03, 2021, 09:25:05 am »

Or on Win64, cause Win64 does not use the FPU.


On Windows 7 / 64 the code is NOT grayed.

Winni

Paolo

  • Sr. Member
  • ****
  • Posts: 499
Re: Arctan2
« Reply #8 on: August 03, 2021, 10:24:12 am »
I am on windows 10 64bit, the code is grayed.

BobDog

  • Sr. Member
  • ****
  • Posts: 394
Re: Arctan2
« Reply #9 on: August 03, 2021, 12:57:30 pm »

Fear not people.
I can test arctan2 from the math unit with the real deal.
Code: Pascal  [Select][+][-]
  1. program compare;
  2.  uses
  3.  math;    // only for built in arctan2
  4.    
  5. Function MyATN(x:Double):Double;
  6.        var
  7.       c:longint=1;
  8.       counter,sign:longint;
  9.       xx,term,temp1,temp2,p,accum:double;
  10.  
  11. procedure TaylorSeries(x:double);
  12.     begin
  13.     xx:=x*x;
  14.     p:=x;
  15.     temp1:=x;
  16.     term:=0;
  17.     repeat
  18.         c:=c+2;
  19.         temp2:=temp1;
  20.         counter:=counter+1;
  21.         p:=p*XX;
  22.         term:=p/c;
  23.         sign:=1;
  24.         If (counter And 1) <>0 Then sign:=-1;
  25.         accum:=temp1+sign*term;
  26.         temp1:=accum
  27.        Until (temp1=temp2)
  28.     end;
  29.  
  30.    var
  31.    y,temp:double;
  32.    limit:longint=16;
  33.    f:longint=65536;
  34.    n:longint;
  35.    begin
  36.      accum:=0.0;
  37.      temp:=x;
  38.      y:=0;
  39.     For n  :=1 To limit do
  40.     begin
  41.         y:=temp/(1+Sqrt(1+temp*temp));
  42.         temp:=y;
  43.     end;
  44.     TaylorSeries(y);
  45.     exit(f*accum);
  46.     end;
  47.  
  48.     Function MyAtan2(x:Double;y:Double):Double ;
  49.  
  50.     function sgn(x:double):integer;
  51.     begin
  52.     if x>0 then exit(1);
  53.     if x<0 then exit(-1);
  54.     exit(0);
  55.     end;
  56.  
  57.     var
  58.     k:longint;
  59.     pi,temp:double;
  60.     begin
  61.     k:=0;
  62.     pi:= 4*myatn(1);
  63.     If (y=0) Then exit (Sgn(x)*pi/2);
  64.     k:=1;
  65.     If ((Sgn(y)=1) Or (Sgn(y)=0)) Then k:=0 ;
  66.     k:=k*Sgn(x) ;
  67.     If ((x=0) And (y=0)) Then exit(0.0);
  68.     If ((Sgn(x)=0) And (Sgn(y)=-1)) Then exit(Pi);
  69.     temp:=MyATN(x/y);
  70.     exit ((temp)+k*Pi);
  71. End;
  72.  
  73.     var
  74.     x,y:double;
  75.     i:integer;
  76.  
  77.     begin
  78.     writeln('Scratch','                                      ','Unit math');
  79.     for i:=1 to 20 do
  80.     begin
  81.     x:=random*100-random*100;
  82.     y:=random*100-random*100;
  83.     writeln(Myatan2(x,y),'               ',arctan2(x,y));
  84.     end;
  85.  
  86.     writeln('Press return to finish');
  87.     readln;
  88.     end.      
  89.      
  90.  

PascalDragon

  • Hero Member
  • *****
  • Posts: 5444
  • Compiler Developer
Re: Arctan2
« Reply #10 on: August 03, 2021, 01:35:57 pm »

Or on Win64, cause Win64 does not use the FPU.


On Windows 7 / 64 the code is NOT grayed.

Grayed or not: the code is not used on Win64, cause it's protected by FPC_HAS_TYPE_EXTENDED which is not set on Win64, thus the generic Pascal implementation will be used there.

Paolo

  • Sr. Member
  • ****
  • Posts: 499
Re: Arctan2
« Reply #11 on: August 04, 2021, 03:35:16 pm »
besides the ASM code, in the pascal version of arctan2 (here shown) I think is not well written, but I can be wrong  O:-)

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    <- this "if" it is uncessary, already checked Y=0 and Y>0
  10.           arctan2:=-pi/2;
  11.       end
  12.     else
  13.       ArcTan2:=ArcTan(y/x);
  14.     if x<0.0 then      <- this part of code is excuted also on the result of "if X=0 then".. that is not necessary, it should be only in the else part within a "Begin..End"
  15.       ArcTan2:=ArcTan2+pi;
  16.     if ArcTan2>pi then
  17.       ArcTan2:=ArcTan2-2*pi;
  18.   end;
  19.  
  20.  

here my version

Code: Pascal  [Select][+][-]
  1. function MyArcTan2(Y, X : double) : double;
  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.  

@Kays: yes the ASM code I showed works ! useful or not I have still  to check...

@BobDog : not tested... but at least in the fragment of your code below, I can say that you never  reach the row n.5 in case y=0 being filtered by row n.1 (PS: you make me crazy having swapped X<->Y as usually intended...)

Code: Pascal  [Select][+][-]
  1.     If (y=0) Then exit (Sgn(x)*pi/2);
  2.     k:=1;
  3.     If ((Sgn(y)=1) Or (Sgn(y)=0)) Then k:=0 ;
  4.     k:=k*Sgn(x) ;
  5.     If ((x=0) And (y=0)) Then exit(0.0);
  6.  

@PascalDragon
sorry for my ignorance by why the code below compile
Code: Pascal  [Select][+][-]
  1. function PartArcTan2(Y, X: double): double;
  2. begin
  3. {$ASMMODE intel}
  4. asm
  5.         movq [rsp-8], xmm0
  6.         movq [rsp-16], xmm1
  7.         fld    qword ptr [rsp-8]
  8.         fld    qword ptr [rsp-16]
  9.         fpatan
  10.         fwait
  11.         fstp   qword ptr [rsp-8]
  12.         movq   xmm0, [rsp-8]
  13. end;
  14. end;
  15.  
..and the code from Kays and Winni does not compile (compiler message "un reconignized op code FLDT or FLD", in my ASM version I have FLD op code !)

thanks for the help to everybody

BobDog

  • Sr. Member
  • ****
  • Posts: 394
Re: Arctan2
« Reply #12 on: August 04, 2021, 05:30:59 pm »
Hi Paolo
Yes you are right, the case 0,0 is handled before that line
Code: Pascal  [Select][+][-]
  1.  
  2.     uses
  3.     math;
  4.  
  5.     Function MyATN(x:Double):Double;
  6.        var
  7.       c:longint=1;
  8.       counter,sign:longint;
  9.       xx,term,temp1,temp2,p,accum:double;
  10.  
  11. procedure TaylorSeries(x:double);
  12.     begin
  13.     xx:=x*x;
  14.     p:=x;
  15.     temp1:=x;
  16.     term:=0;
  17.     repeat
  18.         c:=c+2;
  19.         temp2:=temp1;
  20.         counter:=counter+1;
  21.         p:=p*XX;
  22.         term:=p/c;
  23.         sign:=1;
  24.         If (counter And 1) <>0 Then sign:=-1;
  25.         accum:=temp1+sign*term;
  26.         temp1:=accum
  27.        Until (temp1=temp2)
  28.     end;
  29.  
  30.    var
  31.    y,temp:double;
  32.    limit:longint=16;
  33.    f:longint=65536;
  34.    n:longint;
  35.    begin
  36.      accum:=0.0;
  37.      temp:=x;
  38.      y:=0;
  39.     For n  :=1 To limit do
  40.     begin
  41.         y:=temp/(1+Sqrt(1+temp*temp));
  42.         temp:=y;
  43.     end;
  44.     TaylorSeries(y);
  45.     exit(f*accum);
  46.     end;
  47.  
  48. Function MyAtan2(x:Double;y:Double):Double ;
  49.     function sgn(x:double):integer;
  50.     begin
  51.     if x>0 then exit(1);
  52.     if x<0 then exit(-1);
  53.     exit(0);
  54.     end;
  55.     var
  56.     k:longint=1;
  57.     pi,temp:double;
  58.     begin
  59.     pi:= 4*myatn(1);
  60.     If (y=0) Then exit (Sgn(x)*pi/2);
  61.     If ((Sgn(y)=1) Or (Sgn(y)=0)) Then k:=0 ;
  62.     k:=k*Sgn(x) ;
  63.     If ((Sgn(x)=0) And (Sgn(y)=-1)) Then  exit(pi);
  64.     temp:=MyATN(x/y);
  65.     exit ((temp)+k*Pi);
  66. End;
  67.  
  68.  
  69. function range(mymin:longint;mymax:longint):longint;
  70.          begin
  71.              range:=trunc(int((Random*(Mymax-mymin+1)))+MyMin);
  72.          end;
  73.  
  74.  var
  75.     x,y:longint;
  76.     i:integer;
  77.  
  78.     begin
  79.     writeln('Scratch','                                      ','Unit math');
  80.     for i:=1 to 20 do
  81.     begin
  82.     x:=range(-2,2);
  83.     y:=range(-2,2);
  84.     writeln(Myatan2(x,y),'       ',arctan2(x,y),'  = arctangent2(',x,',',y,')');
  85.     end;
  86.     writeln('Press return to finish');
  87.     readln;
  88.     end.  
  89.      
« Last Edit: August 04, 2021, 06:25:43 pm by BobDog »

avk

  • Hero Member
  • *****
  • Posts: 752
Re: Arctan2
« Reply #13 on: August 05, 2021, 03:14:03 pm »
  ...
...sorry for my ignorance by why the code below compile
...
..and the code from Kays and Winni does not compile (compiler message "un reconignized op code FLDT or FLD", in my ASM version I have FLD op code !)
...

It seems that FPC doesn't really bother with assembly routines and just substitutes the corresponding registers instead of the encountered parameters.
Something like this:
Code: Pascal  [Select][+][-]
  1. {$asmmode intel}
  2. function ArcTan2(Y, X : Double) : Double; assembler;
  3. asm
  4.   fld x  // --> fld xmm1
  5.   fld y  // --> fld xmm0
  6.   fpatan
  7. end;  
  8.  
And since there is no FLD variant with an XMM register, the compiler cheerfully reports an error.

The whole question is whether it is a bug or a feature.
« Last Edit: August 05, 2021, 03:15:46 pm by avk »

PascalDragon

  • Hero Member
  • *****
  • Posts: 5444
  • Compiler Developer
Re: Arctan2
« Reply #14 on: August 06, 2021, 03:23:46 pm »
The whole question is whether it is a bug or a feature.

This is by design. When writing assembly routines you must know how parameters are passed, cause it might make a difference for instructions. The compiler will simply substitute the location of the parameter.

 

TinyPortal © 2005-2018