A sequence S n is defined as:
Where a, b, n, m are positive integers.┌x┐is the ceil of x. For example, ┌3.14┐=4. You are to calculate S n.
You, a top coder, say: So easy!
Input
There are several test cases, each test case in one line contains four positive integers: a, b, n, m. Where 0< a, m < 2 15, (a-1) 2< b < a 2, 0 < b, n < 2 31.The input will finish with the end of file.
Output
For each the case, output an integer S n.
Sample Input
2 3 1 2013 2 3 2 2013 2 2 1 2013
Sample Output
4 14 4
题解:
转载于:https://blog.csdn.net/bigbigship/article/details/45842365
求 f(x) = ceil( (a +sqrt(b))^n )
我们设An = (a +sqrt(b))^n , Bn =(a - sqrt(b))^n;
Cn = An +Bn;
因为An, Bn共轭,所以Cn是一个整数
根据题意, (a-1)^2 < b < a^2 ==> a-1 < sqrt(b) < a;
因此Cn = ceil( An )
Cn * [(a + sqrt(b)) +(a - sqrt(b))]
==> (a +sqrt(b))^(n+1) +(a +sqrt(b))^(n+1) + (a +sqrt(b))*(a -sqrt(b))^(n)+(a -sqrt(b))*(a +sqrt(b))^(n)
==> Cn+1 + (a*a - b)Cn-1
==> Cn+1 = 2*aCn + (b-a*a)*Cn-1;
公式推出来了剩下的就直接用矩阵加速就搞定了
代码 :
#include<stdio.h>
#include<string.h>
#include<iostream>
#include<cmath>
#include<algorithm>
using namespace std;
typedef long long ll;
const int mod=1000000007;
ll mmod;
struct mat
{
int r,c;
ll m[5][5];
void clear()
{
for(int i=1; i<=r; i++)memset(m[i],0,sizeof(m[i]));
}
};
int read()
{
int x=0;
char c=getchar();
while(c>'9'||c<'0')c=getchar();
while(c>='0'&&c<='9')
{
x=x*10+c-'0';
c=getchar();
}
return x;
}
mat MatMul(mat &m1,mat &m2)
{
mat tmp;
tmp.r=m1.r;
tmp.c=m2.c;
int i,j,k;
for(i=1; i<=tmp.r; i++)
{
for(j=1; j<=tmp.c; j++)
{
ll t=0;
for(k=1; k<=m1.c; k++)
{
t=(t+(m1.m[i][k]*m2.m[k][j])%mmod)%mmod;
}
tmp.m[i][j]=t;
}
}
return tmp;
}
mat MatQP(mat &a,int n)
{
mat ans,tmp=a;
ans.r=ans.c=a.r;
memset(ans.m,0,sizeof(ans.m));
for(int i=1; i<=ans.r; i++)
{
ans.m[i][i]=1;
}
while(n)
{
if(n&1)ans=MatMul(ans,tmp);
n>>=1;
tmp=MatMul(tmp,tmp);
}
return ans;
}
ll QP(ll a,ll n) //快速幂。
{
ll tmp=a,ans=1;
while(n)
{
if(n&1)ans=ans*tmp%mod;
tmp=tmp*tmp%mod;
n>>=1;
}
return ans%mod;
}
int main()
{
ll a,b,n;
while(~scanf("%lld%lld%lld%lld",&a,&b,&n,&mmod))
{
mat t,tmp;
t.r=2;
t.c=2;
t.clear();
t.m[1][1]=2*a%mmod;
t.m[1][2]=1;
t.m[2][1]=((b%mmod-a*a%mmod)+mmod)%mmod;
t.m[2][2]=0;
mat q;
q.r=1;
q.c=2;
ll c2=ceil((a+sqrt(b))*(a+sqrt(b)));
ll c1=ceil((a+sqrt(b)));
q.m[1][1]=c2%mmod;
q.m[1][2]=c1%mmod;
if(n==1)
{
cout<<2*a%mmod<<endl;
continue;
}
else
{
tmp=MatQP(t,n-2);
tmp=MatMul(q,tmp);
ll dd=tmp.m[1][1]%mmod;
cout<<dd<<endl;
}
}
}