Fastest Way to calculate Nth fibonacci number modulo M

general
, , , , ,
1

1 ≤ N,M ≤ 1018
code using Fast Doubling

#include 
#include 
#include 
using namespace std;
#define MAXN 60 
#define MAXM 4
long long int M,N ; 
long long F[MAXN][MAXM],FB[MAXN+5];
void readint(long long &a)
{
register int c;
a = 0;
do c = getchar_unlocked(); while(c < '0');
do{
a = (a << 1) + (a <= '0');
}
void printint(long long int a)
{
char s[21];
int t = -1;
 
do{
s[++t] = a % 10 + '0';
a /= 10;
}while(a > 0);
 
while(t >= 0)putchar_unlocked(s[t--]);
putchar_unlocked('\n');
}
long long int mulmod(long long int a,long long int b) {
   long double res = a;
   res *= b;
   long long int c = (long long)(res / M);
   a *= b;
   a -= c * M;
   a %= M;
   if (a < 0) a += M;
   return a;
}
long long int  f(long long int  n,int depth = 0 ) {
	if(n<=MAXN){
		return FB[n]%M ;
	}
	int val = n%4 ;	
	if (F[depth][val] != -1) 
		return F[depth][val];
	long long int  k=n >> 1;
	long long int a,b,c;	
	a = f(k,depth+1) ;
	b = f(k-1,depth+1) ;
	if (n%2==0) { 
		F[depth][val] = (mulmod(a,a) + mulmod(b,b));
	}else { 
		c = f(k+1,depth+1) ;		
		F[depth][val] = (mulmod(a,c) + mulmod(b,a));
	}
	if(F[depth][val] >= M)
		F[depth][val] -= M ;
	return F[depth][val] ;
}
 
int main(){
	int T = 63450 ;
	FB[0] = 1 ;
	FB[1] = 1 ;
	for(int i=2;i<=64;i++){
		FB[i] = FB[i-1] + FB[i-2] ;
	}
	while(T--){
		memset(F,-1,sizeof F) ;
		readint(N) ;
		readint(M) ;
		printint((N==0 ? 0 : f(N-1))) ;
	}	
	return 0 ; 
}
3 Likes
2

Also, I’d like to bring upon this method also:

MATRIX EXPONENTIATION:
(n+1)th fibonacci number is obtained by multiplying the matrix M n times where M is a 2X2 matrix defined as:

M=  [1,1]
    [1,0]

M*Fk= Fk+1

The matrix obtained by multiplying this matrix M n times will give us the n+1th fibonacci number as the entro [0,0] (row=0 and colomn=0). Why such happens can be seen this way:

| 1   1 | x |  f(n)  | = | f(n+1) | // this will represent the current solution
| ----- |   | f(n-1) |   | f(n)   | 
 /|\             /|\
  |               |

M matrix Solution to the previous state

f(n+1) is f(n) + f(n-1) which can come if M[0][0]=1 and M0=1(do matrix multiplication and u’ll get this). f(n) is possible if, M[1][0]=1 and M1=0, that gives us the value of M[][].

Multiplying with M on both sides above, we get:

M x M x |  f(n)  | = M x | f(n+1) | = | f(n+2) |
        | f(n-1) |       |  f(n)  |   | f(n+1) |

Again multiplying repeatedly will gives us :

M^k x |  f(n)  | = | f(n+k) |
      | f(n-1) |   |f(n+k-1)|

From this, we see that a matrix of the first three fibonacci numbers, times a matrix of any three consecutive fibonacci numbers, equals the next. So we know that:

     n
[1 1]   =  [fib(n+1) fib(n)  ]
[1 0]      [fib(n)   fib(n-1)]

which makes the complexity of computation as O(logn) if computed efficiently. For code, you may like seeing this.

Also, wiki says a closed form solution of the same is :

The nth Fibonacci number is given by

       f(n) = Floor(phi^n / sqrt(5) + 1/2)
         where
       phi = (1 + sqrt(5)) / 2

Even this can be computed in O(logn) time. The result is deducible from the very same fact of the above exponentiation method. If the result of the matrix is calculating by finding eigen values, it’ll come this only.

4 Likes
3

@damn_me : thnx damn_me
The method u have suggested is the fastest mehtod known uptill now for calculating fibonacci numbers.

1 Like
4

Is this the fastest?
How many cases does your code pass here?

I got 61400.

3 Likes
5

Thanks a lot.But pls pls explain a bit too, what u did with the “val” in F[depth][val]???
And in - “if (a < 0) a += M;return a;”, when a is going to be negative?
Thanks in advance.

6

could anyone explain what is the logic behind the function long long int f(long long int n, int depth=0) and
long long int mulmod(long long int a, long long int b) . Thankyou in advance.

7 
https://www.codechef.com/problems/CURR3
#include
using namespace std;
#define ll long long 
#define m 1000000007
void multiply(ll f[][2],ll g[][2])
{
    ll a=(f[0][0]*g[0][0]+f[0][1]*g[1][0])%m;
    ll b=(f[0][0]*g[0][1]+f[0][1]*g[1][1])%m;
    ll c=(f[1][0]*g[0][0]+f[1][1]*g[1][0])%m;
    ll d=(f[1][0]*g[0][1]+f[1][1]*g[1][1])%m;
    
    f[0][0]=a;
    f[0][1]=b;
    f[1][0]=c;
    f[1][1]=d;
}
void power(ll f[2][2],ll n)
{
    ll g[2][2]={{1,1},{1,0}};
    if(n==0||n==1)
    return;
    power(f,n/2);
    multiply(f,f);
    
    if(n%2==1)
    multiply(f,g);
}
ll fib(ll n)
{
    ll f[2][2]={{1,1},{1,0}};
    if(n==0)
    return 0;
    power(f,n-1);
    return f[0][0]%m;
}
 
int main() {
	ll n;
	cin>>n;
        cout<<fib(n)%m;
}
1 Like
8

You can add related questions to it, I remember there was a question on spoj based on some GCD thing whose logic was just finding the n’t fibonacci number. Range was large, don’t remember the name or code of that.

9

Method suggested in the thread is even faster than Matrix Exponentiation … FAST DOUBLING

1 Like
10

Fast doubling is the fastest known method, although it is derived from the matrix exponentiation method but many redundant calculations are not there in fast doubling. And thus it differs from the matrix one by a significant constant factor. I mentioned there in the post along with because I myself learnt it the hard way the time I did :stuck_out_tongue:

1 Like
11

@ma5termind I think its better if you change the tag of your thread, include fibonacci also. Search query of People unaware of the method will be “fibonacci”, “large number” or similar. So, its just a suggestion.

12

@damm_me sure i will …

13

Have not tried yet … you can check if you want … and do inform about the result here …

1 Like
14

WA on spoj… LOL
Try it yourself!

15

how did u get 61400?,I got only 30k :/.

1 Like
16

@kislaya- Like in matrix exponentiation, we break n to n/2 size.
Similarly there is a way to break it into n/3.
So my code runs in log (base 3). Thats why its faster. But I wanted to know if anyone knows what algorithm passes 500000 test cases. 2 people have achieved that in the link.

1 Like
17

change limit mentioned above from 55 to 60 or 65 you will get AC on spoj and also consider the case when N = 0 … 60,000 are easily achievable using this code …

1 Like
18

@i_am_what_i_am

can you share your code. we are more happy in learning something new … :smiley:

2 Likes
19

Hmm… Does anyone have any idea how people cam pass 5,00,000 cases?
I have no idea…

20

also remove that printf statement … to get it accepted on spoj …

1 Like