第五节 常系数齐次递推+BM算法
|
132
|
0
|
第四章 多项式与数值积分

1127 字
|
41 分钟

Contents

常系数齐次递推

题目大意：

求一个满足 k 阶齐次线性递推数列 a_i 的第 n 项，即：

a_n=\sum\limits_{i=1}^{k}f_i \times a_{n-i}

第一行两个数 n,k；第二行 k 个数，表示 f_1 \ f_2 \ \cdots \ f_k；第三行 k 个数，表示 a_0 \ a_1 \ \cdots \ a_{k-1}；输出一个数，表示 a_n \bmod 998244353。

#include <bits/stdc++.h>
using namespace std;
#define SZ(x) ((int)(x).size())
typedef vector<int> VI;
typedef long long ll;
const int P = 998244353;

void print(VI a) {
    int n = int(a.size());
    for (int i = 0; i < n - 1; i++) printf("%d ", a[i]);
    printf("%d\n", a[n - 1]);
}
inline void add(int &x, int y) {
    x += y;
    if (x >= P) x -= P;
}
inline void sub(int &x, int y) {
    x -= y;
    if (x < 0) x += P;
}
inline int mul(int x, int y) {
    return 1LL * x * y % P;
}
int ksm(int x, int y) {
    int res = 1;
    for (; y; y >>= 1, x = mul(x, x)) {
        if (y & 1) res = mul(res, x);
    }
    return res;
}
inline int inv(int a) {
    a %= P;
    if (a < 0) a += P;
    int b = P, u = 0, v = 1;
    while (a) {
        int t = b / a;
        b -= t * a; swap(a, b);
        u -= t * v; swap(u, v);
    }
    if (u < 0) u += P;
    return u;
}

namespace NTT {
    int bs = 1, rt = -1, mbs = -1;
    vector<int> rev = {0, 1}, rts = {0, 1};

    void init() {
        int temp = P - 1; mbs = 0;
        while (temp % 2 == 0) {
            temp >>= 1; ++mbs;
        }
        rt = 2;
        while (true) {
            if (ksm(rt, 1 << mbs) == 1 && ksm(rt, 1 << (mbs - 1)) != 1) {
                break;
            }
            ++rt;
        }
    }

    void ensure_base(int nbase) {
        if (mbs == -1) init();
        if (nbase <= bs) return;
        assert(nbase <= mbs);
        rev.resize(1 << nbase);
        for (int i = 0; i < 1 << nbase; ++i) {
            rev[i] = rev[i >> 1] >> 1 | (i & 1) << (nbase - 1);
        }
        rts.resize(1 << nbase);
        while (bs < nbase) {
            int z = ksm(rt, 1 << (mbs - 1 - bs));
            for (int i = 1 << (bs - 1); i < 1 << bs; ++i) {
                rts[i << 1] = rts[i];
                rts[i << 1 | 1] = mul(rts[i], z);
            }
            ++bs;
        }
    }

    void dft(VI &a) {
        int n = SZ(a), zeros = __builtin_ctz(n);
        ensure_base(zeros);
        int shift = bs - zeros;
        for (int i = 0; i < n; ++i) {
            if (i < rev[i] >> shift) {
                swap(a[i], a[rev[i] >> shift]);
            }
        }
        for (int i = 1; i < n; i <<= 1) {
            for (int j = 0; j < n; j += i << 1) {
                for (int k = 0; k < i; ++k) {
                    int x = a[j + k], y = mul(a[j + k + i], rts[i + k]);
                    a[j + k] = (x + y) % P;
                    a[j + k + i] = (x + P - y) % P;
                }
            }
        }
    }

    VI multiply(VI a, VI b) {
        int need = SZ(a) + SZ(b) - 1, nbase = 0;
        while (1 << nbase < need) ++nbase;
        ensure_base(nbase);
        int sz = 1 << nbase;
        a.resize(sz); b.resize(sz);
        bool equal = (a == b);
        dft(a);
        if (equal) b = a; else dft(b);
        int invsz = inv(sz);
        for (int i = 0; i < sz; ++i) {
            a[i] = mul(mul(a[i], b[i]), invsz);
        }
        reverse(a.begin() + 1, a.end());
        dft(a); a.resize(need);
        return a;
    }

