coj 1053 寂寞的素数 （费马小定理 + 矩阵求斐波那契）

Description

Fibonacci 数列是大家的老朋友了 
Fib[1]=1, Fib[2]=1, Fib[3]=2……Fib[n]= Fib[n-1]+ Fib[n-2] 
p是一个不甘寂寞的素数。某日它找到Fibonacci 数列，跳到Fibonacci 数列头上玩了起一个传说中的游戏。于是我们看到： 
Fib[a]^p+ Fib[a+1]^p+ Fib[a+2]^p……….. + Fib[b]^p 
好了，我们的问题来了，求上述表达式的值(即fab[i]的p次方之和，b>=i>=a)，结果模p。

Input

输入三个整数a、b、p，满足： 
1000,000,000>a>0 
1000,000,000>b>a 
1000,000,000 > p

Output

一个整数。题目中要求表达式 模p 的结果

Sample Input

2    4    3

Sample Output

0

一：费马小定理是数论中的一个定理：假如a是一个整数，p是一个质数，那么

当是素数的时候，，所以欧拉定理变为：

   或   

两个整数，，若它们除以正整数所得的余数相等，则称，对于模同余

记作

二： 欧拉定理 （也称 费马-欧拉定理 或 欧拉函数定理 ）是一个关于同余的性质。欧拉定理表明，若 为正 整数 ，且 互素 ， ，则

其中为欧拉函数，为同余关系。

三：斐波那契数列求和 .f(0)+f⑴+f⑵+…+f(n)=f(n+2)-1。

所以：（Fib[a]^p+ Fib[a+1]^p+ Fib[a+2]^p……….. + Fib[b]^p）% p

就是求    (Fib[a]%p+ Fib[a+1]%p+ Fib[a+2]%p……….. + Fib[b]%p )%p

而 Fib[n] %p = (Fib[n-1]%p + Fib[n-2]%p) %p可以通过矩阵相乘来求。

要求a ~ b 的和，可以有斐波那契求和公式先求出 Sum(Fib[0]~Fib[a-1]) 和 Sum(Fib[b]) ,然后用 Sum( Fib[ b ] ) - Sum( Fib[a] ) 得到a ~ b的和

但是 Sum( Fib[ b ] )%p - Sum( Fib[a] )%p的结果可能会小于p所以 (Sum( Fib[ b ] )%p - Sum( Fib[a] ) %p + p )%p 为最后结果。

代码：

#include<iostream>
#include<string.h>
using namespace std;

struct matrics
{
   __int64 a[2][2];
};

matrics mul(matrics x , matrics y , int p)
{
    matrics ans;
    ans.a[0][0]=(x.a[0][0]*y.a[0][0]+x.a[0][1]*y.a[1][0])%p;
	ans.a[0][1]=(x.a[0][0]*y.a[0][1]+x.a[0][1]*y.a[1][1])%p;
	ans.a[1][0]=(x.a[1][0]*y.a[0][0]+x.a[1][1]*y.a[1][0])%p;
	ans.a[1][1]=(x.a[1][0]*y.a[0][1]+x.a[1][1]*y.a[1][1])%p;
   // cout<<ans.a[0][0]+ans.a[0][1]<<endl;//<<ans.a[0][1]<<endl;
 //   cout<<ans.a[1][0]<<"  "<<ans.a[1][1]<<endl; 
    //cout<<"&&&&&&&&&&"<<endl;
    return ans;
}

matrics power(matrics x , int n , int p)
{
    matrics ans , tmp;
    if(n==0) 
    {
        ans.a[0][0]=1; ans.a[0][1]=0;
        ans.a[1][0]=0; ans.a[1][1]=1;
        return ans;
    }
    if(n==1)  return x;
    tmp=power(x , n>>1 ,p);
    ans=mul(tmp , tmp ,p);
    if(n&1) ans=mul(ans, x , p);
    return ans;
}

int main()
{
    int a , b , p;
    __int64 sum ,ansa ,ansb;
    matrics  ans , tmp;
    while(scanf("%d%d%d",&a,&b,&p)!=EOF)
    {
        tmp.a[0][0]=1; tmp.a[0][1]=1;
        tmp.a[1][0]=1; tmp.a[1][1]=0;
        if(a==1 && b==2) { printf("%d\n",(1%p+1%p)%p); continue;}
        if(a==1) ansa=0;
        else if(a==2) ansa=1%p;
        else
        {
            ans = power( tmp , a-1 ,p);
            ansa=ans.a[0][0]+ans.a[0][1]-1;
        } //cout<<"*************************"<<endl;
     //   cout<<"ansa= "<<ansa+1<<endl;
       
        ans = power( tmp , b ,p);  //cout<<"ans.a[0][0]=  "<<ans.a[0][0]<<endl;
        ansb=ans.a[0][0]+ans.a[0][1]-1 ;
        //cout<<"ansb= "<<ansb<<endl;
        sum=((ansb-ansa)%p+p)%p;
        printf("%I64d\n",sum);
    }
    return 0;
}