转自:http://www.klogk.com/posts/hdu4565/
题意非常简单, a,b,n 都是正整数,求
Sn=⌈(a+b√)n⌉%m,(a−1)2<b<a2
这个题目也是2008年Google Codejam Round 1A的C题。
做法其实非常简单,记 (a+b√)n 为 An ,配项
Cn=An+Bn=(a+b√)n+(a−b√)n
两项恰好共轭,所以 Cn 是整数。又根据限制条件
(a−1)2<b<a2⇒0<a−b√<1⇒0<(a−b√)n<1⇒Bn<1
也就是说 Cn=⌈An⌉
Sn=(Cn)%m
求 Cn 的方法是递推。 对 Cn 乘以 (a+b√)+(a−b√)
Cn[(a+b√)+(a−b√)]=[(a+b√)n+(a−b√)n][(a+b√)+(a−b√)]=(a+b√)n+1+(a−b√)n+1+(a+b√)n(a−b√)+(a−b√)n(a+b√)=Cn+1+(a2−b)(a+b√)n−1+(a2−b)(a−b√)n−1=Cn+1+(a2−b)Cn−1
于是
Cn+1=2aCn−(a2−b)Cn−1
把这个递推式写成矩阵形式
[Cn+1Cn]=[2a1−(a2−b)0][CnCn−1]
于是就可以用矩阵快速幂来做了
[Cn+1Cn]=[2a1−(a2−b)0]n[C1C0]
一个需要注意的地方是快速幂的过程中由于取余可能出现负数,要化为正数。
#include<cstdio>
using namespace std;
__int64 A[3][3],s[3][3],tmp[3][3];
void fun(__int64 n,__int64 m)
{
__int64 i,j,k;
for(i=1;i<=2;i++)
for(j=1;j<=2;j++)
{
s[i][j]=A[i][j];
}
--n;
while(n)
{
// printf("n=%d\n",n);
if(n&1)
{
for(i=1;i<=2;i++)
for(j=1;j<=2;j++)
{
tmp[i][j]=A[i][j];
A[i][j]=0;
}
for(i=1;i<=2;i++)
for(j=1;j<=2;j++)
for(k=1;k<=2;k++)
{
A[i][j]=((A[i][j]+tmp[i][k]*s[k][j])%m+m)%m;
}
}
for(i=1;i<=2;i++)
for(j=1;j<=2;j++)
{
tmp[i][j]=s[i][j];
s[i][j]=0;
}
for(i=1;i<=2;i++)
for(j=1;j<=2;j++)
for(k=1;k<=2;k++)
{
s[i][j]=((s[i][j]+tmp[i][k]*tmp[k][j])%m+m)%m;
}
n>>=1;
}
}
int main()
{
__int64 a,b,n,m;
__int64 i,j;
while(scanf("%I64d%I64d%I64d%I64d",&a,&b,&n,&m)==4)
{
if(n==1)
{
printf("%I64d\n",2*a%m);
continue;
}
A[1][1]=2*a; A[1][2]=b-a*a; A[2][1]=1; A[2][2]=0;
fun(n-1,m);
/* for(i=1;i<=2;i++)
{
for(j=1;j<=2;j++)
printf("%I64d ",A[i][j]);
printf("\n\n");
}
*/
/* for(i=1;i<=2;i++)
{
for(j=1;j<=2;j++)
printf("%I64d ",s[i][j]);
printf("\n\n");
}
*/
printf("%I64d\n",((A[1][1]*2*a+A[1][2]*2)%m+m)%m);
}
return 0;
}