首页 > 代码库 > [NOI2011] 兔农 矩阵乘法,矩阵的等比数列求和
[NOI2011] 兔农 矩阵乘法,矩阵的等比数列求和
#include<iostream>#include<cstdio>#include<algorithm>#include<cstring>using namespace std;typedef long long ll;ll n,p;ll k;#define N 4000000ll a[N],fr[N];struct sq{ ll a[3][3]; sq(){memset(a,0,sizeof(a));} sq operator*(sq x)const{ sq ans; for(ll i=1;i<=2;i++)for(ll j=1;j<=2;j++)for(ll z=1;z<=2;z++)(ans.a[i][j]+=a[i][z]*x.a[z][j])%=p; return ans; } void out(){ for(ll i=1;i<=2;i++){ for(ll j=1;j<=2;j++){ cout<<a[i][j]<<" "; } cout<<endl; } cout<<endl; }}tt;ll hash[N];struct sq2{ ll a[5][5]; sq2(){memset(a,0,sizeof(a));} sq2 operator*(sq2 x){ sq2 ans; for(ll i=1;i<=4;i++)for(ll j=1;j<=4;j++)for(ll z=1;z<=4;z++)(ans.a[i][j]+=a[i][z]*x.a[z][j])%=p; return ans; } void out(){ for(ll i=1;i<=4;i++){ for(ll j=1;j<=4;j++){ cout<<a[i][j]<<" "; } cout<<endl; } cout<<endl; }};sq ff(sq x,ll y){ sq ans,k=x; for(ll i=1;i<=2;i++)ans.a[i][i]=1; for(ll i=0;(1LL<<i)<=y;i++){ if((1LL<<i)&y)ans=ans*k; k=k*k; } return ans;}sq cheng(sq x,ll y){ sq2 tr,ans1; for(ll i=1;i<=2;i++)for(ll j=1;j<=2;j++)ans1.a[i][j]=x.a[i][j]; for(ll i=1;i<=2;i++){ for(ll j=1;j<=2;j++){ tr.a[i][j]=x.a[i][j]; } tr.a[i][i+2]=1; tr.a[i+2][i+2]=1; ans1.a[i][i+2]=1; } for(ll i=0;(1LL<<i)<=y;i++){ if(y&(1LL<<i)){ ans1=ans1*tr; } tr=tr*tr; } sq kk; for(ll i=1;i<=2;i++)for(ll j=1;j<=2;j++){ kk.a[i][j]=ans1.a[i][j+2]; } return kk;}ll gcd(ll &p,ll &q,ll x,ll y){ if(x==0){ p=1,q=1; return y; } ll r=gcd(q,p,y%x,x); p-=q*(y/x); return r;}ll gcd(ll x,ll y){ return x==0?y:gcd(y%x,x);}ll bn;ll bb[N],cir;ll ben[N];ll res;ll ni(ll x){ ll p,q; if(gcd(x,k)!=1)return 0; gcd(p,q,x,k); return (p%k+k)%k;}int main(){ tt.a[1][1]=1; tt.a[2][1]=1; tt.a[1][2]=1; sq u; u.a[1][1]=1; a[1]=1;a[2]=1; scanf("%lld%lld%lld",&n,&k,&p); for(ll i=3;i<=k*3;i++){ a[i]=(a[i-1]+a[i-2])%k; if(!fr[a[i]])fr[a[i]]=i; } ll i=0; ll j=1; cir=n+1; ll ss=0; for(;;){ ll u=ni(j); if(!u||!fr[u])break; ll g=a[fr[u]-1]*j%k; if(i+fr[u]>n)break; if(!ben[g]){ ben[g]=1; bb[++bn]=i+fr[u]; hash[g]=bn; }else{ bb[++bn]=i+fr[u]; ss=hash[g]; cir=bb[bn]-bb[ss]; break; } i=i+fr[u]; j=g; } sq yy=ff(tt,cir); res=(u*ff(tt,n)).a[1][2]; for(ll i=1;i<=bn;i++){ if(i<=ss){ (res-=(u*ff(tt,n-bb[i])).a[1][1])%=p; }else{ll uu=(n-bb[i])%cir; (res-=((u*(ff(tt,uu)*cheng(yy,(n-bb[i])/cir)))).a[1][1])%=p; }} res=(res%p+p)%p; cout<<res;}
声明:以上内容来自用户投稿及互联网公开渠道收集整理发布,本网站不拥有所有权,未作人工编辑处理,也不承担相关法律责任,若内容有误或涉及侵权可进行投诉: 投诉/举报 工作人员会在5个工作日内联系你,一经查实,本站将立刻删除涉嫌侵权内容。