    VI inverse(VI a) {
        int n = SZ(a), m = (n + 1) >> 1;
        if (n == 1) {
            return vector<int>(1, inv(a[0]));
        }
        else {
            vector<int> b = inverse(vector<int>(a.begin(), a.begin() + m));
            int need = n << 1, nbase = 0;
            while (1 << nbase < need) ++nbase;
            ensure_base(nbase);
            int sz = 1 << nbase;
            a.resize(sz); b.resize(sz);
            dft(a); dft(b);
            int invsz = inv(sz);
            for (int i = 0; i < sz; ++i) {
                a[i] = mul(mul(P + 2 - mul(a[i], b[i]), b[i]), invsz);
            }
            reverse(a.begin() + 1, a.end());
            dft(a); a.resize(n);
            return a;
        }
    }
}
using NTT::multiply;
using NTT::inverse;

VI& operator *= (VI &a, const VI &b) {
    if (min(a.size(), b.size()) < 128) {
        vector<int> c = a;
        a.assign(a.size() + b.size() - 1, 0);
        for (int i = 0; i < SZ(c); ++i) {
            for (int j = 0; j < SZ(b); ++j) {
                add(a[i + j], mul(c[i], b[j]));
            }
        }
    }
    else a = multiply(a, b);
    return a;
}
VI operator * (const VI &a, const VI &b) {
    VI c = a;
    return c *= b;
}
VI& operator /= (VI &a, const VI &b) {
    int n = SZ(a), m = SZ(b);
    if (n < m) a.clear();
    else {
        VI c = b;
        reverse(a.begin(), a.end());
        reverse(c.begin(), c.end());
        c.resize(n - m + 1);
        a *= inverse(c);
        a.erase(a.begin() + n - m + 1, a.end());
        reverse(a.begin(), a.end());
    }
    return a;
}
VI operator / (const VI &a, const VI &b) {
    vector<int> c = a;
    return c /= b;
}
VI& operator %= (VI &a, const VI &b) {
    int n = SZ(a), m = SZ(b);
    if (n >= m) {
        VI c = (a / b) * b;
        a.resize(m - 1);
        for (int i = 0; i < m - 1; ++i) {
            sub(a[i], c[i]);
        }
    }
    return a;
}
VI operator % (const VI &a, const VI &b) {
    vector<int> c = a;
    return c %= b;
}

VI ksm(VI a, int b, VI c) {
    VI res(1, 1);
    for(; b; b >>= 1, a = a * a % c) {
        if(b & 1) res = res * a % c;
    }
    return res;
}

// a 系数 b 初值 b[n+1]=a[0]*b[n]+a[1]*b[n-1]...
int solve(VI a, VI b, int n) {
    int k = (int)a.size();
    for (int i = 0; i < k; i++) a[i] = (P - a[i]) % P;
    for (int i = 0; i < k; i++) b[i] = (b[i] + P) % P;
    reverse(a.begin(), a.end());
    a.push_back(1);
    a = ksm(VI({0, 1}), n, a);
    a.resize(k);
    int ans = 0;
    for (int i = 0; i < k; i++) add(ans, mul(a[i], b[i]));
    return ans;
}

int n, k;
int main()
{
    cin >> n >> k;
    VI f(k), a(k);
    for(int i = 0; i < k; i++) scanf("%d", &f[i]);
    for(int i = 0; i < k; i++) scanf("%d", &a[i]);
    printf("%d\n", solve(f, a, n));
    return 0;
}

Berlekamp-Massey 算法

朴素 BM 算法

对质数取模；
设原数列长为 n，最短递推数列长为 k，要求第 m 项：
BM 算法复杂度为 O(n^2)，线性递推复杂度为 O(k^2logm)

#include <bits/stdc++.h>
typedef long long ll;
using namespace std;
#define SZ(x) ((int)(x).size())
#define pb push_back
typedef vector<int> VI;
const int P = 998244353;
int ksm(int a, int b) {
    int res = 1;
    for (; b; b >>= 1, a = 1LL * a * a % P) {
        if (b & 1) res = 1LL * res * a % P;
    }
    return res;
}

namespace linear_seq {
    const int N = 40010;
    ll res[N],base[N],_c[N],_d[N],_md[N];

