HDU - 5322 Hope

题目大意

求$$f_i=(i-1)!\sum_{j=0}^{i-1}(i-j)^2{f_j\over j!}$$

简要题解

为求$f_i$我们需要知道$f_0,\cdots,f_{i-1}$，考虑cdq分治，把卷积拆开成关于已知的$f_i$和还没计算出来的部分，发现已知部分还是卷积形式，求出来累加上去就好。

此外还写了个求原根，暴力枚举原根$a$，然后判断是否只有$a^{p-1}=1(\mod p)$而$a^{p-1\over p_i}=1(\mod p)$皆不成立，其中$p_i$是$p$的所有素因子。

#include <bits/stdc++.h>
using namespace std;
namespace my_header {
#define pb push_back
#define mp make_pair
#define pii pair<int, int>
#define vec vector<int>
#define pc putchar
#define clr(t) memset(t, 0, sizeof t)
#define pse(t, v) memset(t, v, sizeof t)
#define bl puts("")
#define wn(x) wr(x), bl
#define ws(x) wr(x), pc(' ')
    const int INF = 0x3f3f3f3f;
    typedef long long LL;
    typedef double DB;
    inline char gchar() {
        char ret = getchar();
        for(; (ret == '\n' || ret == '\r' || ret == ' ') && ret != EOF; ret = getchar());
        return ret; }
    template<class T> inline void fr(T &ret, char c = ' ', int flg = 1) {
        for(c = getchar(); (c < '0' || '9' < c) && c != '-'; c = getchar());
        if (c == '-') { flg = -1; c = getchar(); }
        for(ret = 0; '0' <= c && c <= '9'; c = getchar())
            ret = ret * 10 + c - '0';
        ret = ret * flg; }
    inline int fr() { int t; fr(t); return t; }
    template<class T> inline void fr(T&a, T&b) { fr(a), fr(b); }
    template<class T> inline void fr(T&a, T&b, T&c) { fr(a), fr(b), fr(c); }
    template<class T> inline char wr(T a, int b = 10, bool p = 1) {
        return a < 0 ? pc('-'), wr(-a, b, 0) : (a == 0 ? (p ? pc('0') : p) : 
            (wr(a/b, b, 0), pc('0' + a % b)));
    }
    template<class T> inline void wt(T a) { wn(a); }
    template<class T> inline void wt(T a, T b) { ws(a), wn(b); }
    template<class T> inline void wt(T a, T b, T c) { ws(a), ws(b), wn(c); }
    template<class T> inline void wt(T a, T b, T c, T d) { ws(a), ws(b), ws(c), wn(d); }
    template<class T> inline T gcd(T a, T b) {
        return b == 0 ? a : gcd(b, a % b); }
    template<class T> inline T fpw(T b, T i, T _m, T r = 1) {
        for(; i; i >>= 1, b = b * b % _m)
            if(i & 1) r = r * b % _m;
        return r; }
};
using namespace my_header;

