Numerical evaluation of arctangent functions

Question

I am having difficulties to interpret results of arctangent functions. This behaviour is consistent for all implementations I came across, so I will here limit myself to NumPy and MATLAB.

The idea is to have circle of randomly placed points. The goal is to represent their positions in polar coordinate system and since they are uniformly distributed, I expect the θ angle (which is calculated using atan2 function) to be also distributed randomly over interval -π ... π.

Here is the code for MATLAB:

stp = 2*pi/2^8;
siz = 100;
num = 100000000;

x = randi([-siz, siz], [1, num]);
y = randi([-siz, siz], [1, num]);
m = (x.^2+y.^2) < siz^2;

[t, ~] = cart2pol(x(m), y(m));
figure()
histogram(t, -pi:stp:pi);

And here for Python & NumPy:

import numpy as np
import matplotlib.pyplot as pl

siz = 100
num = 100000000

rng = np.random.default_rng()
x = rng.integers(low=-siz, high=siz, size=num, endpoint=True)
y = rng.integers(low=-siz, high=siz, size=num, endpoint=True)
m = (x**2+y**2) < siz**2

t = np.arctan2(y[m], x[m]);
pl.hist(t, range=[-np.pi, np.pi], bins=2**8)
pl.show()

In both cases I got results looking like this, where one can easily see "steps" for each multiple of π/4.

It looks like some sort of precision error, but strangely for angles where I would not expect that. Also this behaviour is present for ordinary atan function as well.

As Bob said below, it's a discretization issue. You can create non-discrete samples with `x = (rand(1, num)-0.5)*2*siz;`, or you could set the origin of your domain at a random location: `x = x+rand; y = y+rand`. — Cris Luengo, Mar 02 '22 at 20:47

Bob · Answer 1 · 2022-03-03T14:50:58.917

Notice that you are using integers

So for each pair (p,q) you will have floor(sqrt(p**2 + q**2)/gcd(p,q)/r) pairs that give the same angle arctan(p,q). Then for the multiples of (p,q) the gcd(p,q) is 1

Notice also that p**2+q**2 is 1 for the multiples of pi/2 and 2 for the odd multiples of pi/4, with this we can predict that there will be more items that are even multiples of pi/4 than odd mulitples of pi/4. And this agrees with what we see in your plot.

Example

Let's plot the points with integer coordinates that lie in a circle of radius 10.

import numpy as np
import matplotlib.pyplot as plt
from collections import Counter
def gcd(a,b):
    if a == 0 or b == 0:
        return max(a,b)
    while b != 0:
        a,b = b, a%b
    return a;
R = 10
x,y = np.indices((R+1, R+1))
m = (x**2 + y**2) <= R**2
x,y = x[m], y[m]
t = np.linspace(0, np.pi / 2)
plt.figure(figsize=(6, 6))
plt.plot(x, y, 'o')
plt.plot(R * np.cos(t), R * np.sin(t))

lines = Counter((xi / gcd(xi,yi), 
                 yi / gcd(xi,yi)) for xi, yi in zip(x,y))
plt.axis('off')
for (x,y),f in lines.items():
    if f != 1:
        r = np.sqrt(x**2 + y**2)
        plt.plot([0, R*x/r], [0, R*y/r], alpha=0.25)
        plt.text(R*1.03*x/r, R*1.03*y/r, f'{int(y)}/{int(x)}: {f}')

Here you see on the plot a few points that share the same angle as some other. For the 45 degrees there are 7 points, and for multiples of 90 there are 10. Many of the points have a unique angle. Basically you have many angles with few poitns and a few angles that hit many points.

But overall the points are distributed nearly uniformly with respect to angle. Here I plot the cumulative frequency that is nearly a straight line (what it would be if the distribution was unifrom), and the bin frequency form some triangular fractal pattern.

R = 20
x,y = np.indices((R+1, R+1))
m = (x**2 + y**2) <= R**2
x,y = x[m], y[m]
plt.figure(figsize=(6,6))
plt.subplot(211)
plt.plot(np.sort(np.arctan2(x,y))*180/np.pi, np.arange(len(x)), '.', markersize=1)

plt.subplot(212)
plt.plot(np.arctan2(x,y)*180/np.pi, np.gcd(x,y), '.', markersize=4)

If the size of the circle increases and you do a histogram with sufficiently wide bins you will not notice the variations, otherwise you will see this pattern in the histogram.

It looks you are right, thank you! At least if I try floats, it gives me expected results. However, it is difficult for me to follow your reasoning. Would you be so kind and expand it little bit? On what interval is the ```[p, q]``` pair defined (from second paragraph it looks like it must sits on square edge [-1..1;-1..1])? Isn't ```r``` equal to ```sqrt(p^2+q^2)```? How can ```p^2+q^2``` gives you different results for odd and even multiples of ```pi/4``` since squares of ```p``` ang ```q``` must be equal in all cases (for ```pi/4```)? Thank you. — rad, Mar 03 '22 at 10:36

Numerical evaluation of arctangent functions

1 Answers1

Example