Codeforces 1106F Lunar New Year and a Recursive Sequence (数学线性代数线性递推数论BSGS扩展欧几里得算法)

Posted 2021-02-07 suncongbo

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

数论板子X合一？

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

题意: 给一个(n)阶递推序列(f_k=prod^{n}_{i=1} f_{k-i}b_imod 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;
}