numpy.cos在某些数字上的工作时间明显更长

Question

TLDR：

numpy.cos()在特定数字上（例如24000.0），其工作时间延长了30％。加一个小增量（+0.01）可使numpy.cos()照常工作。

我不知道为什么。

---

在与一起工作时，我偶然发现了一个奇怪的问题numpy。我正在检查缓存工作，并意外制作了错误的图形- numpy.cos(X)时间取决于时间X。这是我修改后的代码（从我的Jupyter笔记本复制）：

import numpy as np
import timeit
st = 'import numpy as np'
cmp = []
cmp_list = []
left = 0
right = 50000
step = 1000
# Loop for additional average smoothing
for _ in range(10):
    cmp_list = []
    # Calculate np.cos depending on its argument
    for i in range(left, right, step):
        s=(timeit.timeit('np.cos({})'.format(i), number=15000, setup=st))
        cmp_list.append(int(s*1000)/1000)
    cmp.append(cmp_list)

# Calculate average times
av=[np.average([cmp[i][j] for i in range(len(cmp))]) for j in range(len(cmp[0]))]

# Draw the graph
import matplotlib.pyplot as plt
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
plt.plot(range(left, right, step), av, marker='.')
plt.show()

该图如下所示：

首先，我认为这只是一个随机故障。我重新计算了细胞，但结果几乎相同。因此，我开始使用step参数，包括计算数量和平均列表长度。但是，一切对这个数字都没有影响：

甚至更近：

在此之后，它range是无用的（它无法使用浮点数步进），因此我np.cos手动进行了计算：

print(timeit.timeit('np.cos({})'.format(24000.01),number=5000000,setup=st))
print(timeit.timeit('np.cos({})'.format(24000.00),number=5000000,setup=st))
print(timeit.timeit('np.cos({})'.format(23999.99),number=5000000,setup=st))

结果是：

3.4297256958670914
4.337243619374931
3.4064380447380245

np.cos()比24000.01长30％时精确计算24000.00！

还有另外一个奇怪的数字（大约在500000左右，我不记得了）。

我浏览了numpy文档及其源代码，但与这种效果无关。我知道三角函数使用几种算法，具体取决于值的大小，精度等，但是对我来说，确切的数字可以更长的计算时间感到困惑。

为什么np.cos()会有这种奇怪的效果？这是某种处理器的副作用numpy.cos吗（因为使用依赖于处理器的C函数）？如果有帮助，我已经安装了Intel Core i5和Ubuntu。

---

编辑1：我尝试使用AMD Ryzen 5在另一台计算机上重现它。结果只是不稳定。这是相同代码的三个顺序运行的图形：

import numpy as np
import timeit

s = 'import numpy as np'
times = []
x_ranges = np.arange(23999, 24001, 0.01)
for x in x_ranges:
    times.append(timeit.timeit('np.cos({})'.format(x), number=100000, setup=s))

# ---------------

import matplotlib.pyplot as plt

fig = plt.figure()
ax = fig.add_subplot(111)
plt.plot(x_ranges, times)
plt.show()

嗯，有一些模式（例如大部分一致的左侧部分和不一致的右侧部分），但是它与Intel处理器的运行有很大不同。看起来这实际上只是处理器的特殊方面，而AMD的不确定性行为更容易预测：)

PS @WarrenWeckesser感谢您的``np.arange`''功能。它确实很有用，但它不会像预期的那样改变结果。

Answer 1

对于那些特殊数字，计算结果的缓慢可能与精确舍入和制表者的困境有关。

为说明起见，假设您将一个指数函数表制作到4个地方。然后exp（1.626）= 5.0835。应该四舍五入为5.083还是5.084？如果更仔细地计算exp（1.626），它将变为5.08350。然后是5.083500。然后5.0835000。由于exp是先验的，因此在区分exp（1.626）是5.083500 ... 0ddd还是5.0834999 ... 9ddd之前，它可以任意进行很长时间。

但是，由于这个原因，IEEE标准不要求先验函数要精确四舍五入，但是该math.cos函数的实现可能会受到此问题的困扰，同时会尽力计算出最准确的结果，然后找出效果不值得的努力。

为了证明某些数字是这种情况，X必须math.cos(X)高精度地计算的值并检查其二进制表示-尾数的可表示部分必须跟随以下模式之一：

• 1和长期的0
• 0和长期的1（当计算值的精度低于适应运行中所有1的精度时，此情况显示为第一个）

因此，数字将成为超越函数的慢变量的概率为1/2 ⁿ，其中n是算法看到的上述模式的最大长度，在此之后，它放弃尝试得出精确的舍入结果。

演示突出显示了IEEE 754双精度情况下尾数的可表示部分（尾数有53位）：

In [1]: from mpmath import mp

In [2]: import math

In [3]: def show_mantissa_bits(x, n, k):
   ...:     print(bin(int(mp.floor(abs(x) * 2**n)))[2:])
   ...:     print('^'*k)
   ...:     

In [4]: mp.prec = 100

In [5]: show_mantissa_bits(mp.cos(108), 64, 53)
110000000100001011001011010000110111110010100011000011000000000
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

In [6]: show_mantissa_bits(mp.cos(108.01), 64, 53)
101110111000000110001101110001000010100111000010101100000100110
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

In [7]: show_mantissa_bits(mp.cos(448), 64, 53)
101000101000100111000010111100001011111000001111110001000000000
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

In [8]: show_mantissa_bits(mp.cos(448.01), 64, 53)
101001110110001010010100100000110001111100000001101110111010111
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

In [9]: show_mantissa_bits(mp.cos(495), 64, 53)
11001010100101110110001100110101010011110010000000000011111111
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

In [10]: show_mantissa_bits(mp.cos(495.01), 64, 53)
11010100100111100110000000011000110000001001101100010000001010
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

In [11]: show_mantissa_bits(mp.cos(24000), 64, 53)
11001000100000001100110111011101001101101101000000110011111111
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

In [12]: show_mantissa_bits(mp.cos(24000.01), 64, 53)
10111110011100111001010101100101110001011010101011001010110011
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^