    vector<int> Md;
    void mul(ll *a, ll *b, int k) {
        for (int i = 0; i < k + k; i++) _c[i]=0;
        for (int i = 0; i < k; i++) if (a[i]) {
            for (int j = 0; j < k; j++) {
                (_c[i + j] += a[i] * b[j]) %= P;
            }
        }
        for (int i = k + k - 1; i >= k; i--) if(_c[i]) {
            for (int j = 0; j < SZ(Md); j++) {
                (_c[i - k + Md[j]] -= _c[i] * _md[Md[j]]) %= P;
            }
        }
        for (int i = 0; i < k; i++) a[i]=_c[i];
    }

    // a 系数 b 初值 b[n+1]=a[0]*b[n]+...
    int solve(ll n, VI a, VI b) {
        ll ans = 0, pnt = 0;
        int k = SZ(a);
        assert(SZ(a) == SZ(b));
        for (int i = 0; i < k; i++) _md[k - 1 - i] = -a[i];
        _md[k] = 1; Md.clear();
        for (int i = 0; i < k; i++) {
            if (_md[i] != 0) Md.push_back(i);
            res[i] = base[i] = 0;
        }
        res[0] = 1;
        while ((1LL << pnt) <= n) pnt++;
        for (ll p = pnt; p >= 0; p--) {
            mul(res, res, k);
            if ((n >> p) & 1) {
                for (int i = k - 1; i >= 0; i--) {
                    res[i + 1] = res[i];
                }
                res[0] = 0;
                for (int j = 0; j < SZ(Md); j++) {
                    (res[Md[j]] -= res[k] * _md[Md[j]]) %= P;
                }
            }
        }
        for (int i = 0; i < k; i++) (ans += res[i] * b[i]) %= P;
        if (ans < 0) ans += P;
        return (int)ans;
    }

    VI BM(VI s) {
        VI C(1, 1), B(1, 1);
        int L = 0, m = 1, b = 1;
        for (int n = 0; n < SZ(s); n++) {
            ll d = 0;
            for (int i = 0; i <= L; i++) {
                (d += 1LL * C[i] * s[n - i]) %= P;
            }
            if (d == 0) m++;
            else if (2 * L <= n) {
                VI T = C;
                ll c = P - d * ksm(b, P-2) % P;
                while (SZ(C) < SZ(B) + m) C.push_back(0);
                for (int i = 0; i < SZ(B); i++) {
                    C[i + m] = (C[i + m] + c * B[i]) % P;
                }
                L = n + 1 - L; B = T; b = int(d); m = 1;
            }
            else {
                ll c = P - d * ksm(b, P-2) % P;
                while (SZ(C) < SZ(B) + m) C.push_back(0);
                for (int i = 0; i < SZ(B); i++) {
                    C[i + m] = (C[i + m] + c * B[i]) % P;
                }
                ++m;
            }
        }
        C.erase(C.begin());
        for (int &x : C) x = (P - x) % P;
        return C;
    }

    int gao(VI a, ll n) {
        VI c = BM(a);
        for (int i = 0, s = (int)c.size(); i < s; i++) {
            printf("%d%c", c[i], i == s - 1 ? '\n' : ' ');
        }
        return solve(n, c, VI(a.begin(), a.begin() + SZ(c)));
    }
}

int main() {
    int n, m; scanf("%d%d", &n, &m);
    VI v(n);
    for (int &x : v) scanf("%d", &x);
    printf("%d\n", linear_seq::gao(v, m));
    return 0;
}

多项式加速

与上边的常系数齐次线性递推相结合。

对质数取模；
设原数列长为 n，最短递推数列长为 k，要求第 m 项：
BM 算法复杂度为 O(n^2)，线性递推复杂度为 O(k\cdot logk\cdot logm)

#include <bits/stdc++.h>
using namespace std;
#define SZ(x) ((int)(x).size())
typedef vector<int> VI;
typedef long long ll;
const int P = 998244353;

