zoukankan      html  css  js  c++  java
  • Codeforces 1106F Lunar New Year and a Recursive Sequence (线性代数、线性递推、数论、BSGS、扩展欧几里得算法)

    哎呀大水题。。我写了一个多小时。。好没救啊。。

    数论板子X合一?

    注意: 本文中变量名称区分大小写。

    题意: 给一个(n)阶递推序列(f_k=prod^{n}_{i=1} f_{k-i}^{b_i}mod P)其中(P=998244353), 输入(b_1,b_2,...,b_n)以及已知(f_1,f_2,...,f_{n-1}=1), 再给定一个数(m)和第(m)项的值(f_m), 求出一个合法的(f_n)值使得按照这个值递推出来的序列满足第(m)项的值为给定的(f_m).

    题解: 首先一个显然的结论是(f_m)可以表示成(prod^{n}_{i=1} f_i^{a_i}), 而且由于(i=1,2,...,n-1)(f_i)的任何次幂都为(1), 因此就是(f_m=f_n^{a}). 令(A(m))(f_m)(f_n)的次数,则有(A[1..n]=[0,0,0,0,...,0,1]), (A_m=sum^{n-1}_{i=1} A(m-i)b_i (m>n)), 即(A)数组满足一个常系数线性递推序列。因此可以用矩阵乘法在(O(n^3log m))的时间内求出(A(m)). 注意因为是指数的运算(((a^n)^m=a^{nm})), 根据费马小定理,这个指数应该模(phi(P)=P-1)而不是(P) (((a^n)^mmod P=a^{nmmod (P-1)}mod P))

    求出来(a=A(m))之后这题就变成了,(f_m=f_n^amod P), 已知(f_m, a), 求出一组合法的(f_n).

    根据常识,(998244353)有原根(3), 我们下文令(G=3) (实际上任何一个原根均可). 设(f_m=G^p, f_n=G^q), 则有(G^pequiv (G^q)^a (mod P)), (pequiv qa(mod P-1)), 然后用BSGS求离散对数(p), exgcd解出(q)就可以了啊……

    时间复杂度(O(sqrt Plog P+n^3log P))

    坑: 注意解同余方程的时候那个(P)的系数不要设成负的。

    代码

    #include<cstdio>
    #include<cstdlib>
    #include<cstring>
    #include<map>
    #define llong long long
    using namespace std;
    
    const int N = 100;
    const int G = 3;
    const int P = 998244353;
    llong quickpow(llong x,llong y)
    {
    	llong cur = x,ret = 1ll;
    	for(int i=0; y; i++)
    	{
    		if(y&(1ll<<i)) {y-=(1ll<<i); ret = ret*cur%P;}
    		cur = cur*cur%P;
    	}
    	return ret;
    }
    struct Matrix
    {
    	llong a[N+3][N+3]; int sz1,sz2,sz;
    	void init() {for(int i=1; i<=sz1; i++) for(int j=1; j<=sz2; j++) a[i][j] = 0ll;}
    	Matrix() {}
    	Matrix(int _sz) {sz = sz1 = sz2 = _sz; init();}
    	Matrix(int _sz1,int _sz2) {sz1 = _sz1,sz2 = _sz2; init();}
    	void uinit(int _sz) {sz = sz1 = sz2 = _sz; for(int i=1; i<=sz; i++) for(int j=1; j<=sz; j++) a[i][j] = (i==j)?1:0;}
    	void output() {for(int i=1; i<=sz1; i++) {for(int j=1; j<=sz2; j++) printf("%lld ",a[i][j]); puts("");}}
    };
    Matrix operator *(Matrix x,Matrix y)
    {
    	Matrix ret = Matrix(x.sz1,y.sz2);
    	for(int i=1; i<=x.sz1; i++)
    	{
    		for(int j=1; j<=x.sz2; j++)
    		{
    			for(int k=1; k<=y.sz2; k++)
    			{
    				ret.a[i][j] = (ret.a[i][j]+x.a[i][k]*y.a[k][j])%(P-1ll);
    			}
    		}
    	}
    	return ret;
    }
    Matrix mquickpow(Matrix x,llong y)
    {
    	Matrix cur = x,ret; ret.uinit(x.sz);
    	for(int i=0; y; i++)
    	{
    		if(y&(1ll<<i)) {y-=(1ll<<i); ret = ret*cur;}
    		cur = cur*cur;
    	}
    	return ret;
    }
    namespace BSGS
    {
    	const int B = 31595;
    	map<llong,int> mp;
    	void init()
    	{
    		llong bs = quickpow(G,B); llong j = 1ll;
    		for(int i=0; i<=P; i+=B,j=(j*bs)%P)
    		{
    			mp[j] = i;
    		}
    	}
    	llong Logarithm(llong x)
    	{
    		llong j = 1ll;
    		for(int i=0; i<=B; i++,j=(j*G)%P)
    		{
    			llong tmp = x*j%P;
    			if(mp.count(tmp))
    			{
    				llong ret = (mp[tmp]-i+(P-1))%(P-1);
    				return ret;
    			}
    		}
    		return P-1;
    	}
    }
    Matrix mA,mB,mC;
    llong a[N+3],b[N+3];
    int n; llong m,p,q,lq,lx;
    
    llong gcd(llong x,llong y) {return y==0 ? x : gcd(y,x%y);}
    void exgcd(llong _a,llong _b,llong &_x,llong &_y)
    {
    	if(_b==0ll) {_x = 1ll,_y = 0ll; return;}
    	exgcd(_b,_a%_b,_x,_y);
    	llong tmp = _x; _x = _y; _y = tmp-(_a/_b)*_y;
    }
    llong CongruenceEquation(llong _a,llong _b,llong mod)
    {
    	llong g = gcd(_a,mod),x,y;
    	exgcd(_a/g,mod/g,x,y);
    	return (x*(_b/g)%mod+mod)%mod;
    }
    
    int main()
    {
    	BSGS::init();
    	scanf("%d",&n);
    	for(int i=1; i<=n; i++) scanf("%I64d",&b[i]);
    	scanf("%I64d",&m);
    	mA = Matrix(1,n); mA.a[1][1] = 1ll;
    	mB = Matrix(n); for(int i=1; i<n; i++) mB.a[i][i+1] = 1ll; for(int i=1; i<=n; i++) mB.a[i][1] = b[i];
    	mC = mA*mquickpow(mB,m-n); p = mC.a[1][1]; scanf("%I64d",&q);
    	lq = BSGS::Logarithm(q);
    	if(lq%gcd(P-1,p)!=0) {printf("-1
    "); return 0;}
    	lx = CongruenceEquation(p,lq,P-1);
    	llong ans = quickpow(G,lx);
    	printf("%I64d
    ",ans);
    	return 0;
    }
    
  • 相关阅读:
    DGA域名可以是色情网站域名
    使用cloudflare加速你的网站隐藏你的网站IP
    167. Two Sum II
    leetcode 563. Binary Tree Tilt
    python 多线程
    leetcode 404. Sum of Left Leaves
    leetcode 100. Same Tree
    leetcode 383. Ransom Note
    leetcode 122. Best Time to Buy and Sell Stock II
    天津Uber优步司机奖励政策(12月28日到12月29日)
  • 原文地址:https://www.cnblogs.com/suncongbo/p/10352633.html
Copyright © 2011-2022 走看看