百万分之一斐波那契数 - Numpy Python 实施
Millionth Fibonacci Number - Numpy Python Implementation
我正在尝试实现一个用于检索百万分之一斐波那契数或更高数的代码。我将矩阵乘法与 Numpy 结合使用以加快计算速度。
根据我的理解，它应该花费 O(logN) 时间，并且一百万的最坏情况应该导致：将近 6 秒，这应该没问题。
以下是我的实现：
def fib(n):
import numpy as np
matrix = np.matrix([[1, 1], [1, 0]]) ** abs(n)
if n%2 == 0 and n < 0:
return -matrix[0,1]
return matrix[0, 1]
但是，留下 100 万，它甚至没有为 1000 生成正确的响应
实际回复：
43466557686937456435688527675040625802564660517371780402481729089536555417949051890403879840079255169295922593080322634775209689623239873322471161642996440906533187938298969649928516003704476137795166849228875
我的回复：
817770325994397771
为什么 Python 通常会从文档中截断响应，它甚至应该能够计算 10**1000 个值。我哪里错了？
Numpy 可以以高性能方式处理数字和计算（内存效率和计算时间）。因此，虽然 Python 可以处理足够大的数字，但 Numpy 不能。可以让Python进行计算得到结果，以性能降低作为交换。
示例代码：
import numpy as np
def fib(n):
# the difference is dtype=object, it will let python do the calculation
matrix = np.matrix([[1, 1], [1, 0]], dtype=object) ** abs(n)
if n%2 == 0 and n < 0:
return -matrix[0,1]
return matrix[0, 1]
print(fib(1000))
输出：
43466557686937456435688527675040625802564660517371780402481729089536555417949051890403879840079255169295922593080322634775209689623239873322471161642996440906533187938298969649928516003704476137795166849228875
PS：警告
Milionth Fibonacci 数非常大，您应该确保 Python 可以处理它。如果没有，您将不得不 implement/find 一些模块来处理这些大数字。
我不相信 numpy 在这里有多大帮助，因为它不直接支持 Python 在向量化操作中的非常大的整数。 O(logN) 算法的基本 Python 实现在我的笔记本电脑上在 0.15 秒内获得第 1 个斐波那契数。迭代（慢速）方法在 12 秒内完成。:
def slowfibo(N):
a = 0
b = 1
for _ in range(1,N): a,b = b,a+b
return a
# Nth Fibonacci number (exponential iterations) O(log(N)) time (N>=0)
def fastFibo(N):
a,b = 1,1
f0,f1 = 0,1
r,s = (1,1) if N&1 else (0,1)
N //=2
while N > 0:
a,b = f0*a+f1*b, f0*b+f1*(a+b)
f0,f1 = b-a,a
if N&1: r,s = f0*r+f1*s, f0*s+f1*(r+s)
N //= 2
return r
输出：
f1K = slowFibo(1000) # 0.00009 sec
f1K = fib(1000) # 0.00011 sec (tandat's)
f1K = fastFibo(1000) # 0.00002 sec
43466557686937456435688527675040625802564660517371780402481729089536555417949051890403879840079255169295922593080322634775209689623239873322471161642996440906533187938298969649928516003704476137795166849228875
f1M = slowFibo(1_000_000) # 12.52 sec
f1M = fib(1_000_000) # 0.2769 sec (tandat's)
f1M = fastFibo(1_000_000) # 0.14718 sec
19532821287077577316...68996526838242546875
len(str(f1M)) # 208988 digits
你的函数的核心是
np.matrix([[1, 1], [1, 0]]) ** abs(n)
在 Wiki 文章中讨论
np.matrix 将 ** 实现为 __pwr__，后者又使用 np.linalg.matrix_power。本质上，这是一个重复的 dot 矩阵乘法，通过将乘积按 2 的幂分组进行适度增强。
In [319]: M=np.matrix([[1, 1], [1, 0]])
In [320]: M**10
Out[320]:
matrix([[89, 55],
[55, 34]])
不鼓励使用np.matrix，所以我可以用
做同样的事情
In [321]: A = np.array(M)
In [322]: A
Out[322]:
array([[1, 1],
[1, 0]])
In [323]: np.linalg.matrix_power(A,10)
Out[323]:
array([[89, 55],
[55, 34]])
使用（新的）@ 矩阵乘法运算符，这与：
In [324]: A@A@A@A@A@A@A@A@A@A
Out[324]:
array([[89, 55],
[55, 34]])
matrix_power 做的事情更像是：
In [325]: A2=A@A; A4=A2@A2; A8=A4@A4; A8@A2
Out[325]:
array([[89, 55],
[55, 34]])
和一些比较时间：
In [326]: timeit np.linalg.matrix_power(A,10)
16.2 µs ± 58.1 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
In [327]: timeit M**10
33.5 µs ± 38.8 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
In [328]: timeit A@A@A@A@A@A@A@A@A@A
25.6 µs ± 914 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
In [329]: timeit A2=A@A; A4=A2@A2; A8=A4@A4; A8@A2
10.2 µs ± 97.9 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
numpy 整数在 int64、c 中实现，因此大小有限。因此我们得到适度的 100 溢出：
In [330]: np.linalg.matrix_power(A,100)
Out[330]:
array([[ 1298777728820984005, 3736710778780434371],
[ 3736710778780434371, -2437933049959450366]])
我们可以通过将 dtype 更改为 object 来解决这个问题。这些值然后是 Python 整数，并且可以无限增长：
In [331]: Ao = A.astype(object)
In [332]: Ao
Out[332]:
array([[1, 1],
[1, 0]], dtype=object)
幸运的是 matrix_power 可以干净地处理 object dtype:
In [333]: np.linalg.matrix_power(Ao,100)
Out[333]:
array([[573147844013817084101, 354224848179261915075],
[354224848179261915075, 218922995834555169026]], dtype=object)
通常对象 dtype 的数学运算速度较慢，但在这种情况下不会：
In [334]: timeit np.linalg.matrix_power(Ao,10)
14.9 µs ± 198 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
我猜这是因为数组的小 (2,2) 大小，快速编译的方法在这里没有用。这基本上是一个迭代任务，numpy 没有任何优势。
缩放还不错 - 将 n 增加 10，并且只增加 3-4 倍的时间。
In [337]: np.linalg.matrix_power(Ao,1000)
Out[337]:
array([[70330367711422815821835254877183549770181269836358732742604905087154537118196933579742249494562611733487750449241765991088186363265450223647106012053374121273867339111198139373125598767690091902245245323403501,
43466557686937456435688527675040625802564660517371780402481729089536555417949051890403879840079255169295922593080322634775209689623239873322471161642996440906533187938298969649928516003704476137795166849228875],
[43466557686937456435688527675040625802564660517371780402481729089536555417949051890403879840079255169295922593080322634775209689623239873322471161642996440906533187938298969649928516003704476137795166849228875,
26863810024485359386146727202142923967616609318986952340123175997617981700247881689338369654483356564191827856161443356312976673642210350324634850410377680367334151172899169723197082763985615764450078474174626]],
dtype=object)
In [338]: timeit np.linalg.matrix_power(Ao,1000)
53.8 µs ± 83 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
对象数据类型np.matrix:
In [340]: Mo = M.astype(object)
In [344]: timeit Mo**1000
86.1 µs ± 164 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
对于百万来说，时间并不像我预期的那么糟糕：
In [352]: timeit np.linalg.matrix_power(Ao,1_000_000)
423 ms ± 1.92 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
为了比较，我机器上的 fastFibo 时间是：
In [354]: fastFibo(100)
Out[354]: 354224848179261915075
In [355]: timeit fastFibo(100)
3.91 µs ± 154 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
In [356]: timeit fastFibo(1000)
9.37 µs ± 23.2 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
In [357]: timeit fastFibo(1_000_000)
226 ms ± 12.5 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
我正在尝试实现一个用于检索百万分之一斐波那契数或更高数的代码。我将矩阵乘法与 Numpy 结合使用以加快计算速度。
根据我的理解，它应该花费 O(logN) 时间，并且一百万的最坏情况应该导致：将近 6 秒，这应该没问题。
以下是我的实现：
def fib(n):
import numpy as np
matrix = np.matrix([[1, 1], [1, 0]]) ** abs(n)
if n%2 == 0 and n < 0:
return -matrix[0,1]
return matrix[0, 1]
但是，留下 100 万，它甚至没有为 1000 生成正确的响应
实际回复：
43466557686937456435688527675040625802564660517371780402481729089536555417949051890403879840079255169295922593080322634775209689623239873322471161642996440906533187938298969649928516003704476137795166849228875
我的回复：
817770325994397771
为什么 Python 通常会从文档中截断响应，它甚至应该能够计算 10**1000 个值。我哪里错了？
Numpy 可以以高性能方式处理数字和计算（内存效率和计算时间）。因此，虽然 Python 可以处理足够大的数字，但 Numpy 不能。可以让Python进行计算得到结果，以性能降低作为交换。
示例代码：
import numpy as np
def fib(n):
# the difference is dtype=object, it will let python do the calculation
matrix = np.matrix([[1, 1], [1, 0]], dtype=object) ** abs(n)
if n%2 == 0 and n < 0:
return -matrix[0,1]
return matrix[0, 1]
print(fib(1000))
输出：
43466557686937456435688527675040625802564660517371780402481729089536555417949051890403879840079255169295922593080322634775209689623239873322471161642996440906533187938298969649928516003704476137795166849228875
PS：警告
Milionth Fibonacci 数非常大，您应该确保 Python 可以处理它。如果没有，您将不得不 implement/find 一些模块来处理这些大数字。
我不相信 numpy 在这里有多大帮助，因为它不直接支持 Python 在向量化操作中的非常大的整数。 O(logN) 算法的基本 Python 实现在我的笔记本电脑上在 0.15 秒内获得第 1 个斐波那契数。迭代（慢速）方法在 12 秒内完成。:
def slowfibo(N):
a = 0
b = 1
for _ in range(1,N): a,b = b,a+b
return a
# Nth Fibonacci number (exponential iterations) O(log(N)) time (N>=0)
def fastFibo(N):
a,b = 1,1
f0,f1 = 0,1
r,s = (1,1) if N&1 else (0,1)
N //=2
while N > 0:
a,b = f0*a+f1*b, f0*b+f1*(a+b)
f0,f1 = b-a,a
if N&1: r,s = f0*r+f1*s, f0*s+f1*(r+s)
N //= 2
return r
输出：
f1K = slowFibo(1000) # 0.00009 sec
f1K = fib(1000) # 0.00011 sec (tandat's)
f1K = fastFibo(1000) # 0.00002 sec
43466557686937456435688527675040625802564660517371780402481729089536555417949051890403879840079255169295922593080322634775209689623239873322471161642996440906533187938298969649928516003704476137795166849228875
f1M = slowFibo(1_000_000) # 12.52 sec
f1M = fib(1_000_000) # 0.2769 sec (tandat's)
f1M = fastFibo(1_000_000) # 0.14718 sec
19532821287077577316...68996526838242546875
len(str(f1M)) # 208988 digits
你的函数的核心是
np.matrix([[1, 1], [1, 0]]) ** abs(n)
在 Wiki 文章中讨论
np.matrix 将 ** 实现为 __pwr__，后者又使用 np.linalg.matrix_power。本质上，这是一个重复的 dot 矩阵乘法，通过将乘积按 2 的幂分组进行适度增强。
In [319]: M=np.matrix([[1, 1], [1, 0]])
In [320]: M**10
Out[320]:
matrix([[89, 55],
[55, 34]])
不鼓励使用np.matrix，所以我可以用
In [321]: A = np.array(M)
In [322]: A
Out[322]:
array([[1, 1],
[1, 0]])
In [323]: np.linalg.matrix_power(A,10)
Out[323]:
array([[89, 55],
[55, 34]])
使用（新的）@ 矩阵乘法运算符，这与：
In [324]: A@A@A@A@A@A@A@A@A@A
Out[324]:
array([[89, 55],
[55, 34]])
matrix_power 做的事情更像是：
In [325]: A2=A@A; A4=A2@A2; A8=A4@A4; A8@A2
Out[325]:
array([[89, 55],
[55, 34]])
和一些比较时间：
In [326]: timeit np.linalg.matrix_power(A,10)
16.2 µs ± 58.1 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
In [327]: timeit M**10
33.5 µs ± 38.8 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
In [328]: timeit A@A@A@A@A@A@A@A@A@A
25.6 µs ± 914 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
In [329]: timeit A2=A@A; A4=A2@A2; A8=A4@A4; A8@A2
10.2 µs ± 97.9 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
numpy 整数在 int64、c 中实现，因此大小有限。因此我们得到适度的 100 溢出：
In [330]: np.linalg.matrix_power(A,100)
Out[330]:
array([[ 1298777728820984005, 3736710778780434371],
[ 3736710778780434371, -2437933049959450366]])
我们可以通过将 dtype 更改为 object 来解决这个问题。这些值然后是 Python 整数，并且可以无限增长：
In [331]: Ao = A.astype(object)
In [332]: Ao
Out[332]:
array([[1, 1],
[1, 0]], dtype=object)
幸运的是 matrix_power 可以干净地处理 object dtype:
In [333]: np.linalg.matrix_power(Ao,100)
Out[333]:
array([[573147844013817084101, 354224848179261915075],
[354224848179261915075, 218922995834555169026]], dtype=object)
通常对象 dtype 的数学运算速度较慢，但在这种情况下不会：
In [334]: timeit np.linalg.matrix_power(Ao,10)
14.9 µs ± 198 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
我猜这是因为数组的小 (2,2) 大小，快速编译的方法在这里没有用。这基本上是一个迭代任务，numpy 没有任何优势。
缩放还不错 - 将 n 增加 10，并且只增加 3-4 倍的时间。
In [337]: np.linalg.matrix_power(Ao,1000)
Out[337]:
array([[70330367711422815821835254877183549770181269836358732742604905087154537118196933579742249494562611733487750449241765991088186363265450223647106012053374121273867339111198139373125598767690091902245245323403501,
43466557686937456435688527675040625802564660517371780402481729089536555417949051890403879840079255169295922593080322634775209689623239873322471161642996440906533187938298969649928516003704476137795166849228875],
[43466557686937456435688527675040625802564660517371780402481729089536555417949051890403879840079255169295922593080322634775209689623239873322471161642996440906533187938298969649928516003704476137795166849228875,
26863810024485359386146727202142923967616609318986952340123175997617981700247881689338369654483356564191827856161443356312976673642210350324634850410377680367334151172899169723197082763985615764450078474174626]],
dtype=object)
In [338]: timeit np.linalg.matrix_power(Ao,1000)
53.8 µs ± 83 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
对象数据类型np.matrix:
In [340]: Mo = M.astype(object)
In [344]: timeit Mo**1000
86.1 µs ± 164 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
对于百万来说，时间并不像我预期的那么糟糕：
In [352]: timeit np.linalg.matrix_power(Ao,1_000_000)
423 ms ± 1.92 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
为了比较，我机器上的 fastFibo 时间是：
In [354]: fastFibo(100)
Out[354]: 354224848179261915075
In [355]: timeit fastFibo(100)
3.91 µs ± 154 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
In [356]: timeit fastFibo(1000)
9.37 µs ± 23.2 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
In [357]: timeit fastFibo(1_000_000)
226 ms ± 12.5 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)