void print(VI a) {
    int n = int(a.size());
    for (int i = 0; i < n - 1; i++) printf("%d ", a[i]);
    printf("%d\n", a[n - 1]);
}
inline void add(int &x, int y) {
    x += y;
    if (x >= P) x -= P;
}
inline void sub(int &x, int y) {
    x -= y;
    if (x < 0) x += P;
}
inline int mul(int x, int y) {
    return 1LL * x * y % P;
}
int ksm(int x, int y) {
    int res = 1;
    for (; y; y >>= 1, x = mul(x, x)) {
        if (y & 1) res = mul(res, x);
    }
    return res;
}
inline int inv(int a) {
    a %= P;
    if (a < 0) a += P;
    int b = P, u = 0, v = 1;
    while (a) {
        int t = b / a;
        b -= t * a; swap(a, b);
        u -= t * v; swap(u, v);
    }
    if (u < 0) u += P;
    return u;
}

namespace NTT {
    int bs = 1, rt = -1, mbs = -1;
    vector<int> rev = {0, 1}, rts = {0, 1};

    void init() {
        int temp = P - 1; mbs = 0;
        while (temp % 2 == 0) {
            temp >>= 1; ++mbs;
        }
        rt = 2;
        while (true) {
            if (ksm(rt, 1 << mbs) == 1 && ksm(rt, 1 << (mbs - 1)) != 1) {
                break;
            }
            ++rt;
        }
    }

    void ensure_base(int nbase) {
        if (mbs == -1) init();
        if (nbase <= bs) return;
        assert(nbase <= mbs);
        rev.resize(1 << nbase);
        for (int i = 0; i < 1 << nbase; ++i) {
            rev[i] = rev[i >> 1] >> 1 | (i & 1) << (nbase - 1);
        }
        rts.resize(1 << nbase);
        while (bs < nbase) {
            int z = ksm(rt, 1 << (mbs - 1 - bs));
            for (int i = 1 << (bs - 1); i < 1 << bs; ++i) {
                rts[i << 1] = rts[i];
                rts[i << 1 | 1] = mul(rts[i], z);
            }
            ++bs;
        }
    }

    void dft(VI &a) {
        int n = SZ(a), zeros = __builtin_ctz(n);
        ensure_base(zeros);
        int shift = bs - zeros;
        for (int i = 0; i < n; ++i) {
            if (i < rev[i] >> shift) {
                swap(a[i], a[rev[i] >> shift]);
            }
        }
        for (int i = 1; i < n; i <<= 1) {
            for (int j = 0; j < n; j += i << 1) {
                for (int k = 0; k < i; ++k) {
                    int x = a[j + k], y = mul(a[j + k + i], rts[i + k]);
                    a[j + k] = (x + y) % P;
                    a[j + k + i] = (x + P - y) % P;
                }
            }
        }
    }

    VI multiply(VI a, VI b) {
        int need = SZ(a) + SZ(b) - 1, nbase = 0;
        while (1 << nbase < need) ++nbase;
        ensure_base(nbase);
        int sz = 1 << nbase;
        a.resize(sz); b.resize(sz);
        bool equal = (a == b);
        dft(a);
        if (equal) b = a; else dft(b);
        int invsz = inv(sz);
        for (int i = 0; i < sz; ++i) {
            a[i] = mul(mul(a[i], b[i]), invsz);
        }
        reverse(a.begin() + 1, a.end());
        dft(a); a.resize(need);
        return a;
    }

    VI inverse(VI a) {
        int n = SZ(a), m = (n + 1) >> 1;
        if (n == 1) {
            return vector<int>(1, inv(a[0]));
        }
        else {
            vector<int> b = inverse(vector<int>(a.begin(), a.begin() + m));
            int need = n << 1, nbase = 0;
            while (1 << nbase < need) ++nbase;
            ensure_base(nbase);
            int sz = 1 << nbase;
            a.resize(sz); b.resize(sz);
            dft(a); dft(b);
            int invsz = inv(sz);
            for (int i = 0; i < sz; ++i) {
                a[i] = mul(mul(P + 2 - mul(a[i], b[i]), b[i]), invsz);
            }
            reverse(a.begin() + 1, a.end());
            dft(a); a.resize(n);
            return a;
        }
    }
}
using NTT::multiply;
using NTT::inverse;

