How exactly is atan2 implemented in matlab?

Hello Dear Matlab Community,

I am currently experimenting with the atan2 function built into matlab. According to my internet research there are several ways to implement atan2, see e.g. Wikipedia (https://en.wikipedia.org/wiki/Atan2).

My question now is, how is atan2 implemented in Matlab?

Thank you very much for your help!

[Bjorn Gustavsson]

What do you mean by "implemented"? The numerical algorithm to calculate the angle or the order of the arguments? If the latter you should find the information through:

helpwin atan2

Ig the former you might find some information in the documentation but I remember Cleve Moler bloging about the algorithms for the tan (and atan?) functions a couple of years ago, perhaps you can find it if you search the mathworks-blog-archives?

HTH

[Adam Danz]

Check out the 'More info' section

https://www.mathworks.com/help/matlab/ref/atan2.html#bucvbr8-4

[Stefan Holzinger]

@Bjorn Gustavsson: With "implementation" i mean the numerical algorithm to calculate the angle. I haven't find the blog of Cleve Moler yet, but im goint to search for it. Thank for this hint.

@Adam Danz: I have already looked at the online help. The reference to the IEEE standard is interesting. I will try to replicate atan2 using the IEEE standard.

Answer 1

[Jan]

Which problem do you want to solve?

The easiest way to see the source code to create the atan2 function is to apply for a job at MathWorks and sign the NDA. If you job is concerned with these codes, you can take a look in it. But does this really help to solve your problem?

Some functions of Matlab are calcualted by code from fdlibm library, as far as I remember e.g. acos (this was mentioned in the documentation of older releases), see http://www.netlib.org/ (link). In my tests the output of Matlab's atan2 and e_atan2.c have not been identical for some inputs, but it depends on the compiler also.

If you just want to know, how atan2 can be implemented efficiently, look at the netlib.org source code. The code used in Matlab might differ a little bit.

[Stefan Holzinger]

I'm working on a research project. However, the atan2 function plays an essential role there.

I'll take a look at the code on netlib, maybe I can use it to emulate the exact behavior of the Matlab atan2 implementation.

In my previous attempts I was able to simulate all results of the Matlab atan2 implementation. But one feature seems to be missing in my implementation. The Matlab atan2 implementation provides atan2(sin(-pi), cos(-pi))=-pi for the case. Do you know how Matlab calculates this result? In my implementation I get MyAtan2(sin(-pi), cos(-pi))=+pi.

Here is the code:

function z = MyAtan2(y, x)
if x>0
    z=atan(y/x);
end
if y>=0 && x<0
    z=pi+atan(y/x);
end
if y<0 && x<0
    z=-pi+atan(y/x);
end
if y>0 && x==0
    z=pi/2;
end
if y<0 && x==0
    z=-pi/2;
end
end

[Jan]

@Stefan: I've formatted your code to improve the readability. You can use the tools for formatting also.

Of course I do not know, how atan2 is implemented in Matlab - the licence agreement prohibits a reverse engineering. But I'm convinced that the code used in Matlab is ways more complicated, catchs NaN, Inf, inputs of different classes and arrays of different sizes. The difference to the IEEE suggestions is mentioned in the documentation: Because Matlab treats -0 as 0 consequently, atan2(0,-0) replies 0 instead of pi.

Your code can be accelerated by using elseif 's. After the case x>0 has been caugth already, there is no need to check the other cases:

function z = MyAtan2(y, x)
if x>0
    z = atan(y / x);
elseif y>=0 && x<0
    z = pi + atan(y / x);
elseif y<0 && x<0
    z = -pi + atan(y / x);
elseif y>0 && x==0
    z = pi / 2;
elseif y<0 && x==0
    z = -pi / 2;
end
end

The case x==0, y==0 is not caught. This is a typical programming problem, if several conditions are tested. The general rule is: "No if without an else!" Split the conditions to be sure, that all possibilities are considered:

if x > 0
    z = atan(y / x);
elseif x < 0  % Care about y in a 2nd level
    ...
else  % x == 0 or NaN:
    z = sign(y) * pi / 2;  % Is the correct and exhaustive?
end

Now all cases of x are considered by design.

Note that sin(-pi) is not 0.0, but -1.22e-16. Such rounding effects of the input must be taken into account carefully. Look at the netlib implementation for an exhaustive approach.