arctangent函数的数值计算

Question

我很难解释arctangent函数的结果。对于我遇到的所有实现来说，这种行为都是一致的，所以这里我将限制自己使用NumPy和MATLAB。

这个想法是要有一个随机放置点的圆圈。其目的是表示它们在极坐标系统中的位置，由于它们都是均匀分布的，所以我期望atan2函数计算的θ角也会随机分布在区间-π.π上。

下面是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);

这里是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()

在这两种情况下，我都得到了如下结果，可以很容易地看到π/4的每一个倍数的“步骤”。

​

​

这看起来像是某种精确性的误差，但奇怪的是，在我预料不到的角度上。同样，这种行为也存在于普通的atan函数中。

Answer 1

注意，您使用的是整数。

因此，对于每一对(p，q)，都会有floor(sqrt(p**2 + q**2)/gcd(p,q)/r)对，它们的角度arctan(p,q)相同。对于(p，q)的倍数，gcd(p,q)是1

还请注意，p**2+q**2对于pi/2的倍数是1，对于pi/4的奇数倍数是2，因此我们可以预测，比pi/4的奇数倍数更多的项目是pi/4的偶数倍。这和我们在你的阴谋中看到的是一致的。

示例

让我们用整数坐标绘制半径为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}')

​

​

在这里，你可以看到，在图中，有几个点与其他点有相同的角度。45度有7点，90倍有10点，许多点有一个独特的角度。基本上，你有很多角度，很少的诗句和几个角度，击中许多点。

但总体上，这些点相对于角度几乎是均匀分布的。在这里，我画出了累积频率，它几乎是一条直线(如果分布是一一对应的话)，而bin频率形成了一些三角形的分形图案。

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)

​

​

如果圆圈的大小增加，并且你做的直方图足够宽，你不会注意到变化，否则你会在直方图中看到这个模式。

​

​