VI& operator *= (VI &a, const VI &b) {
    if (min(a.size(), b.size()) < 128) {
        vector<int> c = a;
        a.assign(a.size() + b.size() - 1, 0);
        for (int i = 0; i < SZ(c); ++i) {
            for (int j = 0; j < SZ(b); ++j) {
                add(a[i + j], mul(c[i], b[j]));
            }
        }
    }
    else a = multiply(a, b);
    return a;
}
VI operator * (const VI &a, const VI &b) {
    VI c = a;
    return c *= b;
}
VI& operator /= (VI &a, const VI &b) {
    int n = SZ(a), m = SZ(b);
    if (n < m) a.clear();
    else {
        VI c = b;
        reverse(a.begin(), a.end());
        reverse(c.begin(), c.end());
        c.resize(n - m + 1);
        a *= inverse(c);
        a.erase(a.begin() + n - m + 1, a.end());
        reverse(a.begin(), a.end());
    }
    return a;
}
VI operator / (const VI &a, const VI &b) {
    vector<int> c = a;
    return c /= b;
}
VI& operator %= (VI &a, const VI &b) {
    int n = SZ(a), m = SZ(b);
    if (n >= m) {
        VI c = (a / b) * b;
        a.resize(m - 1);
        for (int i = 0; i < m - 1; ++i) {
            sub(a[i], c[i]);
        }
    }
    return a;
}
VI operator % (const VI &a, const VI &b) {
    vector<int> c = a;
    return c %= b;
}

VI ksm(VI a, ll b, VI c) {
    VI res(1, 1);
    for(; b; b >>= 1, a = a * a % c) {
        if(b & 1) res = res * a % c;
    }
    return res;
}

namespace linear_seq {
    const int N = 40010;
    ll res[N],base[N],_c[N],_d[N],_md[N];
    vector<int> Md;

    // a 系数 b 初值 b[n+1]=a[0]*b[n]+a[1]*b[n-1]...
    int solve(ll n, VI a, VI b) {
        int k = (int)a.size();
        for (int i = 0; i < k; i++) a[i] = (P - a[i]) % P;
        for (int i = 0; i < k; i++) b[i] = (b[i] + P) % P;
        reverse(a.begin(), a.end());
        a.push_back(1);
        a = ksm(VI({0, 1}), n, a);
        a.resize(k);
        int ans = 0;
        for (int i = 0; i < k; i++) add(ans, mul(a[i], b[i]));
        return ans;
    }
    VI BM(VI s) {
        VI C(1, 1), B(1, 1);
        int L = 0, m = 1, b = 1;
        for (int n = 0; n < SZ(s); n++) {
            ll d = 0;
            for (int i = 0; i <= L; i++) {
                (d += 1LL * C[i] * s[n - i]) %= P;
            }
            if (d == 0) m++;
            else if (2 * L <= n) {
                VI T = C;
                ll c = P - d * ksm(b, P-2) % P;
                while (SZ(C) < SZ(B) + m) C.push_back(0);
                for (int i = 0; i < SZ(B); i++) {
                    C[i + m] = (C[i + m] + c * B[i]) % P;
                }
                L = n + 1 - L; B = T; b = int(d); m = 1;
            }
            else {
                ll c = P - d * ksm(b, P-2) % P;
                while (SZ(C) < SZ(B) + m) C.push_back(0);
                for (int i = 0; i < SZ(B); i++) {
                    C[i + m] = (C[i + m] + c * B[i]) % P;
                }
                ++m;
            }
        }
        C.erase(C.begin());
        for (int &x : C) x = (P - x) % P;
        return C;
    }

    int gao(VI a, ll n) {
        VI c = BM(a);
        for (int i = 0, s = (int)c.size(); i < s; i++) {
            printf("%d%c", c[i], i == s - 1 ? '\n' : ' ');
        }
        return solve(n, c, VI(a.begin(), a.begin() + SZ(c)));
    }
}

int main(int argc, const char * argv[]) {
    int n, m; scanf("%d%d", &n, &m);
    VI v(n);
    for (int &x : v) scanf("%d", &x);
    printf("%d\n", linear_seq::gao(v, m));
    return 0;
}

广义 BM

可以对非质数取模。其余复杂度决定于常系数线性递推的复杂度。前面的代码和朴素 BM 一样，这里省略。

例题：2019牛客多校第九场 A

题目大意：

求

\sum_{i=0}^{n}F_i^m\ (mod\ 1000000000)

其中F_i是斐波那契数列的第i项。

