标签：
题意：
求((b+√d)/2)^n的整数部分。b*b<d<10^18，n<10^18，d%4==1，b*b%4==1，模数约等于7*10^18
题解：
神题。由一些性质可以得出一个数列:An=bAn-1+(d-b*b)/4*An-2，且这个数列的通项公式为An=((b+√d)/2)^n+((b-√d)/2)^n，且由题目条件得(d-b*b)/4为正整数，故可以用矩阵乘法求出An，由于(b-√d)/2∈(-1,0]，故答案为(An)-1当且仅当n为偶数且b*b!=d。矩阵递推式：
An-1 An-2 b 1 An An-1
* =
0 0 (d-b*b)/4 0 0 0
反思：蒟蒻不知道矩乘不满足交换律，调了很久样例。同时由于模数太大，除了需要用unsigned long long外，乘法还要用快速乘（就是用快速幂的方法计算乘法）以防溢出。
代码：
1 #include <cstdio> 2 #include <cstring> 3 #include <algorithm> 4 #define ll unsigned long long 5 #define inc(i,j,k) for(int i=j;i<=k;i++) 6 #define mod 7528443412579576937 7 using namespace std; 8 9 struct M{ 10 ll a[5][5]; 11 M(){inc(i,0,1)inc(j,0,1)a[i][j]=0;} 12 }; 13 ll b,d,n; M st,ans; 14 ll cheng(ll a,ll b){ 15 if(b==0)return 0; if(b==1)return a; ll c=cheng(a,b>>1)%mod; 16 if(b&1)return ((c+c)%mod+a)%mod;else return (c+c)%mod; 17 } 18 M mul(M a,M b){ 19 M c; inc(i,0,1)inc(j,0,1)inc(k,0,1) 20 c.a[i][j]=(c.a[i][j]+cheng(a.a[i][k],b.a[k][j]))%mod; 21 return c; 22 } 23 M pow(M a,ll b){ 24 if(b==1)return a; M c=pow(a,b>>1); if(b&1)return mul(mul(c,c),a);else return mul(c,c); 25 } 26 int main(){ 27 scanf("%lld%lld%lld",&b,&d,&n); 28 if(n==0){printf("1"); return 0;} 29 st.a[0][0]=b%mod; st.a[0][1]=2; 30 if(n==1)ans=st;else{ 31 ans.a[0][0]=b%mod; ans.a[0][1]=1; ans.a[1][0]=(d-b*b)/4%mod; 32 ans=pow(ans,n-1); ans=mul(st,ans); 33 } 34 if(b*b!=d&&!(n&1))printf("%lld",(ans.a[0][0]+mod-1)%mod);else printf("%lld",ans.a[0][0]); 35 return 0; 36 }
20160710
标签：
原文地址：http://www.cnblogs.com/YuanZiming/p/5658579.html
