Расчет очень большой мощности

Глядя на ответ на 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) занимает всего в 10 раз больше времени, чем s(10**6).

Поскольку значение 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 в конце, мы не можем брать по модулю на каждом шаге и держать числа, которыми мы манипулируем, маленькими.

Однако использование предложенного подхода с умножением матриц (я не видел ссылки на формулу рекурсии, когда начинал с этого ответа, очень плохо!) Позволит вам это сделать.