4

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?

Mikey Freeman
  • 43
  • 5
  • 3
    google [modpow](https://stackoverflow.com/q/18577076/2521214) ... its much faster to do `modpow(a,b,c)` then `pow(a,b) mod c` because you do not need big numbers ... I saw some people (younger and on newer lanuages like JAVA or python) are using `powmod` name instead however IIRC it should be `modpow` So simply do power by squaring with modular arithmetrics ... just like I did in the linked `modpow` implementation – Spektre Dec 13 '20 at 09:36
  • @Spektre how about float numbers? in my case, (3 + sqrt(17)) ^ n. Is it also modpow? – Mikey Freeman Dec 13 '20 at 09:43
  • hmm that is a problem because it depends on what do you consider to be a mod of float? – Spektre Dec 13 '20 at 09:44
  • you can easily use [power by squaring or log/exp](https://stackoverflow.com/a/30962495/2521214) approach on floats and fixed point thge same way as on integers but mod is a special case ... so if you truncate your float to int and thed mod then yes if you use something like fmod then you need to multiply by number you are modulo to truncate to integer and then restore the result somehow ... – Spektre Dec 13 '20 at 09:46
  • @Spektre I mean, it doesn't need mod, but also I can't calculate very large n for (3 + sqrt(17)) ^ n :-?? – Mikey Freeman Dec 13 '20 at 09:46
  • Follow this. Hope it will be effective https://stackoverflow.com/a/538583/5929910 – mhhabib Dec 13 '20 at 09:47
  • @toRex Thanks I saw it before, and I use Python 3.8, but my numbers are very very large I think that it can't handle it. – Mikey Freeman Dec 13 '20 at 09:48
  • Perhaps you should explain a bit about the circumstances that lead to this formula, because it seems a bit odd: apparently you want to do modular arithmetic, but there is an (necessarily inexact) square root of 17 in here as well, and 17 is not a quadratic residue modulo 1000000007 – harold Dec 13 '20 at 09:49
  • @Spektre Thanks I will search about power by squaring :-? – Mikey Freeman Dec 13 '20 at 09:50
  • @harold it's about [this](https://math.stackexchange.com/a/686374/842598). (Sn) – Mikey Freeman Dec 13 '20 at 09:51
  • the problem is that power function implementation on floats is most likely using log/exp approach and computing that on arbitrary big floats is not easy as LUTs are out of question... if it is ok for the result try to truncate the `(3 + sqrt(17))` into integer `x` and then use `x**n` that should force to use integer power which should be fine ... however I am not coding in python so I do not know for sure how good it is implemented. You can also do do a fixed point computation if you need also few decimal places after dot ... – Spektre Dec 13 '20 at 09:51
  • @Spektre that would give me the wrong answer (I think) – Mikey Freeman Dec 13 '20 at 09:53
  • 1
    @Spektre it's about this https://math.stackexchange.com/a/686374/842598 (The Sn part) – Mikey Freeman Dec 13 '20 at 09:54
  • 1
    Use integer arithmetic throughout: you can compute `xd(n)` by finding the `n`th power of the 2-by-2 integer matrix `[[3, 2], [1, 0]]` and multiplying by the vector `[4, 1]` - that gives you `xd(n+1)` together with `xd(n)`. To compute the nth power of the matrix efficiently modulo `1000000007`, use the standard modular exponentiation algorithm (but applied to 2-by-2 integer matrices rather than plain integers). – Mark Dickinson Dec 13 '20 at 09:54
  • @MikeyFreeman OK I recommend using a different mathematical solution, for example based on a matrix power (this is probably possible but I didn't work it out) which would not involve square roots or divisions – harold Dec 13 '20 at 09:55
  • just try if it will overflow ... if not then do something like `a = (x*1000)**n; b =1000**n;` and then convert to float `y = a/b` which should give you your answer with limited precision but still better than integers alone ... the `1000` can be enlarged .... to improve accuracy – Spektre Dec 13 '20 at 09:55
  • @MarkDickinson Thanks, I think that would be the solution, I will work on it. – Mikey Freeman Dec 13 '20 at 09:57
  • @harold Yeah, I'll work on it. Thanks – Mikey Freeman Dec 13 '20 at 09:57
  • @Spektre I will try it, thank you – Mikey Freeman Dec 13 '20 at 09:57
  • @MarkDickinson Sorry, how did you come up with `[[3, 2], [1, 0]]` or `[4, 1]`? – Mikey Freeman Dec 13 '20 at 10:00
  • @MikeyFreeman: If you express `xd(n+1)` and `xd(n)` in terms of `xd(n)` and `xd(n-1)`, you get the matrix equation: `[xd(n+1), xd(n)] = [[3, 2], [1, 0]] * [xd(n), xd(n-1)]` (awkward to write in a comment - think of the vectors as columns vectors here). The top row of the matrix comes from the recurrence relation `xd(n+1) = 3*xd(n) + 2*xd(n-1)`; the bottom row just gives the identity `xd(n) = 1*xd(n) + 0*xd(n-1)`. The starting vector `[4, 1]` is `[xd(1), xd[0]]`. – Mark Dickinson Dec 13 '20 at 10:06
  • xd(n+1) = xd(n+2) + xd(n+1). This is the Sn (xd=Sn), not the An (https://math.stackexchange.com/a/686374/842598). Or am I wrong? – Mikey Freeman Dec 13 '20 at 10:16
  • @MarkDickinson No, I was wrong, Thanks, I appreciate your help – Mikey Freeman Dec 13 '20 at 10:21

3 Answers3

3

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

Thierry Lathuille
  • 23,663
  • 10
  • 44
  • 50
  • You are absolutely correct that the solution in terms of `a[n]` is much cleaner than dealing with floats. One thing to notice is that any linear recurrence admits a `O(log n)` solution. – user58697 Dec 13 '20 at 19:59
1

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.

rossum
  • 15,344
  • 1
  • 24
  • 38
Deepak Tatyaji Ahire
  • 4,883
  • 2
  • 13
  • 35
  • this will be good only for small `n` in `a^n` for big exponents you need power by squaring ... – Spektre Dec 13 '20 at 09:41
  • Correct @Spektre. – Deepak Tatyaji Ahire Dec 13 '20 at 09:41
  • 1
    However OP clarified that in `a^b` the `a` is `float` !!! so either fixed or floating point power is needed ... I would go for fixed as that one can be done on integers without probems with arbitrary precision log,exp which are in my opinion the cause of overflow – Spektre Dec 13 '20 at 09:58
  • 1
    @Spektre but OP is using a fundamentally wrong approach for the original problem, neither floating point power nor fixed point power would solve it – harold Dec 13 '20 at 10:05
1

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.

Thierry Lathuille
  • 23,663
  • 10
  • 44
  • 50