### Bookstore

 Computer Math and Games in Pascal (preview) Lazarus Handbook

### Author Topic: Arctan2  (Read 10769 times)

#### Paolo

• Full Member
• Posts: 222
##### 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.

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)

win64, laz 2.0.12, fpc 3.2.0

#### winni

• Hero Member
• Posts: 2667
##### 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

• Sr. Member
• Posts: 357
• Whasup!?
##### 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

• Full Member
• Posts: 222
##### 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: 2667
##### 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

• Full Member
• Posts: 222
##### 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: 3335
• 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: 2667
##### 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

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

#### BobDog

• Full Member
• Posts: 132
##### 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.
88.     end.
89.
90.

#### PascalDragon

• Hero Member
• Posts: 3335
• 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

• Full Member
• Posts: 222
##### 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

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

• Full Member
• Posts: 132
##### 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;
88.     end.
89.
« Last Edit: August 04, 2021, 06:25:43 pm by BobDog »

#### avk

• Sr. Member
• Posts: 462
##### 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: 3335
• 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.