1. 程式人生 > >Codeforces 536F Lunar New Year and a Recursive Sequence | BSGS/exgcd/矩陣乘法

Codeforces 536F Lunar New Year and a Recursive Sequence | BSGS/exgcd/矩陣乘法

離散對數 log new 考題 block memset ons ORC tps

我詐屍啦!

高三退役選手好不容易拋棄天利和金考卷打場CF,結果打得和shi一樣……還因為queue太長而unrated了!一個學期不敲代碼實在是忘幹凈了……

沒分該沒分,考題還是要訂正的 =v= 歡迎閱讀本題解!

P.S. 這幾個算法我是一個也想不起來了 TAT

題目鏈接

Codeforces 536F Lunar New Year and a Recursive Sequence 新年和遞推數列

題意描述

某數列\(\{f_i\}\)遞推公式:\[f_i = (\prod_{j=1}^kf_{i-j}^{b_j}) \bmod p\]

其中\(b\)是已知的長度為\(k\)的數列,\(p = 998244353\)

\(f_1 = f_2 = ... = f_{k-1} = 1\)\(f_k\)未知。

給出兩個數\(n, m\),構造一個\(f_k\)使得\(f_n = m\),無解輸出-1。

\(k \le 100, n \le 10^9\)

題解

數論!真令人頭禿!

首先這個數據範圍讓人想到什麽?矩陣乘法!

矩陣乘法想推這個全是乘法和乘方的遞推數列咋辦?取對數!離散對數!

於是這道題關鍵的兩個考點就被你發現啦!

(然而我太菜了,並不能發現 = =)

什麽是離散對數?

眾所周知(眾==學過NTT的人等),這個喜聞樂見的模數\(p = 998244353\)有個原根\(g=3\)\(g^i(0\le i < P - 1)\)

\(1\le x < P\)一一對應。那麽類比我們學過的對數,稱這個\(i\)\(x\)的離散對數。

令數列\(h_i\)\(f_i\)的離散對數。

那麽有遞推式:\[h_i = (\sum_{j=1}^kb_j\cdot h_{i-j}) \bmod (p - 1)\]

其中\(h_1 = h_2 = ... = h_{k-1} = 0\)。註意模數變成了\(p - 1\)(費馬小定理)。

這個就可以用矩陣加速了!如果我們把\(h_k\)設為1帶進去,求得\(f_n = c\),那麽有\(h_n = c \cdot h_k \bmod (p - 1)\)

\(h_n\)即為\(m\)

的離散對數,用BSGS可求;

exgcd解剛才這個同余方程即可得到\(h_k\)

\(f_k = g^{h_k}\),快速冪即可得到\(f_k\)

如果exgcd發現沒有解的話就輸出-1。

是不是思路非常清晰啊~

代碼

#include <cstdio>
#include <cstring>
#include <cmath>
#include <algorithm>
#include <iostream>
#include <cassert>
#define space putchar(' ')
#define enter putchar('\n')
using namespace std;
typedef long long ll;
template <class T>
void read(T &x){
    char c;
    bool op = 0;
    while(c = getchar(), c < '0' || c > '9')
        if(c == '-') op = 1;
    x = c - '0';
    while(c = getchar(), c >= '0' && c <= '9')
        x = x * 10 + c - '0';
    if(op) x = -x;
}
template <class T>
void write(T x){
    if(x < 0) putchar('-'), x = -x;
    if(x >= 10) write(x / 10);
    putchar('0' + x % 10);
}

const int N = 102, P = 998244353, P2 = 998244352, G = 3;
int K;
ll b[N], n, m, C;

namespace BSGS {
    const int S = 32000, M = 2000000;
    int cnt = 0, adj[M + 5], nxt[S + 5];
    ll key[S + 5], val[S + 5];
    void insert(ll K, ll V){
        int p = K % M;
        key[++cnt] = K;
        val[cnt] = V;
        nxt[cnt] = adj[p];
        adj[p] = cnt;
    }
    ll search(ll K){
        for(int u = adj[K % M]; u; u = nxt[u])
            if(key[u] == K) return val[u];
        return -1;
    }
    void init(){
        ll sum = 1;
        for(int i = 1; i <= S; i++)
            sum = sum * G % P;
        ll tot = 1;
        for(int i = 1; (i - 1) * S < P - 1; i++)
            tot = tot * sum % P, insert(tot, i * S);
    }
    ll log(ll x){
        ll sum = 1, ret;
        for(int i = 1; i <= S; i++){
            sum = sum * G % P;
            ret = search(sum * x % P);
            if(~ret && ret < P - 1) return ret - i;
        }
        assert(0);
        return -1;
    }
}

struct matrix {
    ll g[N][N];
    matrix(){
        memset(g, 0, sizeof(g));
    }
    matrix(int x){
        memset(g, 0, sizeof(g));
        for(int i = 1; i <= K; i++)
            g[i][i] = 1;
    }
    matrix operator * (const matrix &b){
        matrix c;
        for(int i = 1; i <= K; i++)
            for(int j = 1; j <= K; j++)
                for(int k = 1; k <= K; k++)
                    c.g[i][j] = (c.g[i][j] + g[i][k] * b.g[k][j]) % P2;
        return c;
    }
};

ll qpow(ll a, ll x){
    ll ret = 1;
    while(x){
        if(x & 1) ret = ret * a % P;
        a = a * a % P;
        x >>= 1;
    }
    return ret;
}
matrix qpow(matrix a, ll x){
    matrix ret(1);
    while(x){
        if(x & 1) ret = ret * a;
        a = a * a;
        x >>= 1;
    }
    return ret;
}
ll calcC(){
    matrix ret, op;
    ret.g[K][1] = 1;
    for(int i = 1; i < K; i++)
        op.g[i][i + 1] = 1;
    for(int i = 1; i <= K; i++)
        op.g[K][i] = b[K - i + 1];
    ret = qpow(op, n - K) * ret;
    return ret.g[K][1];
}
void exgcd(ll a, ll b, ll &g, ll &x, ll &y){
    if(!b) return (void)(x = 1, y = 0, g = a);
    exgcd(b, a % b, g, y, x);
    y -= x * (a / b);
}
ll solve(ll A, ll B){ //Ax % P2 == B, solve x
    ll a = A, b = P2, g, x, y;
    exgcd(a, b, g, x, y);
    if(B % g) return -1;
    x *= B / g, y *= B / g;
    ll t = b / g;
    x = (x % t + t) % t;
    return x;
}

int main(){

    BSGS::init();
    read(K);
    for(int i = 1; i <= K; i++) read(b[i]);
    read(n), read(m);
    C = calcC();
    m = BSGS::log(m);
    ll ans = solve(C, m);
    if(ans == -1) puts("-1");
    else write(qpow(G, ans)), enter;

    return 0;
}

Codeforces 536F Lunar New Year and a Recursive Sequence | BSGS/exgcd/矩陣乘法