zoukankan      html  css  js  c++  java
  • 【BZOJ3328】PYXFIB(矩乘+单位根反演)

    点此看题面

    大致题意:(sum_{i=0}^{lfloorfrac nk floor}C_n^{i imes k} imes F_{i imes k})

    前言

    斐波那契?(nle 10^{18})?矩乘啊!

    然后就不会了。。。

    单位根反演

    考虑枚举(i imes k),其实可以看作是在([0,n])范围内枚举(i)并要求(k|i),即:

    [sum_{i=0}^n[k|i] imes C_n^i imes F_i ]

    看到([k|i]),我们就可以套上单位根反演的公式([k|n]=frac 1ksum_{i=0}^{k-1}omega_k^{ni})得到:

    [sum_{i=0}^n(frac 1ksum_{j=0}^{k-1}omega_k^{ij}) imes C_n^{i} imes F_i ]

    我们调整枚举顺序,改为先枚举(j),就得到:

    [frac 1ksum_{j=0}^{k-1}sum_{i=0}^nC_n^i imes F_i imes(omega_k^j)^i ]

    一个众所周知的事实,(F_i)可以表示为矩阵快速幂,即:(([1,1])表示右下角的格子)

    [F_i=egin{bmatrix}0 & 1\1 & 1end{bmatrix}^i[1,1] ]

    那么代入原式就有:

    [(frac 1ksum_{j=0}^{k-1}sum_{i=0}^nC_n^i imes (egin{bmatrix}0 & 1\1 & 1end{bmatrix} imesomega_k^j)^i)[1,1] ]

    观察这个式子,我们发现,它可以直接套用二项式定理得到:(其中(E)为单位矩阵)

    [frac 1ksum_{j=0}^{k-1}(omega_k^jegin{bmatrix}0&1\1&1end{bmatrix}+E)^i[1,1] ]

    由于(k)很小,直接枚举(k)然后做矩乘就可以了。

    代码

    #include<bits/stdc++.h>
    #define Tp template<typename Ty>
    #define Ts template<typename Ty,typename... Ar>
    #define Reg register
    #define RI Reg int
    #define Con const
    #define CI Con int&
    #define I inline
    #define W while
    #define K 20000
    #define LL long long
    using namespace std;
    LL n;int k,g,X;
    struct M
    {
    	int a[2][2];I M() {a[0][0]=a[0][1]=a[1][0]=a[1][1]=0;}
    	I int* operator [] (CI x) {return a[x];}
    	I M operator + (Con M& o) Con//矩阵加法
    	{
    		M t;t[0][0]=(a[0][0]+o.a[0][0])%X,t[0][1]=(a[0][1]+o.a[0][1])%X,
    		t[1][0]=(a[1][0]+o.a[1][0])%X,t[1][1]=(a[1][1]+o.a[1][1])%X;return t;
    	}
    	I M operator * (CI x) Con//矩阵乘数
    	{
    		M t;t[0][0]=1LL*a[0][0]*x%X,t[0][1]=1LL*a[0][1]*x%X,
    		t[1][0]=1LL*a[1][0]*x%X,t[1][1]=1LL*a[1][1]*x%X;return t;
    	}
    	I M operator * (Con M& o) Con//矩阵乘法
    	{
    		M t;t[0][0]=(1LL*a[0][0]*o.a[0][0]+1LL*a[0][1]*o.a[1][0])%X,
    		t[0][1]=(1LL*a[0][0]*o.a[0][1]+1LL*a[0][1]*o.a[1][1])%X,
    		t[1][0]=(1LL*a[1][0]*o.a[0][0]+1LL*a[1][1]*o.a[1][0])%X,
    		t[1][1]=(1LL*a[1][0]*o.a[0][1]+1LL*a[1][1]*o.a[1][1])%X;return t;
    	}
    	I M operator ^ (LL y) Con//矩阵快速幂
    	{
    		M x=*this,t;t[0][0]=t[1][1]=1;W(y) y&1&&(t=t*x,0),x=x*x,y>>=1;return t;
    	}
    }E,U,S;
    I int QP(RI x,RI y,CI X) {RI t=1;W(y) y&1&&(t=1LL*t*x%X),x=1LL*x*x%X,y>>=1;return t;}
    int p[K+5];I void GetGR(CI x)//求原根
    {
    	RI i,j,t=0,s=x-1;for(i=2;i*i<=s;++i) if(!(s%i)) {p[++t]=i;W(!(s%i)) s/=i;}//求出质因数
    	for(s^1&&(p[++t]=s),i=2;i<=x;++i)//枚举数
    	{
    		for(j=1;j<=t;++j) if(QP(i,(x-1)/p[j],x)==1) break;if(j>t) return (void)(g=i);//判断是否为原根
    	}
    }
    int main()
    {
    	E[0][0]=E[1][1]=U[0][1]=U[1][0]=U[1][1]=1;//初始化单位矩阵以及斐波那契矩阵
    	RI Tt,i,w,t;scanf("%d",&Tt);W(Tt--)
    	{
    		scanf("%lld%d%d",&n,&k,&X),S=M();
    		for(GetGR(X),w=QP(g,(X-1)/k,X),t=1,i=0;i^k;t=1LL*t*w%X,++i) S=S+(((U*t)+E)^n);//枚举k,累加矩阵,w为单位根
    		printf("%d
    ",1LL*QP(k,X-2,X)*S[1][1]%X);//输出答案,注意乘1/k
    	}return 0;
    }
    
  • 相关阅读:
    Linux防火墙--iptables学习
    LVS持久化
    LVS管理工具--ipvsadm
    Linux负载均衡--LVS(IPVS)
    一步步学习python
    驱动工程师需要的技能
    红外图像盲元补偿matlab实现源码与效果验证
    红外图像非均匀矫正——两点矫正
    夏日炎炎 python写个天气预报
    解决OV系列摄像头寄存器读数据无法收到的问题
  • 原文地址:https://www.cnblogs.com/chenxiaoran666/p/BZOJ3328.html
Copyright © 2011-2022 走看看