const int MAXN = 2e6 + 100;
namespace NTT {

const int MOD = 998244353;

int inv[MAXN];


vec getPrime(vec&isPrime, int rng = 1000000) {
    isPrime.resize(rng);
    fill(isPrime.begin(), isPrime.end(), 1);
    vec pri;
    isPrime[0] = isPrime[1] = 0;
    inv[1] = 1;
    for (int i = 2; i < rng; ++i) {
        if (isPrime[i]) {
            pri.pb(i);
            inv[i] = fpw((LL)i, MOD-2LL, (LL)MOD);
        }
        for (int j = 0; j < (int)pri.size() && pri[j] * 1LL * i < rng; ++j) {
            isPrime[i * 1LL * pri[j]] = 0;
            inv[i * 1LL * pri[j]] = 1LL * inv[i] * inv[pri[j]] % MOD;
        }
    }
    return pri;
}

void dfs_enum(vector<pii>&fac, int p, int d, vec&res) {
    if (d == (int)fac.size()) {
        res.pb(p);
        return ;
    }
    int v = fac[d].first, t = fac[d].second, c = 1;
    dfs_enum(fac, p, d + 1, res);
    for (int i = 0; i < t; ++i) {
        c = c * v;
        dfs_enum(fac, p * c, d + 1, res);
    }
}

vec enumerate(vector<pii>&fac) {
    vec res;
    dfs_enum(fac, 1, 0, res);
    return res;
}

vec factorize(LL val) {
    vec isPrime;
    vector<pii> fac;
    vec pri = getPrime(isPrime);
    for (int i = 0; i < (int)pri.size(); ++i) {
        int p = pri[i];
        if (val % p == 0) {
            fac.pb(mp(p, 0));
            while (val % p == 0) {
                val /= p;
                ++fac.back().second;
            }
        }
    }
    if (val != 1)
        fac.pb(mp(val, 1));
    return enumerate(fac);
}

int getPrimeRoot(LL val) {
    int pr;
    vec f = factorize(val - 1);
    for (pr = 2; ; ++pr) {
        if (fpw((LL)pr, val - 1LL, (LL)MOD) == 1) {
            bool ok = true;
            for (int i = 0; i < (int)f.size(); ++i)
                if (f[i] != (val - 1) && fpw((LL)pr, (LL)f[i], (LL)MOD) == 1)
                    ok = false;
            if (ok) return pr;
        }
    }
}

int GN[24], rGN[24];
int getGN(int m, int r) {
    if (r == 1) {
        if (GN[m] != 0) return GN[m];
        return GN[m] = fpw(3LL, (MOD-1LL) / (1<<m), (LL)MOD);
    } else {
        if (rGN[m] != 0) return rGN[m];
        return rGN[m] = fpw((LL)GN[m], MOD-2LL, (LL)MOD);
    }
}


int b[MAXN];

void dft(int *a, int m, int r) {
    if (m == 0) return ;
    int N = 1<<m;
    memcpy(b, a, (sizeof a[0]) * N);
    for (int i = 0; i < (N>>1); ++i)
        a[i] = b[i<<1], a[(N>>1) + i] = b[(i<<1)|1];
    dft(a, m - 1, r), dft(a + (N>>1), m - 1, r);
    int gn = getGN(m, r), c = 1;
    for (int i = 0; i < (N>>1); c = 1LL * c * gn % MOD, ++i) {
        b[i] = (MOD + (a[i] + 1LL * c * a[i + (N>>1)] % MOD) % MOD) % MOD;
        b[i + (N>>1)] = (MOD + (a[i] - 1LL * c * a[i + (N>>1)] % MOD) % MOD) % MOD;
    }
    memcpy(a, b, (sizeof a[0]) * N);
}

};
using namespace NTT;

int j2[MAXN], fac[MAXN], ifac[MAXN], ta[MAXN], tb[MAXN], f[MAXN];

void calc(int l, int r) {
    if (l == r) return;
    int m = (l + r) >> 1;
    calc(l, m);
    int len = r - l + 1, t = 0;
    for (; (1<<t) < len; ++t);
    for (int i = l; i <= m; ++i)
        ta[i - l] = f[i] * 1LL * ifac[i] % MOD;
    for (int i = m + 1; i < m + (1<<t); ++i)
        ta[i - l] = 0;
    for (int i = 0; i < (1<<t); ++i)
        tb[i] = j2[i] % MOD;
    dft(ta, t, 1);
    dft(tb, t, 1);
    for (int i = 0; i < (1<<t); ++i)
        ta[i] = ta[i] * 1LL * tb[i] % MOD;
    dft(ta, t, -1);
    for (int i = m + 1; i <= r; ++i)
        (f[i] += 1LL * fac[i - 1] * ta[i - l] % MOD * inv[1<<t] % MOD) %= MOD;
    calc(m + 1, r);
}

int main() {
#ifdef lol
    freopen("G.in", "r", stdin);
    freopen("G.out", "w", stdout);
#endif

    getPrimeRoot(MOD);

    fac[0] = 1;
    j2[0] = 0;
    for (int i = 1; i < MAXN; ++i) {
        fac[i] = fac[i-1] * 1LL * i % MOD;
        j2[i] = i * 1LL * i % MOD;
    }
    ifac[MAXN-1] = fpw((LL)fac[MAXN-1], MOD - 2LL, (LL)MOD);
    for (int i = MAXN-2; 0 <= i; --i)
        ifac[i] = (i + 1LL) * ifac[i + 1] % MOD;

    
    f[0] = 1;
    calc(0, 100000);
    int x;
    while (scanf("%d", &x) != EOF)
        wt(f[x]);

    return 0;
}