Calculating very large power

Question

I want to calculate for very large numbers like n = 10^15.

Somehow I can't, because of OverflowError.

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

even for n=1000, it wouldn't be calculated.

Though, I should mention that I want the modular of it (1000000007)

What would be the solution?

Answer 1

Looking at the answer on maths.stackexchange where the formula comes from, it appears that the easiest thing to calculate are the a(n).

So, this can be calculated by recurence very simply, and this time, as we only use multiplications and additions, we can take advantage of the rules of modulo arithmetic and keep the numbers we manipulate small:

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)

We get the same results as with the other formulas, (which is a really good thing...). As the numbers used during the calculation keep the same size, at most 5 times the modulo, the calculation time is about O(n) - s(10**7) takes only 10 times more time than s(10**6).

Answer 2

A working way to calculate it with integers only is to develop your expression using the binomial expansion. After rearranging it a bit, we get a rather easy way to calculate it, with an almost identical formula for terms of even and odd power:

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

We can compare it to your original formula to check the results:

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))

Output:

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

It is "rather" fast for not to large values of n (tested on an old machine):

#%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)

For the last test, before taking the modulo, we have a 11033 digits number.

The main problem with this approach is that, as we have to divide by 2**n at the end, we can't take the modulo at each step and keep the numbers we manipulate small.

Using the suggested approach with matrix multiplication (I hadn't seen the link to the recursion formula when I started with this answer, too bad!) will allow you to do this, though.

Answer 3

As the value of n is very high, integer overflow is obvious.

Follow the following rules for modular arithmetic:

1. Addition: (a+b)%m = (a%m + b%m)%m
2. Subtraction: (a-b)%m = (a%m + b%m + m)%m
3. Multiplication: (a*b)%m = (a%m * b%m)%m
4. Exponential: Use loop.

Example: For a^n, use a = (a%m * a%m)%m, n number of times.

For larger values of n, use the python's pow(x, e, m) function to get the modulo calculated which takes a lot less time.