CODEWITHSUNDEEP
python
William Yi
William Yi

Reputation: 69

Fibonacci Sequence Mod 1000000007

Everyone knows that the Fibonacci sequence goes

F[0]=1, F[1]=1, F[2]=2, F[3]=3, F[4]=5, F[5]=8,

with F[n] = F[n-1]+F[n-2].

Now, how do you compute a number in the Fibonacci sequence when taken modulo 1000000007 = 10^9+7?

Needs to run as efficiently as possible, and in Python language :)

For example F[10**15] should take less than a second or so.

I know matrix exponentiation works, but how do you correct Matrix Exponentiation to reflect MODULO? (Another example, see https://www.nayuki.io/page/fast-fibonacci-algorithms)

Upvotes: 4

Views: 3398

Answers (1)

Him
Him

Reputation: 5549

The tricks needed:

1) Use the closed form of Fibonacci numbers, this is much faster than recursion. http://mathworld.wolfram.com/FibonacciNumber.html (formula 6)

2) modulo essentially factors over multiplication and addition, and sort of factors over division (you have to compute the multiplicative inverse in the mod space with the extended euclidean algorithm first) so you can basically just mod as you go. https://en.wikipedia.org/wiki/Modulo_operation#Equivalencies https://en.wikipedia.org/wiki/Extended_Euclidean_algorithm#Modular_integers

code:

def rootiply(a1,b1,a2,b2,c):
    ''' multipy a1+b1*sqrt(c) and a2+b2*sqrt(c)... return a,b'''
    return a1*a2 + b1*b2*c, a1*b2 + a2*b1

def rootipower(a,b,c,n):
    ''' raise a + b * sqrt(c) to the nth power... returns the new a,b and c of the result in the same format'''
    ar,br = 1,0
    while n != 0:
        if n%2:
            ar,br = rootiply(ar,br,a,b,c)
        a,b = rootiply(a,b,a,b,c)
        n /= 2
    return ar,br

def rootipowermod(a,b,c,k,n):
    ''' compute root powers, but modding as we go'''
    ar,br = 1,0
    while k != 0:
        if k%2:
            ar,br = rootiply(ar,br,a,b,c) 
            ar,br = ar%n,br%n
        a,b = rootiply(a,b,a,b,c)
        a,b = a%n, b%n
        k /= 2
    return ar,br

def fib(k):
    ''' the kth fibonacci number'''
    a1,b1 = rootipower(1,1,5,k)
    a2,b2 = rootipower(1,-1,5,k)
    a = a1-a2
    b = b1-b2
    a,b = rootiply(0,1,a,b,5)
    # b should be 0!
    assert b == 0
    return a/2**k/5

def powermod(a,k,n):
    ''' raise a**k, modding as we go by n'''
    r = 1
    while k!=0:
        if k%2:
            r = (a*r)%n
        a = (a**2)%n
        k/=2
    return r

def mod_inv(a,n):
    ''' compute the multiplicative inverse of a, mod n'''
    t,newt,r,newr = 0,1,n,a
    while newr != 0:
        quotient = r / newr
        t, newt = newt, t - quotient * newt
        r, newr = newr, r - quotient * newr
    if r > 1: return "a is not invertible"
    if t < 0: t = t + n
    return t

def fibmod(k,n):
    ''' compute the kth fibonacci number mod n, modding as we go for efficiency'''
    a1,b1 = rootipowermod(1,1,5,k,n)
    a2,b2 = rootipowermod(1,-1,5,k,n)
    a = a1-a2
    b = b1-b2
    a,b = rootiply(0,1,a,b,5)
    a,b = a%n,b%n
    assert b == 0
    return (a*mod_inv(5,n)*mod_inv(powermod(2,k,n),n))%n

if __name__ == "__main__":
    assert rootipower(1,2,3,3) == (37,30) # 1+2sqrt(3) **3 => 13 + 4sqrt(3) => 39 + 30sqrt(3)
    assert fib(10)==55
    #print fib(10**15)%(10**9+7) # takes forever because the integers involved are REALLY REALLY REALLY BIG
    print fibmod(10**15,10**9+7) # much faster because we never deal with integers bigger than 10**9+7

Upvotes: 2

Related Questions