#include <bits/stdc++.h>
#define pb push_back
#define SZ(x) ((int)(x).size())
using namespace std;
const int P = 1e9;
typedef long long ll;
typedef vector<ll> VI;
int ksm(int a, ll b) {
    int res = 1;
    for (; b; b >>= 1, a = 1LL * a * a % P) {
        if (b & 1) res = 1LL * res * a % P;
    }
    return res;
}
namespace linear_seq {
    const int N=10010;
    ll res[N], base[N], _c[N], _md[N];
    vector<int> Md;
    void mul(ll *a, ll *b, int k) {}
    // a 系数 b 初值 b[n+1]=a[0]*b[n]+...
    int solve(ll n, VI a, VI b) {}
    VI BM(VI s) {}

    static void extand(VI &a, size_t d, ll value = 0) {
        if (d <= a.size()) return;
        a.resize(d, value);
    }
    static void exgcd(ll a, ll b, ll &g, ll &x, ll &y) {
        if (!b) {x = 1; y = 0; g = a;}
        else {
            exgcd(b, a%b, g, y, x);
            y -= x * (a / b);
        }
    }
    static ll crt(const VI &c, const VI &m) {
        int n = SZ(c);
        ll M = 1, ans = 0;
        for (int i = 0; i < n; ++i) M *= m[i];
        for (int i = 0; i < n; ++i) {
            ll x, y, g, tm = M / m[i];
            exgcd(tm, m[i], g, x, y);
            ans = (ans + tm * x * c[i] % M) % M;
        }
        return (ans + M) % M;
    }
    static VI ReedsSloane(const VI &s, ll M) {
        auto inverse = [](ll a, ll m) {
            ll d, x, y;
            exgcd(a, m, d, x, y);
            return d == 1 ? (x % m + m) % m : -1;
        };
        auto L = [](const VI &a, const VI &b) {
            auto judge = [](const VI &x) -> bool {
                return x.size() > 1 || (x.size() == 1 && x[0]);
            };
            int da = judge(a) ? SZ(a) - 1 : -1000;
            int db = judge(b) ? SZ(b) : -1000;
            return max(da, db + 1);
        };
        auto prime_power = [&](const VI &s, ll M, ll p, ll e) {
            vector<VI> a(e), b(e), an(e), bn(e), ao(e), bo(e);
            VI t(e), u(e), r(e), to(e, 1), uo(e), pw(e + 1);;
            pw[0] = 1;
            for (int i = pw[0] = 1; i <= e; ++i) pw[i] = pw[i - 1] * p;
            for (ll i = 0; i < e; i++) {
                a[i] = {pw[i]}; an[i] = {pw[i]};
                b[i] = {0}; bn[i] = {s[0] * pw[i] % M};
                t[i] = s[0] * pw[i] % M;
                if (t[i] == 0) {t[i] = 1; u[i] = e;}
                else for(u[i] = 0; t[i] % p == 0; t[i] /= p, ++u[i]);
            }
            for (size_t k = 1; k < SZ(s); k++) {
                for (int g = 0; g < e; ++g) {
                    if (L(an[g], bn[g]) > L(a[g], b[g])) {
                        ao[g] = a[e - 1 - u[g]];
                        bo[g] = b[e - 1 - u[g]];
                        to[g] = t[e - 1 - u[g]];
                        uo[g] = u[e - 1 - u[g]];
                        r[g] = k - 1;
                    }
                }
                a = an; b = bn;
                for (int o = 0; o < e; ++o) {
                    ll d = 0;
                    for (size_t i = 0; i < a[o].size() && i <= k; i++) {
                        d = (d + a[o][i] * s[k - i]) % M;
                    }
                    if (d == 0) {t[o] = 1; u[o] = e; continue;}
                    for (u[o] = 0, t[o] = d; t[o] % p == 0;
                         t[o] /= p, ++u[o]);
                    int g = int(e - 1 - u[o]);
                    if (L(a[g], b[g]) == 0) {
                        extand(bn[o], k + 1);
                        bn[o][k] = (bn[o][k] + d) % M;
                        continue;
                    }
                    ll coef = t[o] * inverse(to[g], M) % M
                              * pw[u[o] - uo[g]] % M;
                    int m = int(k - r[g]);
                    extand(an[o], ao[g].size() + m);
                    extand(bn[o], bo[g].size() + m);
                    for (size_t i = 0; i < ao[g].size(); i++) {
                        an[o][i + m] -= coef * ao[g][i] % M;
                        if (an[o][i + m] < 0) an[o][i + m] += M;
                    }
                    while (an[o].size() && an[o].back() == 0) {
                        an[o].pop_back();
                    }
                    for (size_t i = 0; i < bo[g].size(); i++) {
                        bn[o][i + m] -= coef * bo[g][i] % M;
                        if (bn[o][i + m] < 0) bn[o][i + m] -= M;
                    }
                    while (bn[o].size() && bn[o].back() == 0) {
                        bn[o].pop_back();
                    }
                }
            }
            return make_pair(an[0], bn[0]);
        };
        vector<tuple<ll, ll, int>> fac;
        for (ll i = 2; i * i <= M; ++i) {
            if (M % i == 0) {
                ll cnt = 0, pw = 1;
                while (M % i == 0) {M /= i; ++cnt; pw *= i;}
                fac.emplace_back(pw, i, cnt);
            }
        }
        if (M > 1) fac.emplace_back(M, M, 1);
        vector<VI> as;
        size_t n = 0;
        for (auto &&x: fac) {
            ll M, p, e;
            VI a, b;
            tie(M, p, e) = x;
            auto ss = s;
            for (auto &&x: ss) x %= M;
            tie(a, b) = prime_power(ss, M, p, e);
            as.emplace_back(a);
            n = max(n, a.size());
        }
        VI a(n), c(as.size()), m(as.size());
        for (size_t i = 0; i < n; ++i) {
            for (size_t j = 0; j < as.size(); ++j) {
                m[j] = std::get<0>(fac[j]);
                c[j] = i < as[j].size() ? as[j][i] : 0;
            }
            a[i] = crt(c, m);
        }
        a.erase(a.begin());
        for (ll &x : a) x = (P - x) % P;
        return a;
    }
    ll gao(VI a, ll n, ll mod, bool prime = true) {
        VI c;
        if (prime) c = BM(a);
        else c = ReedsSloane(a, mod);
        return solve(n, c, VI(a.begin(), a.begin()+SZ(c)));
    }
}
int main()
{
    ll n, m; cin >> n >> m;
    vector<ll> v; v.pb(0); v.pb(1);
    for(int i = 2; i < 1000; ++i) v.pb((v[i-1] + v[i-2]) % P);
    for(int i = 1; i < 1000; ++i) v[i] = ksm((int)v[i], m);
    for(int i = 1; i < 1000; ++i) v[i] = (v[i - 1] + v[i]) % P;
    printf("%lld\n", linear_seq::gao(v,n,P,false));
    return 0;
}
暂无评论

