CodeForces - 1106F Lunar New Year and a Recursive Sequence

题意：

已知 $f_1=f_2=\cdots =f_{k-1}=1,f_i=\prod _{j=1}^k{f}_{i-j}^{b_j}(j>k),f_n\equiv m(\mathrm{m}\mathrm{o}\mathrm{d}p)$，给定 $k,b_i(1\le i\le k),n,m$，求满足的一个 $f_k(1\le f_k<p)$。$(p=998244353,k\le 100,k<n\le 10^9,1\le m<p)$

链接：

https://codeforces.com/contest/1106/problem/F

解题思路：

$3$ 为 $998244353$ 的原根，那么 $f(n)$ 可以被唯一表示成 $3^{g(n)mod(p-1)}modp$，取离散对数后，则可化为线性递推式。
$$\begin{gathered} 3^{g(n)mod(p-1)}\equiv \prod _{i=1}^k{3}^{b_{i}g(n-i)mod(p-1)}(\mathrm{m}\mathrm{o}\mathrm{d}p) \\ g(n)\equiv \sum _{i=1}^k{b}_{i}g(n-i)(\mathrm{m}\mathrm{o}\mathrm{d}p-1) \\ 矩阵快速幂可得g(n)\equiv cg(k)(\mathrm{m}\mathrm{o}\mathrm{d}p-1) \\ exgcd解出g(k),则f(k)\equiv 3^{g(k)}(\mathrm{m}\mathrm{o}\mathrm{d}p) \end{gathered}$$

参考代码：

#include<bits/stdc++.h>
 
using namespace std;
typedef long long ll;
typedef pair<int, int> pii;
#define sz(a) ((int)a.size())
#define pb push_back
#define lson (rt << 1)
#define rson (rt << 1 | 1)
#define gmid (l + r >> 1)
const int maxn = 1e2 + 5;
const int inf = 0x3f3f3f3f;
const int mod = 998244353;

struct Mat{
    
    const static int N = 1e2 + 5;
    const static int mod = 998244353 - 1;
    ll a[N][N]; int n;
    Mat(int n) : n(n) {memset(a, 0, sizeof a);}
    void setI(){
        memset(a, 0, sizeof a);
        for(int i = 0; i < n; ++i) a[i][i] = 1;
    }
    Mat operator * (const Mat &o) const{
        Mat ret(n);
        for(int i = 0; i < n; ++i)
            for(int j = 0; j < n; ++j)
                for(int k = 0; k < n; ++k)
                    ret.a[i][j] = (ret.a[i][j] + a[i][k] * o.a[k][j]) % mod;
        return ret;
    }
};

Mat qpow(Mat a, ll b){

    Mat ret(a.n); ret.setI();
    while(b){

        if(b & 1) ret = ret * a;
        a = a * a;
        b >>= 1;
    }
    return ret;
}

ll qpow(ll a, ll b, ll p){

    ll ret = 1;
    while(b){

        if(b & 1) ret = ret * a % p;
        a = a * a % p;
        b >>= 1;
    }
    return ret;
}

void exgcd(ll a, ll b, ll &x, ll &y, ll &d){

    if(!b) x = 1, y = 0, d = a;
    else exgcd(b, a % b, y, x, d), y -= a / b * x;
}

ll exgcd(ll a, ll b, ll c){

    ll x, y, d;
    exgcd(a, b, x, y, d);
    if(c % d) return -1;
    x *= c / d;
    ll tmp = b / d;
    x = (x % tmp + tmp) % tmp;
    return x;
}

ll inv(ll a, ll p){

    return exgcd(a, p, 1);
}

map<ll, ll> mp;
ll BSGS(ll a, ll b, ll p){

    a %= p, b %= p;
    if(!a && b) return -1;
    mp.clear();
    ll lim = ceil(sqrt(p));
    ll cur = 1;
    for(int i = 0; i < lim; ++i){

        if(!mp.count(cur)) mp[cur] = i;
        cur = cur * a % p;
    } 
    ll inva = inv(cur, p);
    for(int i = 0; i <= lim; ++i){

        if(mp.count(b)) return i * lim + mp[b];
        b = b * inva % p;
    }
    return -1;
}

ll b[maxn];
ll n, k, m;

int main(){
 
    ios::sync_with_stdio(0); cin.tie(0);
    cin >> k;
    for(int i = 1; i <= k; ++i) cin >> b[i];
    cin >> n >> m;
    Mat A(k);
    for(int i = 0; i < k; ++i) A.a[0][i] = b[i + 1];
    for(int i = 1; i < k; ++i) A.a[i][i - 1] = 1;
    A = qpow(A, n - k);
    ll g_n = BSGS(3, m, mod);
    ll g_k = exgcd(A.a[0][0], mod - 1, g_n);
    if(g_k == -1){

        cout << "-1" << endl;
        return 0;
    }
    ll f_k = qpow(3, g_k, mod);
    cout << f_k << endl;
    return 0;
}