计算非常大的指数

Question

我想计算像 n = 10^15 这样非常大的数字。

但是由于溢出错误，我无法进行计算。

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。

那么，应该怎么解决呢？

- Mikey Freeman

---

使用 modpow(a,b,c) 比起 pow(a,b) mod c 更快，因为你不需要使用大数。我看到一些人（年轻人或使用新语言如 JAVA 或 Python）使用 powmod 名称来代替，但据我所知应该是 modpow。所以就像链接中的 modpow 实现一样，使用平方法和模算术即可进行幂运算。请注意，不要改变原来的意思。 - Spektre

@Spektre，浮点数怎么办？在我的情况下，(3 + sqrt(17)) ^ n。这也是modpow吗？ - Mikey Freeman

嗯，这是一个问题，因为它取决于您认为浮点数的模是什么？ - Spektre

@Spektre，这是关于这个的 https://math.stackexchange.com/a/686374/842598（Sn部分） - Mikey Freeman

始终使用整数算术：您可以通过找到2x2整数矩阵[[3,2]，[1,0]]的第n个幂并乘以向量[4,1]来计算xd（n）-这会给您xd（n + 1）和xd（n）。为了高效地计算模1000000007的矩阵的第n个幂，请使用标准模指数算法（但应用于2x2整数矩阵而不是普通整数）。 - Mark Dickinson

Answer 1

查看数学堆栈交换上的答案，可以得知计算最简单的是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) 多10倍。

- Thierry Lathuille

---

你说的没错，用 a[n] 的解决方案比使用浮点数更加简洁。需要注意的一点是，任何线性递归都有一个 O(log n) 的解决方案。 - user58697

Answer 2

由于 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) 函数来计算模数，这需要更少的时间。

- Deepak Tatyaji Ahire

---

对于 a^n 中的较小指数 n，简单的乘法即可计算幂，但对于大幂次，则需要采用“平方法”（power by squaring）进行求解。 - Spektre

正确的@Spektre。 - Deepak Tatyaji Ahire

然而，OP澄清了a^b中的a是float！因此需要固定或浮点幂...我会选择固定幂，因为可以在整数上进行计算，而不会出现任意精度log、exp引起溢出的问题。 - Spektre

@Spektre 但是原问题的解决方法是基本错误的，无论是浮点幂还是定点幂都无法解决它。 - harold

Answer 3

用整数计算的可行方法是使用二项式展开式来发展您的表达式。稍微重新排列一下后，我们可以得到一种相当简单的计算方法，对于偶数和奇数幂的项几乎有相同的公式:

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，因此我们无法在每个步骤中取模并保持我们操作的数字较小。

使用矩阵乘法的建议方法（当我开始回答时，我还没有看到递归公式的链接，太糟糕了！）将允许您这样做。

- Thierry Lathuille

---

谢谢，我感激您的时间和专注。 - Mikey Freeman

我将选择矩阵乘法。 - Mikey Freeman