发送评论 编辑评论 


				

			

						

				

					

						

							

								
							

							
						

					

				

				

					

						

							

								
							

							
						

					

				

				

					

						

							

								
							

							
							
													

					

				

			

			
			

				

					
				

			

				

											

							
							
						

																				
					
											
						

	

		
|´・ω・)ノ
ヾ(≧∇≦*)ゝ
(☆ω☆)
（╯‵□′）╯︵┴─┴
￣﹃￣
(/ω＼)
∠( ᐛ 」∠)＿
(๑•̀ㅁ•́ฅ)
→_→
୧(๑•̀⌄•́๑)૭
٩(ˊᗜˋ*)و
(ノ°ο°)ノ
(´இ皿இ｀)
⌇●﹏●⌇
(ฅ´ω`ฅ)
(╯°A°)╯︵○○○
φ(￣∇￣o)
ヾ(´･ ･｀｡)ノ"
( ง ᵒ̌皿ᵒ̌)ง⁼³₌₃
(ó﹏ò｡)
Σ(っ °Д °;)っ
( ,,´･ω･)ﾉ"(´っω･｀｡)
╮(╯▽╰)╭
o(*////▽////*)q
＞﹏＜
( ๑´•ω•) "(ㆆᴗㆆ)
😂
😀
😅
😊
🙂
🙃
😌
😍
😘
😜
😝
😏
😒
🙄
😳
😡
😔
😫
😱
😭
💩
👻
🙌
🖕
👍
👫
👬
👭
🌚
🌝
🙈
💊
😶
🙏
🍦
🍉
😣
Source: github.com/k4yt3x/flowerhd
	

	

		
颜文字
Emoji
小恐龙
花!
	


									

			

			
			
		

	






上一篇
下一篇