2025年hdu4565（矩阵快速幂）

科技前沿 • 2025-02-18 16:16 • 阅读 54

大家好，我是讯享网，很高兴认识大家。

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

求 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=; 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; } } }

2025年hdu4565（矩阵快速幂）

相关推荐