计算非常大的功率

Calculating very large power

我想计算非常大的数字,例如 n = 10^15。

不知何故我不能,因为 OverflowError。

xd = lambda n : ((((5+ sqrt(17)) * ((3 + sqrt(17)) ** n)) - ((5-sqrt(17))* ((3 - sqrt(17)) ** n)))/((2 ** (n+1)) * sqrt(17)))

即使n=1000,也不会计算

不过,我应该提到我想要它的模块化 (1000000007)

解决方案是什么?

由于n的值很高,整数溢出很明显

模运算遵循以下规则:

  1. 加法: (a+b)%m = (a%m + b%m)%m
  2. 减法: (a-b)%m = (a%m + b%m + m)%m
  3. 乘法: (a*b)%m = (a%m * b%m)%m
  4. 指数:使用循环。

示例:对于a^n,使用a = (a%m * a%m)%m,n次。

对于较大的 n 值,使用 python 的 pow(x, e, m) 函数计算取模时间少了很多。

仅使用整数计算它的一种有效方法是使用二项式展开来展开表达式。稍微重新排列后,我们得到了一种相当简单的计算方法,偶数和奇数幂项的公式几乎相同:

def form(n, mod):
    cnk = 1
    total = 0
    for k in range(n+1):
        term = cnk * 3**k * 17**((n-k)//2)
        if (n-k) % 2 == 1:
            term *= 5
        total += term
        cnk *= (n-k)
        cnk //= (k+1)

        
    return (total // (2**n)) #% mod

我们可以将其与您的原始公式进行比较以检查结果:

from math import sqrt

def orig(n):
    return ((((5+ sqrt(17)) * ((3 + sqrt(17)) ** n)) - ((5-sqrt(17))* ((3 - sqrt(17)) ** n)))/((2 ** (n+1)) * sqrt(17)))

for n in range(20):
    print(n, orig(n), form(n, mod))

输出:

0 1.0000000000000002 1
1 4.0 4
2 14.000000000000002 14
3 50.0 50
4 178.0 178
5 634.0000000000001 634
6 2258.0 2258
7 8042.0 8042
8 28642.000000000004 28642
9 102010.00000000001 102010
10 363314.0 363314
11 1293962.0000000002 1293962
12 4608514.0 4608514
13 16413466.000000004 16413466
14 58457426.00000001 58457426
15 208199210.00000003 208199210
16 741512482.0000001 741512482
17 2640935866.000001 2640935866
18 9405832562.0 9405832562
19 33499369418.000004 33499369418

对于不太大的 n 值来说“相当”快(在旧机器上测试):

#%timeit form(1000, mod)
# 9.34 ms ± 87.6 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

#%timeit form(10000, mod)
# 3.79 s ± 14.8 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

#%timeit form(20000, mod)
# 23.6 s ± 37.3 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

最后一次测试,在取模之前,我们有一个 11033 位数字。

这种方法的主要问题是,由于我们必须在最后除以 2**n,我们不能在每一步都取模并保持我们操作的数字很小。

使用矩阵乘法的建议方法(当我开始这个答案时我没有看到 link 到递归公式,太糟糕了!)将允许你这样做,但是。

查看 the answer on maths.stackexchange 公式的来源,似乎最容易计算的是 a(n)。

所以,这可以很简单地通过递归来计算,而这次,由于我们只使用乘法和加法,我们可以利用模运算的规则,并使我们操作的数字保持较小:

def s(n, mod):
    a1 = 1
    a2 = 3
    for k in range(n-1):
        a1, a2 = a2, (3*a2 + 2* a1) % mod
    return (a1 + a2) % mod


mod = 1000000007

print(s(10, mod))
# 363314, as with the other formulas...

print(s(10**6, mod))
# 982192189

%timeit s(10**6, mod)
# 310 ms ± 6.46 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

%timeit s(10**7, mod)
# 3.39 s ± 93.8 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

我们得到与其他公式相同的结果(这真是一件好事...)。由于计算中使用的数字保持相同大小,最多模数的5倍,计算时间约为O(n) - s(10**7)只比s(10**6).[=14多10倍的时间=]