标签: 多项式

  • CFGYM103415A CCPC2021广州A Math Ball

    题意: 有$n$种球,每种有无限个,同时第$i$种球有一个代价$c_i$,你要拿不超过$w$个球。如果最后第$i$种球你拿了$k_i$个,那么你会获得$\prod_{1\leq i\leq n}{k_i^{c_i}}$的权值,求所有合法方案的权值和。$n\leq 1e5,\sum{c_i}\leq 1e5,w\leq 10^{18}$

    $$  \text{考虑对于价值是}c_i\text{的球,构造生成函数}  \\  F_{c_i}\left( x \right) =\sum_{n\ge 0}{n^{c_i}x^n}  \\  \text{这样}\frac{\prod_i{F_{c_i}\left( x \right)}}{1-x}\text{的}w\text{次项即为答案}  \\  \text{设}F_k\left( x \right) =\sum_{n\ge 0}{n^kx^n},\text{显然可得}F_k\left( x \right) =x\frac{\mathrm{d}F_{k-1}\left( x \right)}{\mathrm{d}x}  \\  \text{进一步递推可得}, F_k\left( x \right) =\frac{\sum_{0\le i\le k}{T\left( k,i \right) x^i}}{\left( 1-x \right) ^{k+1}},\text{其中}T\left( k,i \right) \text{表示欧拉数}  \\  \text{考虑如何快速计算欧拉数}  \\  \text{首先由具体数学可得}  \\  \sum_{0\le i\le k}{T\left( k,i \right) \times \left( z+1 \right) ^i}=\sum_{0\le i\le k}{S_2\left( k,i \right) \times i!\times z^{k-i}}  \\  \text{进一步推导得欧拉数的通项公式}T\left( n,k \right) =\sum_{0\le i\le k}{\begin{array}{c}    \binom{n+1}{i}\left( -1 \right) ^i\left( k+1-i \right) ^n\\  \end{array}}  \\  \text{构造卷积形式}T\left( n,k \right) =\left( n+1 \right) !\sum_{0\le i\le k}{\frac{\left( -1 \right) ^i}{i!\left( n+1-i \right) !}}\times \left( k+1-i \right) ^n  \\  \text{卷积即可求出一行欧拉数。}  \\  \text{现在我们可以对每一个}c_i\text{计算出}F_{c_i}\left( x \right) \text{的分子,并且得到他们分子的乘积,设为}F\left( x \right)   \\  \text{设他们分母的乘积,再乘个}1-x\text{为}\left( 1-x \right) ^k  \\  \text{则答案即为}\frac{F\left( x \right)}{\left( 1-x \right) ^k}\text{的}w\text{次项系数}  \\  \text{即}\sum_{0\le i\le n}{\left[ x^i \right] F\left( x \right) \binom{w-i+k-1}{w-i}}\left( \text{广义二项式定理展开分母} \right)   \\  \text{令}a_i=\left[ x^i \right] F\left( x \right)   \\  \text{即}\sum_{0\le i\le n}{a_i\binom{w-i+k-1}{k-1}}  \\  \text{后者可以通过递推快速转移。}  \\    \\    \\  \frac{x^a}{\left( 1-x \right) ^b}=x^a\left( \sum_{n\ge 0}{\binom{-b}{n}\left( -1 \right) ^n}x^n \right)   \\  =x^a\left( \sum_{n\ge 0}{\frac{\left( -b \right) ^{\underline{n}}}{n!}\left( -1 \right) ^n}x^n \right)   \\  =x^a\left( \sum_{n\ge 0}{\frac{\left( -b-0 \right) \left( -b-1 \right) \left( -b-2 \right) ..\left( -b-\left( n-1 \right) \right)}{n!}\left( -1 \right) ^n}x^n \right)   \\  =x^a\left( \sum_{n\ge 0}{\frac{\left( b+0 \right) \left( b+1 \right) \left( b+2 \right) ..\left( b+\left( n-1 \right) \right)}{n!}}x^n \right)   \\  =x^a\left( \sum_{n\ge 0}{\frac{\left( b+n-1 \right) ^{\underline{n}}}{n!}}x^n \right)   \\  =\sum_{n\ge 0}{\frac{\left( b+n-1 \right) ^{\underline{n}}}{n!}x^{n+a}}  \\  =\sum_{n\ge 0}{\binom{b+n-1}{n}x^{n+a}}  \\  =\sum_{n\ge a}{\binom{b+n-a-1}{n-a}x^n}  \\  =\sum_{n\ge a}{\binom{b+n-a-1}{b-1}x^n}  \\  =\sum_{n\ge a}{\begin{array}{c}     \frac{\left( b-a+n-1 \right) ^{\underline{b-1}}}{\left( b-1 \right) !}x^n\\  \end{array}}  \\  \text{令}t=b-1  \\  \text{则}  \\  \sum_{n\ge a}{\begin{array}{c}      \frac{\left( t-a+n \right) ^{\underline{t}}}{t!}x^n\\  \end{array}}  \\  \text{从}a\text{推进到}a+1  \\  \frac{\left( t-a+n \right) ^{\underline{t}}}{\left( t-a+n-1 \right) ^{\underline{t}}}=\frac{t-a+n}{n-a}  \\    \\  10*9*8/\left( 9*8*7 \right) =10/7  \\    \\  \left( \frac{x^a}{\left( 1-x \right) ^b}\text{的}n\text{次项} \right)   \\    \\  \binom{n+b-1}{n}=\binom{n+b-1}{b-1}  \\  \binom{n+b-1}{n}\rightarrow \binom{n+b}{n+1}  \\    \\  \frac{\binom{n+b}{n+1}}{\binom{n+b-1}{n}}=\frac{n+b}{\left( n+1 \right)}  \\    \\  \frac{\left( n-2+b-i \right) !\left( b-1 \right) !\left( n-i \right) !}{\left( b-1 \right) !\left( n-i-1 \right) !\left( n-i+b-1 \right) !}  \\    \\  \frac{\left( n-i+b-2 \right) !\left( n-i \right)}{\left( n-i+b-1 \right) !}  \\  \frac{\left( n-i \right)}{\left( n-i+b-1 \right)}  \\    \\    $$
    #include <cinttypes>
    #include <cstring>
    #include <fstream>
    #include <iostream>
    #include <random>
    using int_t = int;
    using std::cin;
    using std::cout;
    using std::endl;
    
    const int_t LARGE = 2e6;
    
    using i64 = int64_t;
    #ifdef NTTCNT
    std::ofstream nttcnt("nttcnt.txt");
    #endif
    const int_t mod = 998244353;
    const int_t g = 3;
    
    int_t power(int_t b, int_t i) {
        int_t r = 1;
        if (i < 0)
            i = ((i64)i % (mod - 1) + mod - 1) % (mod - 1);
        while (i) {
            if (i & 1)
                r = (i64)r * b % mod;
            b = (i64)b * b % mod;
            i >>= 1;
        }
        return r;
    }
    
    void makeflip(int_t* arr, int_t size2) {
        int_t len = (1 << size2);
        arr[0] = 0;
        for (int_t i = 1; i < len; i++) {
            arr[i] = (arr[i >> 1] >> 1) | ((i & 1) << (size2 - 1));
        }
    }
    
    int_t upper2n(int_t x) {
        int_t r = 0;
        while ((1 << r) < x)
            r++;
        return r;
    }
    template <int_t arg = 1>
    void transform(int_t* A, int_t size2, int_t* flip) {
        // #define int_t i64
        int_t len = (1 << size2);
    #ifdef NTTCNT
        nttcnt << len << endl;
    #endif
        for (int_t i = 0; i < len; i++) {
            int_t r = flip[i];
            if (r > i)
                std::swap(A[i], A[r]);
        }
        for (int_t i = 2; i <= len; i *= 2) {
            int_t mr = power(g, (i64)arg * (mod - 1) / i);
            for (int_t j = 0; j < len; j += i) {
                int_t curr = 1;
                for (int_t k = 0; k < i / 2; k++) {
                    int_t u = A[j + k], v = (i64)curr * A[j + k + i / 2] % mod;
                    A[j + k] = ((i64)u + v) % mod;
                    A[j + k + i / 2] = ((i64)u - v + mod) % mod;
                    curr = (i64)curr * mr % mod;
                }
            }
        }
        // #undef int_t
    }
    void poly_mul(const int_t* A, int_t n, const int_t* B, int_t m, int_t* C) {
        /*
           计算n次多项式A与m次多项式B的乘积
        */
        int_t size2 = upper2n(n + m + 1);
        int_t len = 1 << size2;
        static int_t T1[LARGE], T2[LARGE];
        for (int_t i = 0; i < len; i++) {
            if (i <= n)
                T1[i] = A[i];
            else
                T1[i] = 0;
            if (i <= m)
                T2[i] = B[i];
            else
                T2[i] = 0;
        }
        static int_t fliparr[LARGE];
        makeflip(fliparr, size2);
        transform(T1, size2, fliparr);
        transform(T2, size2, fliparr);
        for (int_t i = 0; i < len; i++)
            T1[i] = (i64)T1[i] * T2[i] % mod;
        transform<-1>(T1, size2, fliparr);
        int_t inv = power(len, -1);
        for (int_t i = 0; i <= n + m; i++)
            C[i] = (i64)T1[i] * inv % mod;
    }
    
    const int_t FLARGE = 2e5 + 10;
    
    int_t fact[FLARGE + 1], inv[FLARGE + 1], factInv[FLARGE + 1];
    void euler_num(int_t n, int_t* out) {
        static int_t P1[LARGE];
        static int_t P2[LARGE];
        for (int_t i = 0; i <= n; i++) {
            P1[i] = (i64)((i % 2) ? (mod - 1) : 1) * factInv[i] % mod *
                    factInv[n + 1 - i] % mod;
            P2[i] = power(i, n);
        }
        poly_mul(P1, n, P2, n, out);
        for (int_t i = 0; i <= n; i++)
            out[i] = (i64)out[i] * fact[n + 1] % mod;
    }
    
    using poly_t = std::vector<int_t>;
    poly_t poly_dcmul(const poly_t* P, int_t left, int_t right) {
        if (left == right)
            return P[left];
        int_t mid = (left + right) / 2;
        auto lret = poly_dcmul(P, left, mid);
        auto rret = poly_dcmul(P, mid + 1, right);
        poly_t ret;
        ret.resize(lret.size() - 1 + rret.size() - 1 + 1);
        poly_mul(&lret[0], lret.size() - 1, &rret[0], rret.size() - 1, &ret[0]);
        while (!ret.empty() && ret.back() == 0)
            ret.pop_back();
        return ret;
    }
    int_t C(int_t n, int_t m) {
        return (i64)fact[n] * factInv[m] % mod * factInv[n - m] % mod;
    }
    i64 getC(i64 n, int_t m) {
        if (n < 0 || n < m)
            return 0;
        i64 prod = 1;
        for (int_t i = 0; i < m; i++)
            prod = prod * (n % mod - i + mod) % mod;
        return prod * factInv[m] % mod;
    }
    int main() {
        std::ios::sync_with_stdio(false);
        cin.tie(0), cout.tie(0);
        {
            fact[0] = fact[1] = inv[1] = factInv[0] = factInv[1] = 1;
            for (int_t i = 2; i <= FLARGE; i++) {
                fact[i] = (i64)fact[i - 1] * i % mod;
                inv[i] = (i64)(mod - mod / i) * inv[mod % i] % mod;
                factInv[i] = (i64)factInv[i - 1] * inv[i] % mod;
            }
        }
        int_t n;
        i64 w;
        cin >> n >> w;
    
        static poly_t up[int_t(1e5 + 10)];
        static poly_t cache[int_t(1e5 + 10)];
        int_t sum = 0;
        for (int_t i = 1; i <= n; i++) {
            int_t c;
            cin >> c;
            sum += c + 1;
            poly_t& curr = up[i];
            if (!cache[c].empty()) {
                curr = cache[c];
                continue;
            }
            curr.resize(2 * c + 3);
            euler_num(c, &curr[0]);
            curr.resize(c + 1);
            cache[c] = curr;
    #ifdef DEBUG
            cout << "up " << i << " = ";
            for (auto t : curr)
                cout << t << " ";
            cout << endl;
    #endif
        }
        sum++;
        // for (int_t i = 0; i <= sum; i++) {
        //     down[i] = (i64)((i % 2) ? (mod - 1) : 1) * C(sum, i) % mod;
        // }
        auto upprod = poly_dcmul(up, 1, n);
        // cout << upprod.size() << endl;
        upprod.resize(sum + 1);
    #ifdef DEBUG
        cout << "up prod = ";
        for (auto t : upprod)
            cout << t << " ";
        cout << endl;
        cout << "sum = " << sum << endl;
    #endif
        // cout << "down size = " << sum << endl;
        // w -= cutcount;
        // if (w < 0) {
        //     cout << 0 << endl;
        //     return 0;
        // }
        // int_t ret = poly_divat(&upprod[0], down, upprod.size() - 1, sum, w);
        i64 ret = 0;
        {
            i64 b = sum;
            i64 n = w;
            i64 curr = getC(n - 0 + b - 1, b - 1);
    // C(n-i+b-1,b-1)到C(n-(i+1)+b-1,b-1)
    // getC(n - i + b - 1, b - 1)
    #ifdef DEBUG
            cout << "init curr = " << curr << endl;
    #endif
            for (i64 i = 0; i < upprod.size(); i++) {
                ret = (ret + upprod[i] * curr % mod) % mod;
                curr = curr * (n % mod - i + mod) % mod *
                       power((n % mod - i + b - 1 + (i64)3 * mod) % mod, -1) % mod;
                if (curr == 0 || n <= i - n + 1)
                    break;
    #ifdef DEBUG
                cout << "curr = " << curr << " at i = " << i << " with n = " << n
                     << endl;
    #endif
            }
            // int_t curr = getC(w - b + b - 1, b - 1);
            // for (i64 i = std::max<i64>(0, w - sum); i < upprod.size(); i++) {
            //     ret = (ret + curr * upprod[i] % mod) % mod;
            //     curr = curr * (n % mod + n % mod - sum % mod + mod + i) % mod *
            //            power(n + 1, -1) % mod;
            // }
            // i64 n = w;
            // i64 t = sum - 1;
            // i64 prod = 1;
            // // i64 curr = (t + n) % mod;
            // for (int_t i = 0; i < t; i++)
            //     prod = (prod * (t + n % mod - i) % mod);
            // for (int_t a = 0; a < upprod.size(); a++) {
            //     if (n >= a) {
            //         ret = (ret + prod * factInv[t] % mod * upprod[a] % mod) %
            //         mod;
            //     }
            //     prod = prod * (t - a + n % mod) % mod *
            //            power((n % mod - a + mod) % mod, -1) % mod;
            // }
        }
        cout << ret << endl;
        return 0;
    }

     

  • 跑的比较快的多项式板子

    #include <cstring>
    #include <iostream>
    #include <random>
    using int_t = long long int;
    using std::cin;
    using std::cout;
    using std::endl;
    
    const int_t LARGE = 2e6;
    
    const int_t mod = 998244353;
    const int_t g = 3;
    int_t power(int_t b, int_t i) {
        int_t r = 1;
        if (i < 0)
            i = (i % (mod - 1) + mod - 1) % (mod - 1);
        while (i) {
            if (i & 1)
                r = r * b % mod;
            b = b * b % mod;
            i >>= 1;
        }
        return r;
    }
    
    void makeflip(int_t* arr, int_t size2) {
        int_t len = (1 << size2);
        arr[0] = 0;
        for (int_t i = 1; i < len; i++) {
            arr[i] = (arr[i >> 1] >> 1) | ((i & 1) << (size2 - 1));
        }
    }
    
    int_t upper2n(int_t x) {
        int_t r = 0;
        while ((1 << r) < x)
            r++;
        return r;
    }
    template <int_t arg = 1>
    void transform(int_t* A, int_t size2, int_t* flip) {
        int_t len = (1 << size2);
        for (int_t i = 0; i < len; i++) {
            int_t r = flip[i];
            if (r > i)
                std::swap(A[i], A[r]);
        }
        for (int_t i = 2; i <= len; i *= 2) {
            int_t mr = power(g, arg * (mod - 1) / i);
            for (int_t j = 0; j < len; j += i) {
                int_t curr = 1;
                for (int_t k = 0; k < i / 2; k++) {
                    int_t u = A[j + k], v = curr * A[j + k + i / 2] % mod;
                    A[j + k] = (u + v) % mod;
                    A[j + k + i / 2] = (u - v + mod) % mod;
                    curr = curr * mr % mod;
                }
            }
        }
    }
    void poly_inv(const int_t* A, int_t n, int_t* result) {
        /*
        这里的n是模x^n
        计算B(x)*A(x) = 1 mod x^n, 其中A(x)已知
        假设已知A(x)*B(x) = 1 mod x^{ceil(n/2)}
        假设C(x)*A(x) = 1 mod x^n
        (A(x)B(x)-1)^2 = A^2(x)B^2(x)-2A(x)B(x)+1= 0
        A(x)B^2(x)-2B(x)+C(x) = 0
        C(x) = 2B(x)-A(x)B^2(x)
        */
        // #ifdef DEBUG
        //     cout << "at " << n << endl;
        // #endif
        if (n == 1) {
            result[0] = power(A[0], -1);
            return;
        }
        int_t next = n / 2 + n % 2;
        poly_inv(A, next, result);
        //次数不要选错了,应该用n次的A和B去卷
        int_t size2 = upper2n(n + 2 * next + 1);
        static int_t X[LARGE];
        static int_t Y[LARGE];
        int_t len = (1 << size2);
        // 写错设置范围了!
        memset(X + n, 0, sizeof(X[0]) * (len - n));
        memset(Y + next, 0, sizeof(Y[0]) * (len - next));
        memcpy(X, A, sizeof(A[0]) * n);
        memcpy(Y, result, sizeof(result[0]) * next);
        static int_t fliparr[LARGE];
        makeflip(fliparr, size2);
        transform<1>(X, size2, fliparr);
        transform<1>(Y, size2, fliparr);
    
        for (int_t i = 0; i < len; i++) {
            X[i] = (2 * Y[i] - X[i] * Y[i] % mod * Y[i] % mod + mod) % mod;
        }
        transform<-1>(X, size2, fliparr);
        const int_t inv = power(len, -1);
        for (int_t i = 0; i < n; i++)
            result[i] = X[i] * inv % mod;
    #ifdef DEBUG
        cout << "poly inv at " << n << endl;
        cout << "result = ";
        for (int_t i = 0; i < next; i++)
            cout << result[i] << " ";
        cout << endl;
    
    #endif
    }
    int_t poly_divat(const int_t* F, const int_t* G, int_t n, int_t k) {
        /*
            n次多项式F和G
            计算F(x)/G(x)的k次项前系数
            考虑F(x)*G(-x)/G(x)*G(-x),分母只有偶数次项,写为C(x^2);分子写成xA(x^2)+B(x^2),如果k是奇数,那么递归(A,C,n,k/2),如果k是偶数,那么递归(B,C,n,k/2)
            到k<=n时直接计算
        */
        int_t size2 = upper2n(2 * n + 1);
        int_t len = 1 << size2;
        static int_t fliparr[LARGE];
        makeflip(fliparr, size2);
        static int_t T1[LARGE], T2[LARGE], T3[LARGE];
        for (int_t i = 0; i < len; i++) {
            if (i <= n)
                T1[i] = F[i], T2[i] = G[i];
            else
                T1[i] = T2[i] = T3[i] = 0;
        }
        const int_t inv = power(len, -1);
        while (k >= n) {
    #ifdef DEBUG
            cout << "curr at k = " << k << endl;
    #endif
            for (int_t i = 0; i < len; i++) {
                if (i <= n) {
                    T3[i] = T2[i] * (i % 2 ? (mod - 1) : 1);
                } else
                    T3[i] = 0;
            }
            transform(T1, size2, fliparr);
            transform(T2, size2, fliparr);
            transform(T3, size2, fliparr);
            for (int_t i = 0; i < len; i++) {
                T1[i] = T1[i] * T3[i] % mod;
                T2[i] = T2[i] * T3[i] % mod;
            }
            transform<-1>(T1, size2, fliparr);
            transform<-1>(T2, size2, fliparr);
    #ifdef DEBUG
            cout << "prod T1 = ";
            for (int_t i = 0; i < len; i++)
                cout << T1[i] * inv % mod << " ";
            cout << endl;
            cout << "prod T2 = ";
            for (int_t i = 0; i < len; i++)
                cout << T2[i] * inv % mod << " ";
            cout << endl;
    
    #endif
            for (int_t i = 0; i < len; i++) {
                if (i * 2 < len) {
                    T2[i] = T2[i * 2] * inv % mod;
                } else
                    T2[i] = 0;
            }
            int_t b = k % 2;
            for (int_t i = 0; i < len; i++) {
                if (i % 2 == b) {
                    T1[i / 2] = T1[i] * inv % mod;
    #ifdef DEBUG
                    cout << "at " << i << " assgin " << T1[i] * inv << " to "
                         << i / 2 << endl;
    #endif
                }
    #ifdef DEBUG
                cout << "assign 0 to " << i << endl;
    #endif
                if (i > 0)
                    T1[i] = 0;
            }
    #ifdef DEBUG
            cout << "finished T1 = ";
            for (int_t i = 0; i < len; i++)
                cout << T1[i] % mod << " ";
            cout << endl;
            cout << "finished T2 = ";
            for (int_t i = 0; i < len; i++)
                cout << T2[i] % mod << " ";
            cout << endl;
    
    #endif
            k >>= 1;
        }
    
        poly_inv(T2, k + 1, T3);
        // #ifdef DEBUG
    
        //     cout << "finished k = " << k << endl;
        //     cout << "T1 = ";
        //     for (int_t i = 0; i < len; i++)
        //         cout << T1[i] % mod << " ";
        //     cout << endl;
        //     cout << "T2 = ";
        //     for (int_t i = 0; i < len; i++)
        //         cout << T2[i] % mod << " ";
        //     cout << endl;
        //     cout << "T2 inv = ";
        //     for (int_t i = 0; i < len; i++)
        //         cout << T3[i] % mod << " ";
        //     cout << endl;
    
        // #endif
        int_t result = 0;
        //计算结果的k次项
        for (int_t i = 0; i <= k; i++) {
            result = (result + T1[i] * T3[k - i] % mod) % mod;
        }
        return result;
    }
    void poly_mul(const int_t* A, int_t n, const int_t* B, int_t m, int_t* C) {
        /*
           计算n次多项式A与m次多项式B的乘积
        */
        int_t size2 = upper2n(n + m + 1);
        int_t len = 1 << size2;
        static int_t T1[LARGE], T2[LARGE];
        for (int_t i = 0; i < len; i++) {
            if (i <= n)
                T1[i] = A[i];
            else
                T1[i] = 0;
            if (i <= m)
                T2[i] = B[i];
            else
                T2[i] = 0;
        }
        static int_t fliparr[LARGE];
        makeflip(fliparr, size2);
        transform(T1, size2, fliparr);
        transform(T2, size2, fliparr);
        for (int_t i = 0; i < len; i++)
            T1[i] = T1[i] * T2[i] % mod;
        transform<-1>(T1, size2, fliparr);
        int_t inv = power(len, -1);
        for (int_t i = 0; i <= n + m; i++)
            C[i] = T1[i] * inv % mod;
    }
    int_t poly_linear_rec(const int_t* A0, const int_t* F0, int_t n, int_t k) {
        /*
            计算线性递推
            F[1],F[2]...F[k] 递推系数
            A[0],A[1]...A[k-1] 首项
    
            A[m]=A[m-1]*F[1]+A[m-2]*F[2]+...+A[m-k]*F[k]
        */
    
        static int_t T1[LARGE], T2[LARGE];
        static int_t Ax[LARGE], Fx[LARGE];
        Fx[0] = 0;
        Ax[k] = 0;
        for (int i = 1; i <= k; i++) {
            Fx[i] = F0[i];
            Ax[i - 1] = A0[i - 1];
        }
        poly_mul(Ax, k, Fx, k, T1);
        T1[0] = Ax[0];
        for (int i = 1; i <= k - 1; i++) {
            T1[i] = (Ax[i] - T1[i] + mod) % mod;
        }
        for (int i = k; i <= 2 * k; i++)
            T1[i] = 0;
        T2[0] = 1;
        for (int i = 1; i <= k; i++)
            T2[i] = (mod - Fx[i]) % mod;
    
        // #ifdef DEBUG
        //     cout << "T1 = ";
        //     for (int_t i = 0; i <= k; i++) {
        //         cout << T1[i] << " ";
        //     }
        //     cout << endl;
        //     cout << "T2 = ";
        //     for (int_t i = 0; i <= k; i++) {
        //         cout << T2[i] << " ";
        //     }
        //     cout << endl;
    
        // #endif
        return poly_divat(T1, T2, k, n);
    }
    
    void poly_log(const int_t* A, int_t n, int_t* out) {
        /*
            计算log(A(x)), A(x)为n次多项式
            DlogF(x) =DF(x)/F(x)
        */
        static int_t Ad[LARGE];
        static int_t Finv[LARGE];
        static int_t R[LARGE];
        const int_t m = n - 1;
        for (int_t i = 0; i <= m; i++) {
            Ad[i] = A[i + 1] * (i + 1) % mod;
        }
        Ad[n] = 0;
        poly_inv(A, n + 1, Finv);
        poly_mul(Ad, m, Finv, n, R);
        for (int_t i = 1; i <= n; i++) {
            out[i] = R[i - 1] * power(i, -1) % mod;
        }
    }
    void poly_exp(const int_t* A, int_t n, int_t* out) {
        /*
            计算exp(A(x)), A(x)为n次多项式
            H(x)=G(x)(1-logG(x)+A(x))
            H(x)为当次递推项,G(x)为上一次递推项
        */
        if (n == 1) {
            out[0] = 1;
            return;
        }
        int_t r = n / 2 + n % 2;
        poly_exp(A, r, out);
        static int_t G[LARGE];
        static int_t G2[LARGE];
        static int_t logG[LARGE];
        static int_t R[LARGE];
        for (int_t i = 0; i < r; i++)
            G[i] = out[i];
        for (int_t i = r; i < n; i++)
            G[i] = 0;
        poly_log(G, n - 1, logG);
        // for (int_t i = r; i < n; i++)
        //     logG[i] = 0;
        for (int_t i = 0; i < n; i++) {
            G2[i] = (mod - logG[i] + A[i]) % mod;
        }
        G2[0] = (G2[0] + 1) % mod;
        poly_mul(G, n - 1, G2, n - 1, R);
        for (int_t i = 0; i < n; i++)
            out[i] = R[i];
    #ifdef DEBUG
        cout << "at " << n << endl;
        cout << "A = ";
        for (int_t i = 0; i < n; i++)
            cout << A[i] << " ";
        cout << endl;
        cout << "G = ";
        for (int_t i = 0; i < n; i++)
            cout << G[i] << " ";
        cout << endl;
        cout << "logG = ";
        for (int_t i = 0; i < n; i++)
            cout << logG[i] << " ";
        cout << endl;
        cout << "G2 = ";
        for (int_t i = 0; i < n; i++)
            cout << G2[i] << " ";
        cout << endl;
        cout << "R = ";
        for (int_t i = 0; i < n; i++)
            cout << R[i] << " ";
        cout << endl;
        cout << endl;
    #endif
    }
    int main() {
        std::ios::sync_with_stdio(false);
        cin.tie(0), cout.tie(0);
    
        int_t n;
        cin >> n;
        static int_t A[LARGE], B[LARGE];
        static int_t C[LARGE];
        for (int_t i = 0; i < n; i++)
            cin >> A[i];
        poly_exp(A, n, B);
        for (int_t i = 0; i < n; i++)
            cout << B[i] << " ";
    
        return 0;
    }

     

  • 常系数线性递推的新做法 – 计算[x^k]P(x)/Q(x)

    $$ a_n=\sum_{1\le i\le k}{f_ia_{n-i}} \\ \left\{ a_0,a_1….a_{k-1} \right\} \text{已知} \\ \\ \text{设}F\left( x \right) \text{表示从第}k\text{项开始的该数列的生成函数} \\ \sum_{i\ge k}{a_ix^i}=\sum_{i\ge k}{\sum_{1\le j\le k}{f_ja_{i-j}x^i}} \\ =\sum_{1\le j\le k}{f_j\sum_{i\ge k}{a_{i-j}x^i}} \\ =\sum_{1\le j\le k}{f_j\sum_{i\ge k-j}{a_ix^{i+j}}} \\ =\sum_{1\le j\le k}{f_jx^j\sum_{i\ge k-j}{a_ix^i}} \\ =\sum_{1\le j\le k}{f_jx^j\left( F\left( x \right) +\sum_{k-j\le i\le k-1}{a_ix^i} \right)} \\ \\ F\left( x \right) =F\left( x \right) \sum_{1\le j\le k}{f_jx^j}+\sum_{1\le j\le k}{f_jx^j\sum_{k-j\le i\le k-1}{a_ix^i}} \\ F\left( x \right) =\frac{\sum_{1\le j\le k}{f_jx^j\sum_{k-j\le i\le k-1}{a_ix^i}}}{1-\sum_{1\le j\le k}{f_jx^j}} \\ =\frac{\sum_{1\le j\le k}{f_jx^j\sum_{k-j\le i-j\le k-1}{a_{i-j}x^{i-j}}}}{1-\sum_{1\le j\le k}{f_jx^j}} \\ =\frac{\sum_{1\le j\le k}{f_j\sum_{k\le i\le k-1+j}{a_{i-j}x^i}}}{1-\sum_{1\le j\le k}{f_jx^j}} \\ =\frac{\sum_{k\le i\le 2k-1}{\begin{array}{c} x^i\sum_{i-k+1\le j\le k}{a_{i-j}f_j}\\ \end{array}}}{1-\sum_{1\le j\le k}{f_jx^j}} \\ F\left( x \right) +\sum_{0\le i\le k-1}{a_ix^i}=\frac{\left( 1-\sum_{1\le j\le k}{f_jx^j} \right) \left( \sum_{0\le i\le k-1}{a_ix^i} \right) +\sum_{1\le j\le k}{f_j\sum_{k\le i\le k-1+j}{a_{i-j}x^i}}}{1-\sum_{1\le j\le k}{f_jx^j}} \\ \frac{\sum_{0\le i\le k-1}{a_ix^i}-\sum_{1\le j\le k}{f_j\sum_{0\le i\le k-1}{a_ix^{i+j}}}+\sum_{1\le j\le k}{f_j\sum_{k\le i\le k-1+j}{a_{i-j}x^i}}}{1-\sum_{1\le j\le k}{f_jx^j}} \\ =\frac{\sum_{0\le i\le k-1}{a_ix^i}-\sum_{1\le j\le k}{f_j\sum_{j\le i\le j+k-1}{a_{i-j}x^i}}+\sum_{1\le j\le k}{f_j\sum_{k\le i\le k-1+j}{a_{i-j}x^i}}}{1-\sum_{1\le j\le k}{f_jx^j}} \\ \frac{\sum_{0\le i\le k-1}{a_ix^i}+\sum_{1\le j\le k}{f_j\left( \sum_{k\le i\le k-1+j}{a_{i-j}x^i}-\sum_{j\le i\le j+k-1}{a_{i-j}x^i} \right)}}{1-\sum_{1\le j\le k}{f_jx^j}} \\ \frac{\sum_{0\le i\le k-1}{a_ix^i}+\sum_{1\le j\le k}{f_j\left( \sum_{k\le i\le k-1+j}{a_{i-j}x^i}-\sum_{k\le i\le j+k-1}{a_{i-j}x^i}-\sum_{j\le i\le k-1}{a_{i-j}x^i} \right)}}{1-\sum_{1\le j\le k}{f_jx^j}} \\ \frac{\sum_{0\le i\le k-1}{a_ix^i}+\sum_{1\le j\le k}{f_j\left( -\sum_{j\le i\le k-1}{a_{i-j}x^i} \right)}}{1-\sum_{1\le j\le k}{f_jx^j}} \\ =\frac{\sum_{0\le i\le k-1}{a_ix^i}-\sum_{1\le j\le k}{f_j\left( \sum_{j\le i\le k-1}{a_{i-j}x^i} \right)}}{1-\sum_{1\le j\le k}{f_jx^j}} \\ =\frac{\sum_{0\le i\le k-1}{a_ix^i}-\sum_{1\le j\le k-1}{f_j\left( \sum_{j\le i\le k-1}{a_{i-j}x^i} \right)}}{1-\sum_{1\le j\le k}{f_jx^j}} \\ =\frac{\sum_{0\le i\le k-1}{a_ix^i}-\sum_{1\le i\le k-1}{x^i\sum_{1\le j\le i}{a_{i-j}f_j}}}{1-\sum_{1\le j\le k}{f_jx^j}} \\ =\frac{a_0+\sum_{1\le i\le k-1}{x^ia_i}-\sum_{1\le i\le k-1}{x^i\sum_{1\le j\le i}{a_{i-j}f_j}}}{1-\sum_{1\le j\le k}{f_jx^j}} \\ =\frac{a_0+\sum_{1\le i\le k-1}{x^i\left( a_i-\sum_{1\le j\le i}{a_{i-j}f_j} \right)}}{1-\sum_{1\le j\le k}{f_jx^j}} \\ \text{我们得到了原数列}\left\{ a_i \right\} \text{的生成函数}G\left( x \right) =\frac{P\left( x \right)}{Q\left( x \right)} \\ \text{考虑计算这个有理函数的}k\text{次项} \\ \text{令}G_k\left( x \right) =\frac{P\left( x \right)}{Q\left( x \right)}=\frac{P\left( x \right) P\left( -x \right)}{Q\left( x \right) Q\left( -x \right)}=\frac{xA\left( x^2 \right) +B\left( x^2 \right)}{C\left( x^2 \right)} \\ \text{即分母只有偶数次方项},\text{分子的奇数次方项和偶数次方项拆开} \\ \text{如果}k\text{为奇数},\text{那么}x^k\text{之可能存在于}x\frac{A\left( x^2 \right)}{C\left( x^2 \right)}\text{中}\left( \text{只有这边才有奇数次方} \right) \\ \text{同理},\text{如果}k\text{为偶数,那么}x^k\text{只能出现在}\frac{B\left( x^2 \right)}{C\left( x^2 \right)}\text{中,} \\ \text{然后根据}k\text{的奇偶性,继续递归} \\ \text{如果}k\text{是奇数},\text{那么计算}G_{\frac{k-1}{2}}\left( x \right) =\frac{A\left( x \right)}{C\left( x \right)}\left( \text{项的次数除以}2 \right) ,\text{答案为}\left[ x^{\frac{k-1}{2}} \right] G_{\frac{k-1}{2}}\left( x \right) \\ \text{如果}k\text{是偶数},\text{那么计算}G_{\frac{k}{2}}\left( x \right) =\frac{B\left( x \right)}{C\left( x \right)}\text{,答案为}\left[ x^k \right] G_{\frac{k}{2}}\left( x \right) \\ \text{然后有两种选择}:\text{递归到}k=0\text{时,计算}\frac{P\left( 0 \right)}{G\left( 0 \right)}\mathrm{mod}p\text{,此时不需要多项式求逆运算} \\ \text{但常数较大}\left( \text{递归次数多} \right) \\ \text{递归到}k<n\text{时,直接求逆计算。} \\ \\ \\ \\ \\ $$

    代码1: 递归到$k=0$时执行整数运算,较慢

    #include <cstring>
    #include <iostream>
    using int_t = long long int;
    using std::cin;
    using std::cout;
    using std::endl;
    
    const int_t LARGE = 2e6;
    
    const int_t mod = 998244353;
    const int_t g = 3;
    int_t power(int_t b, int_t i) {
        int_t r = 1;
        if (i < 0)
            i = (i % (mod - 1) + mod - 1) % (mod - 1);
        while (i) {
            if (i & 1)
                r = r * b % mod;
            b = b * b % mod;
            i >>= 1;
        }
        return r;
    }
    
    void makeflip(int_t* arr, int_t size2) {
        int_t len = (1 << size2);
        arr[0] = 0;
        for (int_t i = 1; i < len; i++) {
            arr[i] = (arr[i >> 1] >> 1) | ((i & 1) << (size2 - 1));
        }
    }
    
    int_t upper2n(int_t x) {
        int_t r = 0;
        while ((1 << r) < x)
            r++;
        return r;
    }
    template <int_t arg = 1>
    void transform(int_t* A, int_t size2, int_t* flip) {
        int_t len = (1 << size2);
        for (int_t i = 0; i < len; i++) {
            int_t r = flip[i];
            if (r > i)
                std::swap(A[i], A[r]);
        }
        for (int_t i = 2; i <= len; i *= 2) {
            int_t mr = power(g, arg * (mod - 1) / i);
            for (int_t j = 0; j < len; j += i) {
                int_t curr = 1;
                for (int_t k = 0; k < i / 2; k++) {
                    int_t u = A[j + k], v = curr * A[j + k + i / 2] % mod;
                    A[j + k] = (u + v) % mod;
                    A[j + k + i / 2] = (u - v + mod) % mod;
                    curr = curr * mr % mod;
                }
            }
        }
    }
    void poly_inv(const int_t* A, int_t n, int_t* result) {
        /*
        计算B(x)*A(x) = 1 mod x^n, 其中A(x)已知
        假设已知A(x)*B(x) = 1 mod x^{ceil(n/2)}
        假设C(x)*A(x) = 1 mod x^n
        (A(x)B(x)-1)^2 = A^2(x)B^2(x)-2A(x)B(x)+1= 0
        A(x)B^2(x)-2B(x)+C(x) = 0
        C(x) = 2B(x)-A(x)B^2(x)
        */
        if (n == 1) {
            result[0] = power(A[0], -1);
            return;
        }
        int_t next = n / 2 + n % 2;
        poly_inv(A, next, result);
        //次数不要选错了,应该用n次的A和B去卷
        int_t size2 = upper2n(n + 2 * next + 1);
        static int_t X[LARGE];
        static int_t Y[LARGE];
        int_t len = (1 << size2);
        memset(X + next, 0, sizeof(X[0]) * (len - n));
        memset(Y + next, 0, sizeof(Y[0]) * (len - next));
        memcpy(X, A, sizeof(A[0]) * n);
        memcpy(Y, result, sizeof(result[0]) * next);
        static int_t fliparr[LARGE];
        makeflip(fliparr, size2);
        transform<1>(X, size2, fliparr);
        transform<1>(Y, size2, fliparr);
    
        for (int_t i = 0; i < len; i++) {
            X[i] = (2 * Y[i] - X[i] * Y[i] % mod * Y[i] % mod + mod) % mod;
        }
        transform<-1>(X, size2, fliparr);
        const int_t inv = power(len, -1);
        for (int_t i = 0; i < n; i++)
            result[i] = X[i] * inv % mod;
    }
    int_t poly_divat(const int_t* F, const int_t* G, int_t n, int_t k) {
        /*
            n次多项式F和G
            计算F(x)/G(x)的k次项前系数
            考虑F(x)*G(-x)/G(x)*G(-x),分母只有偶数次项,写为C(x^2);分子写成xA(x^2)+B(x^2),如果k是奇数,那么递归(A,C,n,k/2),如果k是偶数,那么递归(B,C,n,k/2)
            到k<=n时直接计算
        */
        int_t size2 = upper2n(2 * n + 1);
        int_t len = 1 << size2;
        static int_t fliparr[LARGE];
        makeflip(fliparr, size2);
        static int_t T1[LARGE], T2[LARGE], T3[LARGE];
        for (int_t i = 0; i < len; i++) {
            if (i <= n)
                T1[i] = F[i], T2[i] = G[i];
            else
                T1[i] = T2[i] = T3[i] = 0;
        }
        const int_t inv = power(len, -1);
        while (k != 0) {
            for (int_t i = 0; i < len; i++) {
                if (i <= n) {
                    T3[i] = T2[i] * (i % 2 ? (mod - 1) : 1);
                } else
                    T3[i] = 0;
            }
            transform(T1, size2, fliparr);
            transform(T2, size2, fliparr);
            transform(T3, size2, fliparr);
            for (int_t i = 0; i < len; i++) {
                T1[i] = T1[i] * T3[i] % mod;
                T2[i] = T2[i] * T3[i] % mod;
            }
            transform<-1>(T1, size2, fliparr);
            transform<-1>(T2, size2, fliparr);
            for (int_t i = 0; i < len; i++) {
                if (i * 2 < len) {
                    T2[i] = T2[i * 2] * inv % mod;
                } else
                    T2[i] = 0;
            }
            int_t b = k % 2;
            for (int_t i = 0; i < len; i++) {
                if (i % 2 == b) {
                    T1[i / 2] = T1[i] * inv % mod;
                }
    
                if (i > 0)  //防止把T1[0]改为0
                    T1[i] = 0;
            }
            k >>= 1;
        }
    
        return T1[0] * power(T2[0], -1) % mod;
    }
    void poly_mul(const int_t* A, int_t n, const int_t* B, int_t m, int_t* C) {
        int_t size2 = upper2n(n + m + 1);
        int_t len = 1 << size2;
        static int_t T1[LARGE], T2[LARGE];
        for (int_t i = 0; i < len; i++) {
            if (i <= n)
                T1[i] = A[i];
            else
                T1[i] = 0;
            if (i <= m)
                T2[i] = B[i];
            else
                T2[i] = 0;
        }
        static int_t fliparr[LARGE];
        makeflip(fliparr, size2);
        transform(T1, size2, fliparr);
        transform(T2, size2, fliparr);
        for (int_t i = 0; i < len; i++)
            T1[i] = T1[i] * T2[i] % mod;
        transform<-1>(T1, size2, fliparr);
        int_t inv = power(len, -1);
        for (int_t i = 0; i <= n + m; i++)
            C[i] = T1[i] * inv % mod;
    }
    int main() {
        std::ios::sync_with_stdio(false);
        static int_t A[LARGE], F[LARGE];
        static int_t T1[LARGE], T2[LARGE];
        int_t n, k;
        cin >> n >> k;
        for (int_t i = 1; i <= k; i++) {
            cin >> F[i];
            F[i] = (F[i] % mod + mod) % mod;
        }
        F[0] = 0;
        for (int_t i = 0; i < k; i++) {
            cin >> A[i];
            A[i] = (A[i] % mod + mod) % mod;
        }
        A[k] = 0;
        poly_mul(A, k, F, k, T1);
        T1[0] = A[0];
        for (int_t i = 1; i <= k - 1; i++) {
            T1[i] = (A[i] - T1[i] + mod) % mod;
        }
        for (int_t i = k; i <= 2 * k; i++)
            T1[i] = 0;
        T2[0] = 1;
        for (int_t i = 1; i <= k; i++)
            T2[i] = (mod - F[i]) % mod;
        int_t r = poly_divat(T1, T2, k, n);
        cout << r << endl;
        return 0;
    }
    /*
    (1+2*x)/(1+x+x^2)
    2 5
    
    1 2 0
    1 1 1
    
    
    ans = -2 = 998244351
    
    */

    代码2: 递归到$k<n$时直接求逆运算,较快,大约是代码1的一半时间

    #include <cstring>
    #include <iostream>
    using int_t = long long int;
    using std::cin;
    using std::cout;
    using std::endl;
    
    const int_t LARGE = 2e6;
    
    const int_t mod = 998244353;
    const int_t g = 3;
    int_t power(int_t b, int_t i) {
        int_t r = 1;
        if (i < 0)
            i = (i % (mod - 1) + mod - 1) % (mod - 1);
        while (i) {
            if (i & 1)
                r = r * b % mod;
            b = b * b % mod;
            i >>= 1;
        }
        return r;
    }
    
    void makeflip(int_t* arr, int_t size2) {
        int_t len = (1 << size2);
        arr[0] = 0;
        for (int_t i = 1; i < len; i++) {
            arr[i] = (arr[i >> 1] >> 1) | ((i & 1) << (size2 - 1));
        }
    }
    
    int_t upper2n(int_t x) {
        int_t r = 0;
        while ((1 << r) < x)
            r++;
        return r;
    }
    template <int_t arg = 1>
    void transform(int_t* A, int_t size2, int_t* flip) {
        int_t len = (1 << size2);
        for (int_t i = 0; i < len; i++) {
            int_t r = flip[i];
            if (r > i)
                std::swap(A[i], A[r]);
        }
        for (int_t i = 2; i <= len; i *= 2) {
            int_t mr = power(g, arg * (mod - 1) / i);
            for (int_t j = 0; j < len; j += i) {
                int_t curr = 1;
                for (int_t k = 0; k < i / 2; k++) {
                    int_t u = A[j + k], v = curr * A[j + k + i / 2] % mod;
                    A[j + k] = (u + v) % mod;
                    A[j + k + i / 2] = (u - v + mod) % mod;
                    curr = curr * mr % mod;
                }
            }
        }
    }
    void poly_inv(const int_t* A, int_t n, int_t* result) {
        /*
        计算B(x)*A(x) = 1 mod x^n, 其中A(x)已知
        假设已知A(x)*B(x) = 1 mod x^{ceil(n/2)}
        假设C(x)*A(x) = 1 mod x^n
        (A(x)B(x)-1)^2 = A^2(x)B^2(x)-2A(x)B(x)+1= 0
        A(x)B^2(x)-2B(x)+C(x) = 0
        C(x) = 2B(x)-A(x)B^2(x)
        */
        if (n == 1) {
            result[0] = power(A[0], -1);
            return;
        }
        int_t next = n / 2 + n % 2;
        poly_inv(A, next, result);
        //次数不要选错了,应该用n次的A和B去卷
        int_t size2 = upper2n(n + 2 * next + 1);
        static int_t X[LARGE];
        static int_t Y[LARGE];
        int_t len = (1 << size2);
        memset(X + next, 0, sizeof(X[0]) * (len - n));
        memset(Y + next, 0, sizeof(Y[0]) * (len - next));
        memcpy(X, A, sizeof(A[0]) * n);
        memcpy(Y, result, sizeof(result[0]) * next);
        static int_t fliparr[LARGE];
        makeflip(fliparr, size2);
        transform<1>(X, size2, fliparr);
        transform<1>(Y, size2, fliparr);
    
        for (int_t i = 0; i < len; i++) {
            X[i] = (2 * Y[i] - X[i] * Y[i] % mod * Y[i] % mod + mod) % mod;
        }
        transform<-1>(X, size2, fliparr);
        const int_t inv = power(len, -1);
        for (int_t i = 0; i < n; i++)
            result[i] = X[i] * inv % mod;
    }
    int_t poly_divat(const int_t* F, const int_t* G, int_t n, int_t k) {
        /*
            n次多项式F和G
            计算F(x)/G(x)的k次项前系数
            考虑F(x)*G(-x)/G(x)*G(-x),分母只有偶数次项,写为C(x^2);分子写成xA(x^2)+B(x^2),如果k是奇数,那么递归(A,C,n,k/2),如果k是偶数,那么递归(B,C,n,k/2)
            到k<=n时直接计算
        */
        int_t size2 = upper2n(2 * n + 1);
        int_t len = 1 << size2;
        static int_t fliparr[LARGE];
        makeflip(fliparr, size2);
        static int_t T1[LARGE], T2[LARGE], T3[LARGE];
        for (int_t i = 0; i < len; i++) {
            if (i <= n)
                T1[i] = F[i], T2[i] = G[i];
            else
                T1[i] = T2[i] = T3[i] = 0;
        }
        const int_t inv = power(len, -1);
        while (k >= n) {
            for (int_t i = 0; i < len; i++) {
                if (i <= n) {
                    T3[i] = T2[i] * (i % 2 ? (mod - 1) : 1);
                } else
                    T3[i] = 0;
            }
            transform(T1, size2, fliparr);
            transform(T2, size2, fliparr);
            transform(T3, size2, fliparr);
            for (int_t i = 0; i < len; i++) {
                T1[i] = T1[i] * T3[i] % mod;
                T2[i] = T2[i] * T3[i] % mod;
            }
            transform<-1>(T1, size2, fliparr);
            transform<-1>(T2, size2, fliparr);
            for (int_t i = 0; i < len; i++) {
                if (i * 2 < len) {
                    T2[i] = T2[i * 2] * inv % mod;
                } else
                    T2[i] = 0;
            }
            int_t b = k % 2;
            for (int_t i = 0; i < len; i++) {
                if (i % 2 == b) {
                    T1[i / 2] = T1[i] * inv % mod;
                }
    
                if (i > 0)  //防止把T1[0]改为0
                    T1[i] = 0;
            }
            k >>= 1;
        }
    
        poly_inv(T2, k + 1, T3);
        int_t result = 0;
        //计算结果的k次项
        for (int_t i = 0; i <= k; i++) {
            result = (result + T1[i] * T3[k - i] % mod) % mod;
        }
        return result;
    }
    void poly_mul(const int_t* A, int_t n, const int_t* B, int_t m, int_t* C) {
        int_t size2 = upper2n(n + m + 1);
        int_t len = 1 << size2;
        static int_t T1[LARGE], T2[LARGE];
        for (int_t i = 0; i < len; i++) {
            if (i <= n)
                T1[i] = A[i];
            else
                T1[i] = 0;
            if (i <= m)
                T2[i] = B[i];
            else
                T2[i] = 0;
        }
        static int_t fliparr[LARGE];
        makeflip(fliparr, size2);
        transform(T1, size2, fliparr);
        transform(T2, size2, fliparr);
        for (int_t i = 0; i < len; i++)
            T1[i] = T1[i] * T2[i] % mod;
        transform<-1>(T1, size2, fliparr);
        int_t inv = power(len, -1);
        for (int_t i = 0; i <= n + m; i++)
            C[i] = T1[i] * inv % mod;
    }
    int main() {
        std::ios::sync_with_stdio(false);
        static int_t A[LARGE], F[LARGE];
        static int_t T1[LARGE], T2[LARGE];
        int_t n, k;
        cin >> n >> k;
        for (int_t i = 1; i <= k; i++) {
            cin >> F[i];
            F[i] = (F[i] % mod + mod) % mod;
        }
        F[0] = 0;
        for (int_t i = 0; i < k; i++) {
            cin >> A[i];
            A[i] = (A[i] % mod + mod) % mod;
        }
        A[k] = 0;
        poly_mul(A, k, F, k, T1);
        T1[0] = A[0];
        for (int_t i = 1; i <= k - 1; i++) {
            T1[i] = (A[i] - T1[i] + mod) % mod;
        }
        for (int_t i = k; i <= 2 * k; i++)
            T1[i] = 0;
        T2[0] = 1;
        for (int_t i = 1; i <= k; i++)
            T2[i] = (mod - F[i]) % mod;
        int_t r = poly_divat(T1, T2, k, n);
        cout << r << endl;
        return 0;
    }
    /*
    (1+2*x)/(1+x+x^2)
    2 5
    
    1 2 0
    1 1 1
    
    
    ans = -2 = 998244351
    
    */

     

  • CFGYM 102823 CCPC2018 桂林 B题 Array Modify

    题倒是还行,达标找规律找出来结论,然后就想到了一个多项式快速幂的做法。

    一开始写的两个log的快速幂,TLE。

    然后尝试改成exp+log,跑得看起来快了,但是遇到多组数据(n的大小递增的时候)就挂掉了?

    仔细调查后发现源于自己一直以来exp函数内一直存在的错误:

    exp内调用了log,log内调用了inv,exp内给G0(用以存储log返回值的数组)预留了upper2n(2*n)的空间,但是inv内却使用了upper2n(3*n+1)的空间!

    (调试方法:数组切换为动态分配并开启内存检查)

    对于只计算exp一次而言,这个错误不会造成影响,但是计算多次的时候由于空间污染会导致错误结果。

    解决方法:在exp内也预留upper2n(3*n)的空间,通过本题。

    另外小测试了一下,两个log的多项式快速幂跑本题极限数据,NTT变换的多项式总长度为“27262976”,而exp+log实现的快速幂,NTT变换的多项式总长度为 19922408,我原本以为因为常数原因,exp+log会跑得更慢,没想到事实并非如此。

     

    // luogu-judger-enable-o2
    #include 
    #include 
    #include 
    #include 
    #include 
    #include 
    #include 
    using int_t = long long int;
    using std::cin;
    using std::cout;
    using std::endl;
    const int mod = 998244353;
    const int g = 3;
    const int LARGE = 1 << 20;
    int revs[21][LARGE + 1];
    
    int power(int base, int index);
    void transform(int* A, int len, int arg);
    void poly_inv(const int* A, int n, int* result);
    void poly_log(const int* A, int n, int* result);
    void poly_exp(const int* A, int n, int* result);
    
    std::vector transform_count;
    #ifdef TDEBUG
    #define TRANSFORM_DEBUG(n) \
        { transform_count.push_back(n); }
    #else
    #define TRANSFORM_DEBUG(n)
    #endif
    int bitRev(int base, int size2) {
        int result = 0;
        for (int i = 1; i < size2; i++) {
            result |= (base & 1);
            base >>= 1;
            result <<= 1;
        }
        result |= (base & 1);
        return result;
    }
    int upper2n(int x) {
        int result = 1;
        while (result < x)
            result *= 2;
        return result;
    }
    int main() {
        std::ios::sync_with_stdio(false);
    
        for (int i = 0; (1 << i) <= LARGE; i++) {
            for (int j = 0; j < LARGE; j++) {
                revs[i][j] = bitRev(j, i);
            }
        }
        static int F[LARGE + 1], G[LARGE + 1];
        static int arr[LARGE + 1];
        int T;
        cin >> T;
        for (int_t _i = 1; _i <= T; _i++) {
            memset(F, 0, sizeof(F));
            memset(G, 0, sizeof(G));
            memset(arr, 0, sizeof(arr));
            int_t n, L, m;
            cin >> n >> L >> m;
    
            for (int_t i = 1; i <= n; i++)
                cin >> arr[i];
            int size2 = upper2n(4 * (n + 1));
            for (int i = 0; i < size2; i++) {
                if (i < L)
                    F[i] = 1;
                else
                    F[i] = 0;
                G[i] = 0;
            }
    #ifdef DEBUG
            cout << "init = " << endl;
            cout << "F = ";
            for (int_t i = 0; i < 2 * n; i++)
                cout << F[i] << " ";
            cout << endl;
    
            cout << "G = ";
            for (int_t i = 0; i < 2 * n; i++)
                cout << G[i] << " ";
            cout << endl;
    #endif
    
            poly_log(F, n, G);
    #ifdef DEBUG
            cout << "after log = " << endl;
            cout << "G = ";
            for (int_t i = 0; i < 2 * n; i++)
                cout << G[i] << " ";
            cout << endl;
    #endif
            for (int i = 0; i < n; i++)
                G[i] = G[i] * m % mod;
            for (int i = 0; i < size2; i++)
                F[i] = 0;
    #ifdef DEBUG
            cout << "before exp = " << endl;
            cout << "G = ";
            for (int_t i = 0; i < 2 * n; i++)
                cout << G[i] << " ";
            cout << endl;
    #endif
            poly_exp(G, n, F);
    #ifdef DEBUG
            cout << "after exp = " << endl;
            cout << "F = ";
            for (int_t i = 0; i < 2 * n; i++)
                cout << F[i] << " ";
            cout << endl;
    #endif
    
            for (int i = 0; i < size2; i++) {
                if (i >= n)
                    F[i] = 0;
                if (i >= n && i < 2 * n)
                    G[i] = arr[i - n + 1];
                else
                    G[i] = 0;
            }
            std::reverse(G + n, G + 2 * n);
    
    #ifdef DEBUG
            cout << "F = ";
            for (int_t i = 0; i < 2 * n; i++)
                cout << F[i] << " ";
            cout << endl;
    
            cout << "G = ";
            for (int_t i = 0; i < 2 * n; i++)
                cout << G[i] << " ";
            cout << endl;
    #endif
    
            transform(F, size2, 1);
            transform(G, size2, 1);
            for (int i = 0; i < size2; i++)
                F[i] = (int_t)F[i] * G[i] % mod;
            transform(F, size2, -1);
            int inv = power(size2, -1);
            std::reverse(F + n, F + 2 * n);
            cout << "Case " << _i << ": ";
    
            for (int i = n; i < 2 * n; i++) {
                cout << (int_t)F[i] * inv % mod << " ";
            }
            cout << endl;
    #ifdef DEBUG
            for (int i = 0; i < n; i++) {
                cout << "F " << i << " = " << F[i] << endl;
            }
    #endif
    #ifdef TDEBUG
            for (auto x : transform_count) {
                cout << x << " ";
            }
            cout << endl;
    #endif
        }
        return 0;
    }
    
    int power(int base, int index) {
        const int phi = mod - 1;
        index = (index % phi + phi) % phi;
        base = (base % mod + mod) % mod;
        int result = 1;
        while (index) {
            if (index & 1)
                result = (int_t)result * base % mod;
            base = (int_t)base * base % mod;
            index >>= 1;
        }
        return result;
    }
    void transform(int* A, int len, int arg) {
        TRANSFORM_DEBUG(len);
        int size2 = int(log2(len) + 0.5);
        for (int i = 0; i < len; i++) {
            int x = revs[size2][i];
            if (x > i)
                std::swap(A[i], A[x]);
        }
        for (int i = 2; i <= len; i *= 2) {
            int mr = power(g, (int_t)arg * (mod - 1) / i);
            for (int j = 0; j < len; j += i) {
                int curr = 1;
                for (int k = 0; k < i / 2; k++) {
                    int u = A[j + k];
                    int t = (int_t)A[j + k + i / 2] * curr % mod;
                    A[j + k] = (u + t) % mod;
                    A[j + k + i / 2] = (u - t + mod) % mod;
                    curr = (int_t)curr * mr % mod;
                }
            }
        }
    }
    //计算多项式A在模x^n下的逆元
    // C(x)<-2B(x)-A(x)B^2(x)
    void poly_inv(const int* A, int n, int* result) {
        if (n == 1) {
            result[0] = power(A[0], -1);
            return;
        }
        poly_inv(A, n / 2 + n % 2, result);
        static int temp[LARGE + 1];
        int size2 = upper2n(3 * n + 1);
        for (int i = 0; i < size2; i++) {
            if (i < n)
                temp[i] = A[i];
            else
                temp[i] = 0;
        }
        transform(temp, size2, 1);
        transform(result, size2, 1);
        for (int i = 0; i < size2; i++) {
            result[i] =
                ((int_t)2 * result[i] % mod -
                 (int_t)temp[i] * result[i] % mod * result[i] % mod + 2 * mod) %
                mod;
        }
        transform(result, size2, -1);
        for (int i = 0; i < size2; i++) {
            if (i < n) {
                result[i] = (int_t)result[i] * power(size2, -1) % mod;
            } else {
                result[i] = 0;
            }
        }
    }
    void poly_log(const int* A, int n, int* result) {
        int size2 = upper2n(3 * n + 1);
        static int Ad[LARGE + 1];
        for (int i = 0; i < size2; i++) {
            if (i < n) {
                Ad[i] = (int_t)(i + 1) * A[i + 1] % mod;
            } else {
                Ad[i] = 0;
            }
            result[i] = 0;
        }
        transform(Ad, size2, 1);
        poly_inv(A, n, result);
        transform(result, size2, 1);
        for (int i = 0; i < size2; i++) {
            result[i] = (int_t)result[i] * Ad[i] % mod;
        }
        transform(result, size2, -1);
        for (int i = 0; i < size2; i++)
            if (i < n)
                result[i] = (int_t)result[i] * power(size2, -1) % mod;
            else
                result[i] = 0;
        for (int i = n - 1; i >= 1; i--) {
            result[i] = (int_t)result[i - 1] * power(i, -1) % mod;
        }
        result[0] = 0;
    }
    void poly_exp(const int* A, int n, int* result) {
        if (n == 1) {
            assert(A[0] == 0);
            result[0] = 1;
            return;
        }
        poly_exp(A, n / 2 + n % 2, result);
        static int Ax[LARGE + 1];
        // memset(G0, 0, sizeof(G0));
        // memset(Ax,0,sizeof(Ax));
        int size2 = upper2n(3 * n + 1);
        // int* G0 = new int[size2];
        static int G0[LARGE + 1];
        for (int_t i = 0; i < size2; i++) {
            if (i < n)
                Ax[i] = A[i];
            else
                Ax[i] = 0;
            G0[i] = 0;
        }
        poly_log(result, n, G0);
    
        transform(Ax, size2, 1);
        transform(G0, size2, 1);
        transform(result, size2, 1);
        for (int i = 0; i < size2; i++)
            result[i] =
                (result[i] - (int_t)result[i] * (G0[i] - Ax[i] + mod) % mod + mod) %
                mod;
        transform(result, size2, -1);
        for (int i = 0; i < size2; i++)
            if (i < n)
                result[i] = (int_t)result[i] * power(size2, -1) % mod;
            else
                result[i] = 0;
        // delete[] G0;
    }
    
  • noi.ac 154 Curiosity

    $$ \text{设}f\left( i,j \right) \text{表示第}i\text{次扔骰子是否取到}j\text{这个数,显然这个东西的期望} \\ \text{是}\frac{1}{K}\left( \text{一次扔骰子只能出现一个数} \right) \\ \text{所以}a_i=\sum_{1\le k\le n}{f\left( k,i \right)} \\ \text{所以}E\left( \prod_{1\le i\le L}{a_i^F} \right) =\prod_{1\le i\le L}{\left( \sum_{1\le j\le N}{f\left( j,i \right)} \right) ^F} \\ \text{考虑}\prod_{1\le i\le L}{\left( \sum_{1\le j\le N}{f\left( j,i \right)} \right) ^F}\text{展开后的一项,如果这一项中存在}f\left( i,a \right) \text{和}f\left( i,b \right) ,\text{并且}a\ne b, \\ \text{那么这一整项的期望必须是}0\left( \text{否则出现了不合法情况} \right) \\ \text{对于其他的项,比如}f\left( \text{1,}a \right) ^2f\left( \text{2,}b \right) ^2….,\text{假设这一项由}p\text{个不同的变量构成} \\ \text{那么这一项的期望为}\frac{1}{K^p}\left( \text{独立变量期望可乘} \right) \\ \text{现在我们的问题在于如何计算出}\prod_{1\le i\le L}{\left( \sum_{1\le j\le N}{f\left( j,i \right)} \right) ^F}\text{中由}p\text{个独立变量构成的项的个数。} \\ \text{首先考虑}\left( \sum_{1\le j\le N}{f\left( j,i \right)} \right) ^F,\text{显然这个式子展开后不存在期望为0的项。} \\ \,\,\text{我们忽略项中是由哪些变量}\left( \text{比如}f\left( \text{1,}i \right) ,f\left( \text{2,}i \right) \text{等等} \right) \text{构成的,也忽略掉项之间的顺序关系} \\ \left( \begin{array}{c} \text{在}\left( f\left( \text{1,}i \right) +f\left( \text{2,}i \right) +f\left( \text{3,}i \right) .. \right) \times \left( f\left( \text{1,}i \right) +f\left( \text{2,}i \right) +f\left( \text{3,}i \right) .. \right) \times \left( f\left( \text{1,}i \right) +f\left( \text{2,}i \right) +f\left( \text{3,}i \right) .. \right) \text{中,}f\left( \text{1,}i \right) ^2\times f\left( \text{2,}i \right)\\ \text{会被计算6次}\\ \end{array} \right) \\ \text{则}\left( \sum_{1\le j\le N}{f\left( j,i \right)} \right) ^F\text{展开后忽略上述限制存在}i\text{个无关变量的项的系数和为}S_2\left( F,i \right) ,\text{其中}S_2\text{为第二类斯特林数。} \\ \text{则}S_2\left( F,i \right) \times \left( \begin{array}{c} F\\ i\\ \end{array} \right) \times i!\text{为不忽略限制的符合条件的项的系数和。} \\ \text{构造生成函数}F\left( x \right) =\sum_{i\ge 0}{S_2\left( F,i \right) x^i} \\ \text{则}F\left( x \right) ^L\text{展开后}i\text{次项前的系数即为忽略上述限制后,存在}i\text{个无关变量的项的系数和。} \\ \text{为了还原至原式,将这个系数乘上}\left( \begin{array}{c} n\\ i\\ \end{array} \right) i!\text{后即为}\prod_{1\le i\le L}{\left( \sum_{1\le j\le N}{f\left( j,i \right)} \right) ^F}\text{展开后存在}i\text{个无关变量的项的系数和。} \\ \text{显然只有不超过}L\times F\text{个无关变量的项才是有意义的}. \\ \text{复杂度}O\left( LF\log LF\times \log L \right) \\ \text{注意}F=\text{1的时候需要特判,否则太慢了跑不过去。} \\ $$

    #include <assert.h>
    #include <algorithm>
    #include <cmath>
    // #include <complex>
    #include <cstring>
    #include <iostream>
    using real_t = double;
    
    using int_t = long long int;
    using std::cin;
    using std::cout;
    using std::endl;
    using cpx_t = struct complex;
    
    struct complex {
        real_t real, imag;
        complex(real_t real = 0, real_t imag = 0) {
            this->real = real;
            this->imag = imag;
        }
        complex operator+(const complex& rhs) const {
            return complex{real + rhs.real, imag + rhs.imag};
        }
        complex operator-(const complex& rhs) const {
            return complex{real - rhs.real, imag - rhs.imag};
        }
        complex operator*(const complex& rhs) const {
            return complex{real * rhs.real - imag * rhs.imag,
                           real * rhs.imag + imag * rhs.real};
        }
        complex& operator*=(const complex& rhs) {
            *this = (*this) * rhs;
            return *this;
        }
    };
    
    const int_t LARGE = 300000;
    const int_t mod = 2333;
    int_t S2[mod][mod];
    int_t fact[LARGE + 1], factInv[LARGE + 1], inv[LARGE + 1];
    const real_t PI = acosl(-1);
    int_t flips[LARGE + 1];
    void transform(cpx_t* A, int_t size2, int_t arg) {
        for (int_t i = 0; i < (1 << size2); i++) {
            int_t x = flips[i];
            if (x > i) std::swap(A[i], A[x]);
        }
        for (int i = 2; i <= (1 << size2); i *= 2) {
            auto mr = cpx_t(cos(2 * PI / i), arg * sin(2 * PI / i));
            for (int j = 0; j < (1 << size2); j += i) {
                auto curr = cpx_t(1, 0);
                for (int k = 0; k < i / 2; k++) {
                    auto u = A[j + k], t = A[j + k + i / 2] * curr;
                    A[j + k] = u + t;
                    A[j + k + i / 2] = u - t;
                    curr *= mr;
                }
            }
        }
    }
    void poly_mul(const int_t* Ax, const int_t* Bx, int_t size2, int_t n,
                  int_t* Cx) {
        static cpx_t A[LARGE * 2], B[LARGE * 2];
        for (int_t i = 0; i < (1 << size2); i++) {
            if (i < n) {
                A[i] = Ax[i];
                B[i] = Bx[i];
            } else {
                A[i] = B[i] = 0;
            }
        }
        transform(A, size2, 1);
        transform(B, size2, 1);
        for (int i = 0; i < (1 << size2); i++) A[i] *= B[i];
        transform(A, size2, -1);
        for (int_t i = 0; i < n; i++) {
            Cx[i] = (int_t((A[i]).real / (1 << size2) + 0.5) % mod);
            assert(Cx[i] >= 0);
        }
    }
    int_t process() {
        static int_t A[LARGE + 1], B[LARGE + 1];
    
        int_t n, k, l, f;
        cin >> n >> k >> l >> f;
        int_t size2 = 0;
        int_t nx = l * f + 1;
        while ((1 << size2) < (nx * 2)) size2++;
        for (int_t i = 0; i < (1 << size2); i++) {
            const auto flip = [&](int_t x) {
                int_t result = 0;
                for (int_t i = 0; i < size2 - 1; i++) {
                    result |= (x & 1);
                    x >>= 1;
                    result <<= 1;
                }
                return result | (x & 1);
            };
            flips[i] = flip(i);
        }
        std::fill(A, A + (1 << size2), 0);
        std::fill(B, B + (1 << size2), 0);
    
        if (k % mod == 0) return 0;
        for (int_t i = 1; i <= f; i++) {
            A[i] = S2[f][i];
        }
        B[0] = 1;
    #ifdef DEBUG
        cout << "size2 = " << size2 << " len = " << (1 << size2) << endl;
        cout << "nx  = " << nx << endl;
        cout << "A = ";
        for (int_t i = 0; i <= f; i++) cout << A[i] << " ";
        cout << endl;
    
    #endif
    
        if (f == 1) {
            B[l] = 1;
        } else {
            while (l) {
                assert(l >= 0);
                if (l & 1) poly_mul(A, B, size2, nx, B);
                poly_mul(A, A, size2, nx, A);
                l >>= 1;
            }
        }
    #ifdef DEBUG
        cout << "mul ok B = ";
        for (int_t i = 0; i < nx; i += 1) cout << B[i] << " ";
        cout << endl;
    
    #endif
        int_t result = 0;
        int_t kinv = inv[k % mod];
        int_t under = kinv;
        int_t prod = n;
        for (int_t i = 1; i < nx; i++) {
            result = (result + B[i] % mod * under * prod % mod) % mod;
    #ifdef DEBUG
            cout << "at " << i << " count on " << B[i] * under * prod % mod << endl;
    #endif
            under = under * kinv % mod;
            prod = prod * (n - i) % mod;
            if (n <= i) break;
        }
        return result;
    }
    int main() {
        S2[0][0] = 1;
        for (int_t i = 1; i < mod; i++) {
            for (int_t j = 1; j <= i; j++) {
                S2[i][j] = (S2[i - 1][j - 1] + S2[i - 1][j] * j) % mod;
            }
        }
        fact[0] = fact[1] = factInv[0] = factInv[1] = inv[1] = 1;
        for (int_t i = 2; i < mod; i++) {
            inv[i] = (mod - mod / i) * inv[mod % i] % mod;
            factInv[i] = factInv[i - 1] * inv[i] % mod;
            fact[i] = fact[i - 1] * i % mod;
        }
        int_t T;
        cin >> T;
        while (T--) cout << process() << endl;
    
        return 0;
    }

     

  • 洛谷4705 玩游戏

    $$ \text{原式}=\frac{1}{nm}\sum_{1\le k\le t}{x^k\sum_{1\le i\le n}{\sum_{1\le j\le m}{\left( a_i+b_j \right) ^k}}} \\ =\frac{1}{nm}\sum_{1\le k\le t}{x^k\sum_{1\le i\le n}{\sum_{1\le j\le m}{\sum_{0\le y\le k}{\left( \begin{array}{c} k\\ y\\ \end{array} \right) a_{i}^{y}b_{j}^{k-y}}}}} \\ =\frac{1}{nm}\sum_{1\le k\le t}{x^kk!\sum_{0\le y\le k}{\frac{\sum_{1\le i\le n}{a_{i}^{y}}}{y!}\times \frac{\sum_{1\le j\le m}{b_{j}^{k-y}}}{\left( k-y \right) !}}} \\ \\ \text{考虑如何计算 }\sum_{1\le i\le n}{a_{i}^{y}} \\ \text{令}F\left( x \right) =\prod_{1\le i\le n}{\left( 1+a_ix \right)},\text{显然}F\left( x \right) \text{可以在}O\left( n\log ^2n \right) \text{的时间内计算出。} \\ \ln F\left( x \right) =\sum_{1\le i\le n}{\ln \left( 1+a_ix \right)},\text{把}\ln \left( x \right) \text{在}x=\text{1处展开可得} \\ \ln F\left( x \right) =\sum_{1\le i\le n}{\sum_{j\ge 1}{\frac{\left( -1 \right) ^{j-1}}{j}\left( a_ix \right) ^j}} \\ =\sum_{j\ge 1}{x^j\frac{\left( -1 \right) ^{j-1}}{j}\sum_{1\le i\le n}{a_i^j}} \\ \text{即}i\text{次项前的系数除以}\frac{\left( -1 \right) ^{i-1}}{i}\text{即}i\text{次幂的和}. $$

    #include <cstring>
    #include <iostream>
    #include <vector>
    using int_t = long long int;
    using std::cin;
    using std::cout;
    using std::endl;
    
    const int_t mod = 998244353;
    const int_t g = 3;
    const int_t LARGE = (1 << 20);
    
    void transform(int_t* A, int_t size2, int_t arg);
    void poly_inv(const int_t* A, int_t n, int_t* result);
    std::vector<int> process(const std::vector<int>& vec);
    void poly_log(const int_t* A, int_t n, int_t* result);
    int_t power(int_t base, int_t index) {
        const auto phi = mod - 1;
        if (index < 0) index = (index % phi + phi) % phi;
        int_t result = 1;
        while (index) {
            if (index & 1) result = result * base % mod;
            index >>= 1;
            base = base * base % mod;
        }
        return result;
    }
    std::vector<int> process(const std::vector<int>& vec) {
        if (vec.size() == 1) {
            return std::vector<int>{1, vec[0]};
        }
        std::vector<int> left, right;
        for (int i = 0; i < vec.size(); i++) {
            if (i < vec.size() / 2)
                left.push_back(vec[i]);
            else
                right.push_back(vec[i]);
        }
        auto leftRet = process(std::move(left));
        auto rightRet = process(std::move(right));
        static int_t A[LARGE], B[LARGE];
        int size2 = 0;
        while ((1 << size2) < leftRet.size() + rightRet.size()) size2++;
        std::copy(leftRet.begin(), leftRet.end(), A);
        std::copy(rightRet.begin(), rightRet.end(), B);
        transform(A, size2, 1);
        transform(B, size2, 1);
        for (int i = 0; i < (1 << size2); i++) A[i] = A[i] * B[i] % mod;
        transform(A, size2, -1);
        const auto inv = power(1 << size2, -1);
        std::vector<int> result;
        for (int i = 0; i < leftRet.size() + rightRet.size() - 1; i++)
            result.push_back(A[i] * inv % mod);
        std::fill(A, A + (1 << size2), 0);
        std::fill(B, B + (1 << size2), 0);
        return result;
    }
    int main() {
        std::vector<int> polyA;
        std::vector<int> polyB;
        static int_t fact[LARGE + 1], factInv[LARGE + 1];
        int n, m, t;
        scanf("%d%d", &n, &m);
        for (int i = 1; i <= n; i++) {
            polyA.push_back(0);
            scanf("%d", &polyA.back());
        }
        for (int i = 1; i <= m; i++) {
            polyB.push_back(0);
            scanf("%d", &polyB.back());
        }
        scanf("%d", &t);
        fact[0] = factInv[0] = fact[1] = factInv[1] = 1;
        for (int_t i = 2; i <= LARGE; i++) {
            fact[i] = fact[i - 1] * i % mod;
            factInv[i] = (mod - mod / i) * factInv[mod % i] % mod;
        }
        for (int_t i = 2; i <= LARGE; i++)
            factInv[i] = factInv[i - 1] * factInv[i] % mod;
        const auto retA = process(polyA);
        const auto retB = process(polyB);
        static int_t A[LARGE], B[LARGE], Ax[LARGE], Bx[LARGE];
        std::copy(retA.begin(), retA.end(), A);
        std::copy(retB.begin(), retB.end(), B);
        poly_log(A, t + 1, Ax);
        poly_log(B, t + 1, Bx);
        memset(A, 0, sizeof(A));
        memset(B, 0, sizeof(B));
        A[0] = n;
        B[0] = m;
        for (int i = 1; i <= t; i++) {
            A[i] =
                Ax[i] * i % mod * ((i % 2) ? 1 : mod - 1) % mod * factInv[i] % mod;
            B[i] =
                Bx[i] * i % mod * ((i % 2) ? 1 : mod - 1) % mod * factInv[i] % mod;
        }
        int size2 = 0, size = 1;
        while ((1 << size2) < 2 * t) size2++;
        size = (1 << size2);
        transform(A, size2, 1);
        transform(B, size2, 1);
        for (int i = 0; i < size; i++) A[i] = A[i] * B[i] % mod;
        transform(A, size2, -1);
        const auto inv = power((int_t)size * n % mod * m % mod, -1);
        for (int i = 1; i <= t; i++) {
            printf("%lld\n", A[i] * inv % mod * fact[i] % mod);
        }
    
        return 0;
    }
    
    void poly_log(const int_t* A, int_t n, int_t* result) {
        int_t size = 1, size2 = 0;
        while ((1 << size2) < n * 2) size2++;
        size = (1 << size2);
        static int_t under[LARGE + 1], Ax[LARGE + 1];
        for (int_t i = 0; i < size; i++) {
            if (i < n)
                Ax[i] = A[i + 1] * (i + 1) % mod;
            else
                under[i] = Ax[i] = 0;
            under[i] = 0;
        }
        poly_inv(A, n, under);
        transform(under, size2, 1);
        transform(Ax, size2, 1);
        for (int i = 0; i < size; i++) Ax[i] = Ax[i] * under[i] % mod;
        transform(Ax, size2, -1);
        const int_t inv = power(size, -1);
        for (int i = 1; i < n; i++)
            result[i] = Ax[i - 1] * power(i, -1) % mod * inv % mod;
        result[0] = 0;
    }
    // C(x)=2B(x)-B(x)^2A(x)
    void poly_inv(const int_t* A, int_t n, int_t* result) {
        if (n == 1) {
            result[0] = power(A[0], -1);
            return;
        }
        poly_inv(A, n / 2 + n % 2, result);
        static int_t Ax[LARGE + 1];
        int_t size2 = 0;
        while ((1 << size2) < 3 * n) size2++;
        for (int_t i = 0; i < (1 << size2); i++) {
            if (i < n)
                Ax[i] = A[i];
            else
                Ax[i] = 0;
        }
    
        transform(Ax, size2, 1);
        transform(result, size2, 1);
        for (int_t i = 0; i < (1 << size2); i++)
            result[i] = (2 * result[i] % mod -
                         result[i] * result[i] % mod * Ax[i] % mod + mod) %
                        mod;
        transform(result, size2, -1);
        const int_t inv = power(1 << size2, -1);
        for (int_t i = 0; i < (1 << size2); i++)
            if (i < n)
                result[i] = result[i] * inv % mod;
            else
                result[i] = 0;
    }
    
    void transform(int_t* A, int_t size2, int_t arg) {
        const auto bitRev = [&](int_t x) {
            int_t result = 0;
            for (int_t i = 1; i < size2; i++) {
                result |= (x & 1);
                result <<= 1;
                x >>= 1;
            }
            return result | (x & 1);
        };
        for (int_t i = 0; i < (1 << size2); i++) {
            int_t x = bitRev(i);
            if (x > i) std::swap(A[i], A[x]);
        }
        for (int_t i = 2; i <= (1 << size2); i *= 2) {
            int_t mr = power(g, (mod - 1) / i * arg);
            for (int_t j = 0; j < (1 << size2); j += i) {
                int_t curr = 1;
                for (int_t k = 0; k < i / 2; k++) {
                    int_t u = A[j + k], t = A[j + k + i / 2] * curr % mod;
                    A[j + k] = (u + t) % mod;
                    A[j + k + i / 2] = (u - t + mod) % mod;
                    curr = curr * mr % mod;
                }
            }
        }
    }

     

  • 多项式反三角函数

    $$ \text{由欧拉公式得}e^{i\theta}=\cos \left( \theta \right) +i\sin \left( \theta \right) \\ e^{-i\theta}=\cos \left( \theta \right) -i\sin \left( \theta \right) \\ \text{故}\sin \left( \theta \right) =\frac{e^{i\theta}-e^{-i\theta}}{2i} \\ \text{设}\sin \left( \theta \right) =x,\arcsin \left( x \right) =\theta \\ \text{设}y=e^{i\theta} \\ \text{则}x=\frac{y-\frac{1}{y}}{2i} \\ 2ix=y-\frac{1}{y} \\ y^2-2ixy-1=0 \\ y=\frac{2ix\pm \sqrt{-4x^2+4}}{2}=\frac{2ix\pm 2\sqrt{1-x^2}}{2}=ix\pm \sqrt{1-x^2} \\ \text{令}\theta =\text{0,可得}x=\text{0,}y=\text{1,故}y=ix+\sqrt{1-x^2} \\ \text{所以}e^{i\theta}=ix+\sqrt{1-x^2} \\ i\theta =\log \left( ix+\sqrt{1-x^2} \right) \\ \theta =\arcsin \left( x \right) =\frac{1}{i}\log \left( ix+\sqrt{1-x^2} \right) =-i\log \left( ix+\sqrt{1-x^2} \right) \\ \arcsin \left( F\left( x \right) \right) =-i\log \left( iF\left( x \right) +\sqrt{1-F\left( x \right) ^2} \right) \\ \cos \left( \theta \right) =\frac{e^{i\theta}+e^{-i\theta}}{2} \\ \text{同上可证}\arccos \left( x \right) =-i\log \left( x\pm \sqrt{x^2-1} \right) \\ \text{同上可证}\arctan \left( x \right) =\frac{i}{2}\log \left( \frac{1-iz}{1+iz} \right) $$

    #include <assert.h>
    #include <cmath>
    #include <cstring>
    #include <iostream>
    using int_t = long long int;
    
    using std::cin;
    using std::cout;
    using std::endl;
    
    const int_t mod = 998244353;
    const int_t g = 3;
    const int_t LARGE = 1 << 19;
    template <typename T>
    T power(T base, int_t index);
    void poly_inv(const struct complex* A, int_t n, struct complex* result);
    int_t bitRev(int_t base, int_t size2);
    int_t upper2n(int_t x);
    void poly_log(const struct complex* A, int_t n, struct complex* result);
    void transform(complex* A, int_t len, int_t arg);
    void poly_exp(const struct complex* A, int_t n, struct complex* result);
    void poly_sqrt(const struct complex* A, int_t n, struct complex* result);
    struct complex {
        int_t real, imag;
        complex operator+(const complex& another) const {
            return complex{(real + another.real) % mod,
                           (imag + another.imag) % mod};
        }
        complex operator-(const complex& another) const {
            return complex{(real - another.real + mod) % mod,
                           (imag - another.imag + mod) % mod};
        }
        complex operator*(const complex& another) const {
            return complex{
                (real * another.real % mod - imag * another.imag % mod + mod) % mod,
                (real * another.imag % mod + imag * another.real % mod) % mod};
        }
        complex operator/(const complex& another) const {
            int_t a = real, b = imag, c = another.real, d = another.imag;
            int_t inv = power(c * c % mod + d * d % mod, -1);
            return complex{(a * c % mod + b * d % mod) % mod * inv % mod,
                           (b * c % mod - a * d % mod + mod) % mod * inv % mod};
        }
        complex operator%(const int_t& x) const {
            return complex{real % x, imag % x};
        }
        complex(int_t a = 0, int_t b = 0) {
            this->real = a;
            this->imag = b;
        }
        complex& operator=(const int_t& rhs) {
            this->real = rhs;
            this->imag = 0;
            return *this;
        }
    };
    std::ostream& operator<<(std::ostream& os, const complex& comp) {
        os << "complex{real=" << comp.real << ",imag=" << comp.imag << "}";
        return os;
    }
    int revs[20][LARGE + 1];
    
    int main() {
        for (int i = 0; (1 << i) <= LARGE; i++) {
            for (int j = 0; j < LARGE; j++) {
                revs[i][j] = bitRev(j, i);
            }
        }
        static complex A[LARGE], Ax[LARGE], B[LARGE];
        int_t n, type;
        scanf("%lld%lld", &n, &type);
        if (type == 0) {
            for (int_t i = 0; i < n; i++) {
                scanf("%lld", &A[i].real);
                Ax[i] = A[i];
            }
            int_t size2 = upper2n(n + n);
            transform(A, size2, 1);
            for (int_t i = 0; i < size2; i++) A[i] = A[i] * A[i];
            transform(A, size2, -1);
            const auto inv = complex{1, 0} / complex{size2, 0};
            for (int_t i = 0; i < size2; i++) {
                if (i < n)
                    A[i] = A[i] * inv;
                else
                    A[i] = 0;
            }
            A[0].real = 1;
            for (int_t i = 1; i < n; i++) A[i] = A[i] * complex{mod - 1, 0};
            poly_sqrt(A, n, B);
            for (int_t i = 0; i < n; i++) B[i] = B[i] + Ax[i] * complex{0, 1};
            memset(A, 0, sizeof(A));
            poly_log(B, n, A);
            for (int_t i = 0; i < n; i++) {
                A[i] = A[i] * complex{0, mod - 1};
            }
            for (int_t i = 0; i < n; i++) {
                printf("%lld ", A[i].real);
            }
        } else {
            for (int_t i = 0; i < n; i++) {
                int_t x;
                scanf("%lld", &x);
                A[i].imag = (mod - x) % mod;
                Ax[i].imag = x % mod;
            }
            Ax[0].real = A[0].real = 1;
            poly_inv(Ax, n, B);
            int_t size2 = upper2n(2 * n);
            transform(B, size2, 1);
            transform(A, size2, 1);
            for (int_t i = 0; i < size2; i++) B[i] = A[i] * B[i];
            transform(B, size2, -1);
            for (int_t i = 0; i < size2; i++) {
                if (i < n)
                    B[i] = B[i] / complex{size2, 0};
                else
                    B[i] = 0;
            }
            memset(A, 0, sizeof(A));
            poly_log(B, n, A);
            for (int_t i = 0; i < n; i++) {
                A[i] = A[i] / complex{2, 0} * complex{0, 1};
                printf("%lld ", A[i].real);
            }
        }
        return 0;
    }
    
    void transform(complex* A, int_t len, int_t arg) {
        int_t size2 = log2(len);
        for (int_t i = 0; i < len; i++) {
            int_t x = revs[size2][i];
            if (x > i) std::swap(A[i], A[x]);
        }
        for (int_t i = 2; i <= len; i *= 2) {
            complex mr = power(complex{g, 0}, arg * (mod - 1) / i);
            for (int_t j = 0; j < len; j += i) {
                complex curr{1, 0};
                for (int_t k = 0; k < i / 2; k++) {
                    complex u = A[j + k];
                    complex t = A[j + k + i / 2] * curr;
                    A[j + k] = (u + t);
                    A[j + k + i / 2] = (u - t);
                    curr = curr * mr;
                }
            }
        }
    }
    //计算多项式A在模x^n下的逆元
    // C(x)<-2B(x)-A(x)B^2(x)
    void poly_inv(const complex* A, int_t n, complex* result) {
        if (n == 1) {
            result[0] = power(A[0], -1);
            return;
        }
        poly_inv(A, n / 2 + n % 2, result);
        static complex temp[LARGE + 1];
        int_t size2 = upper2n(3 * n + 1);
        for (int_t i = 0; i < size2; i++) {
            if (i < n)
                temp[i] = A[i];
            else
                temp[i] = complex{0, 0};
        }
        transform(temp, size2, 1);
        transform(result, size2, 1);
        for (int_t i = 0; i < size2; i++) {
            result[i] =
                (complex{2, 0} * result[i] - temp[i] * result[i] * result[i]);
        }
        transform(result, size2, -1);
        int_t inv = power(size2, -1);
        for (int_t i = 0; i < size2; i++) {
            if (i < n) {
                result[i] = result[i] * complex{inv, 0};
            } else {
                result[i] = complex{0, 0};
            }
        }
    }
    
    void poly_log(const complex* A, int_t n, complex* result) {
        int_t size2 = upper2n(2 * n);
        static complex Ad[LARGE + 1];
        for (int_t i = 0; i < size2; i++) {
            if (i < n) {
                Ad[i] = complex{(i + 1), 0} * A[i + 1];
    
            } else {
                Ad[i] = complex{0, 0};
            }
            result[i] = complex{0, 0};
        }
        transform(Ad, size2, 1);
        poly_inv(A, n, result);
        transform(result, size2, 1);
        for (int_t i = 0; i < size2; i++) {
            result[i] = result[i] * Ad[i];
        }
        transform(result, size2, -1);
        int_t inv = power(size2, -1);
        for (int_t i = 0; i < size2; i++)
            if (i < n)
                result[i] = result[i] * complex{inv, 0};
            else
                result[i] = complex{0, 0};
        for (int_t i = n - 1; i >= 1; i--) {
            result[i] = result[i - 1] / complex{i, 0};
        }
        result[0] = complex{0, 0};
    }
    void poly_exp(const complex* A, int_t n, complex* result) {
        if (n == 1) {
            // assert(A[0] == 0);
            result[0] = complex{1, 0};
            return;
        }
        poly_exp(A, n / 2 + n % 2, result);
        static complex G0[LARGE + 1], Ax[LARGE + 1];
        int_t size2 = upper2n(2 * n);
        poly_log(result, n, G0);
        for (int_t i = 0; i < size2; i++) {
            if (i < n)
                Ax[i] = A[i];
            else
                Ax[i] = complex{0, 0};
        }
        transform(Ax, size2, 1);
        transform(G0, size2, 1);
        transform(result, size2, 1);
        for (int_t i = 0; i < size2; i++)
            result[i] = (result[i] - result[i] * (G0[i] - Ax[i]));
        transform(result, size2, -1);
        int_t inv = power(size2, -1);
        for (int_t i = 0; i < size2; i++)
            if (i < n)
                result[i] = result[i] * complex{inv, 0};
            else
                result[i] = complex{0, 0};
    }
    void poly_sqrt(const complex* A, int_t n, complex* result) {
        if (n == 1) {
            assert(A[0].real <= 1);
            result[0] = A[0];
            return;
        }
        poly_sqrt(A, n / 2 + n % 2, result);
        int_t size2 = upper2n(3 * n);
        static complex Ax[LARGE + 1], pInv[LARGE];
        for (int_t i = 0; i < size2; i++) {
            if (i < n) {
                Ax[i] = A[i];
            } else {
                Ax[i] = 0;
            }
            pInv[i] = 0;
        }
        poly_inv(result, n, pInv);
        transform(Ax, size2, 1);
        transform(result, size2, 1);
        transform(pInv, size2, 1);
        const int_t inv2 = power(2ll, -1);
        for (int_t i = 0; i < size2; i++) {
            result[i] = (result[i] -
                         (result[i] * result[i] % mod - Ax[i] + mod) % mod * inv2 %
                             mod * pInv[i] % mod +
                         mod) %
                        mod;
        }
        transform(result, size2, -1);
        for (int_t i = 0; i < size2; i++) {
            if (i < n)
                result[i] = result[i] * power(size2, -1) % mod;
            else
                result[i] = 0;
        }
    }
    template <typename T>
    T power(T base, int_t index) {
        const int_t phi = mod - 1;
        index = (index % phi + phi) % phi;
        T result = 1;
        while (index) {
            if (index & 1) result = result * base % mod;
            base = base * base % mod;
            index >>= 1;
        }
        return result;
    }
    int_t bitRev(int_t base, int_t size2) {
        int_t result = 0;
        for (int_t i = 1; i < size2; i++) {
            result |= (base & 1);
            base >>= 1;
            result <<= 1;
        }
        result |= (base & 1);
        return result;
    }
    int_t upper2n(int_t x) {
        int_t result = 1;
        while (result < x) result *= 2;
        return result;
    }

     

  • 多项式三角函数 & 多项式幂函数

    多项式三角函数

    根据欧拉公式,$e^{iF(x)}=\cos(F(x))+i\sin(F(x))$,令$F(x)$为模意义下的复多项式,直接计算即可。

    #include <cmath>
    #include <iostream>
    using int_t = long long int;
    
    using std::cin;
    using std::cout;
    using std::endl;
    
    const int_t mod = 998244353;
    const int_t g = 3;
    const int_t LARGE = 1 << 19;
    template <typename T>
    T power(T base, int_t index);
    void poly_inv(const struct complex* A, int_t n, struct complex* result);
    int_t bitRev(int_t base, int_t size2);
    int_t upper2n(int_t x);
    void poly_log(const struct complex* A, int_t n, struct complex* result);
    void transform(complex* A, int_t len, int_t arg);
    void poly_exp(const struct complex* A, int_t n, struct complex* result);
    struct complex {
        int_t real, imag;
        complex operator+(const complex& another) const {
            return complex{(real + another.real) % mod,
                           (imag + another.imag) % mod};
        }
        complex operator-(const complex& another) const {
            return complex{(real - another.real + mod) % mod,
                           (imag - another.imag + mod) % mod};
        }
        complex operator*(const complex& another) const {
            return complex{
                (real * another.real % mod - imag * another.imag % mod + mod) % mod,
                (real * another.imag % mod + imag * another.real % mod) % mod};
        }
        complex operator/(const complex& another) const {
            int_t a = real, b = imag, c = another.real, d = another.imag;
            int_t inv = power(c * c % mod + d * d % mod, -1);
            return complex{(a * c % mod + b * d % mod) % mod * inv % mod,
                           (b * c % mod - a * d % mod + mod) % mod * inv % mod};
        }
        complex operator%(const int_t& x) const {
            return complex{real % x, imag % x};
        }
        complex(int_t a = 0, int_t b = 0) {
            this->real = a;
            this->imag = b;
        }
        complex& operator=(const int_t& rhs) {
            this->real = rhs;
            this->imag = 0;
            return *this;
        }
    };
    std::ostream& operator<<(std::ostream& os, const complex& comp) {
        os << "complex{real=" << comp.real << ",imag=" << comp.imag << "}";
        return os;
    }
    int_t revs[20][LARGE + 1];
    
    int main() {
        for (int i = 0; (1 << i) <= LARGE; i++) {
            for (int j = 0; j < LARGE; j++) {
                revs[i][j] = bitRev(j, i);
            }
        }
        static complex A[LARGE], B[LARGE];
        int_t n, type;
        scanf("%lld%lld", &n, &type);
        for (int_t i = 0; i < n; i++) scanf("%lld", &A[i].imag);
        poly_exp(A, n, B);
        if (type == 0) {
            for (int_t i = 0; i < n; i++) printf("%lld ", B[i].imag);
        } else {
            for (int_t i = 0; i < n; i++) printf("%lld ", B[i].real);
        }
        return 0;
    }
    
    void transform(complex* A, int_t len, int_t arg) {
        int_t size2 = log2(len);
        for (int_t i = 0; i < len; i++) {
            int_t x = revs[size2][i];
            if (x > i) std::swap(A[i], A[x]);
        }
        for (int_t i = 2; i <= len; i *= 2) {
            complex mr = power(complex{g, 0}, arg * (mod - 1) / i);
            for (int_t j = 0; j < len; j += i) {
                complex curr{1, 0};
                for (int_t k = 0; k < i / 2; k++) {
                    complex u = A[j + k];
                    complex t = A[j + k + i / 2] * curr;
                    A[j + k] = (u + t);
                    A[j + k + i / 2] = (u - t);
                    curr = curr * mr;
                }
            }
        }
    }
    //计算多项式A在模x^n下的逆元
    // C(x)<-2B(x)-A(x)B^2(x)
    void poly_inv(const complex* A, int_t n, complex* result) {
        if (n == 1) {
            result[0] = power(A[0], -1);
            return;
        }
        poly_inv(A, n / 2 + n % 2, result);
        static complex temp[LARGE + 1];
        int_t size2 = upper2n(3 * n + 1);
        for (int_t i = 0; i < size2; i++) {
            if (i < n)
                temp[i] = A[i];
            else
                temp[i] = complex{0, 0};
        }
        transform(temp, size2, 1);
        transform(result, size2, 1);
        for (int_t i = 0; i < size2; i++) {
            result[i] =
                (complex{2, 0} * result[i] - temp[i] * result[i] * result[i]);
        }
        transform(result, size2, -1);
        int_t inv = power(size2, -1);
        for (int_t i = 0; i < size2; i++) {
            if (i < n) {
                result[i] = result[i] * complex{inv, 0};
            } else {
                result[i] = complex{0, 0};
            }
        }
    }
    
    void poly_log(const complex* A, int_t n, complex* result) {
        int_t size2 = upper2n(2 * n);
        static complex Ad[LARGE + 1];
        for (int_t i = 0; i < size2; i++) {
            if (i < n) {
                Ad[i] = complex{(i + 1), 0} * A[i + 1];
    
            } else {
                Ad[i] = complex{0, 0};
            }
            result[i] = complex{0, 0};
        }
        transform(Ad, size2, 1);
        poly_inv(A, n, result);
        transform(result, size2, 1);
        for (int_t i = 0; i < size2; i++) {
            result[i] = result[i] * Ad[i];
        }
        transform(result, size2, -1);
        int_t inv = power(size2, -1);
        for (int_t i = 0; i < size2; i++)
            if (i < n)
                result[i] = result[i] * complex{inv, 0};
            else
                result[i] = complex{0, 0};
        for (int_t i = n - 1; i >= 1; i--) {
            result[i] = result[i - 1] / complex{i, 0};
        }
        result[0] = complex{0, 0};
    }
    void poly_exp(const complex* A, int_t n, complex* result) {
        if (n == 1) {
            // assert(A[0] == 0);
            result[0] = complex{1, 0};
            return;
        }
        poly_exp(A, n / 2 + n % 2, result);
        static complex G0[LARGE + 1], Ax[LARGE + 1];
        int_t size2 = upper2n(2 * n);
        poly_log(result, n, G0);
        for (int_t i = 0; i < size2; i++) {
            if (i < n)
                Ax[i] = A[i];
            else
                Ax[i] = complex{0, 0};
        }
        transform(Ax, size2, 1);
        transform(G0, size2, 1);
        transform(result, size2, 1);
        for (int_t i = 0; i < size2; i++)
            result[i] = (result[i] - result[i] * (G0[i] - Ax[i]));
        transform(result, size2, -1);
        int_t inv = power(size2, -1);
        for (int_t i = 0; i < size2; i++)
            if (i < n)
                result[i] = result[i] * complex{inv, 0};
            else
                result[i] = complex{0, 0};
    }
    template <typename T>
    T power(T base, int_t index) {
        const int_t phi = mod - 1;
        index = (index % phi + phi) % phi;
        T result = 1;
        while (index) {
            if (index & 1) result = result * base % mod;
            base = base * base % mod;
            index >>= 1;
        }
        return result;
    }
    int_t bitRev(int_t base, int_t size2) {
        int_t result = 0;
        for (int_t i = 1; i < size2; i++) {
            result |= (base & 1);
            base >>= 1;
            result <<= 1;
        }
        result |= (base & 1);
        return result;
    }
    int_t upper2n(int_t x) {
        int_t result = 1;
        while (result < x) result *= 2;
        return result;
    }

    多项式幂函数

    求$F(x)^m$。

    令$F(x)=cx^pG(x)$,其中$G(x)$为常数项为1多项式。

    则$F(x)^m=(cx^p)^mG(x)^m=c^m\times x^{pm}\times e^{m\log G(x)}$。

    // luogu-judger-enable-o2
    #include <assert.h>
    #include <algorithm>
    #include <cmath>
    #include <cstring>
    #include <ctime>
    #include <iostream>
    using int_t = long long int;
    using std::cin;
    using std::cout;
    using std::endl;
    const int mod = 998244353;
    const int g = 3;
    const int LARGE = 1 << 19;
    int revs[20][LARGE + 1];
    
    int power(int base, int index);
    void transform(int* A, int len, int arg);
    void poly_inv(const int* A, int n, int* result);
    void poly_sqrt(const int* A, int n, int* result);
    void poly_log(const int* A, int n, int* result);
    void poly_exp(const int* A, int n, int* result);
    void poly_power(const int* A, int n, int index, int* result);
    int modSqrt(int a);
    int bitRev(int base, int size2) {
        int result = 0;
        for (int i = 1; i < size2; i++) {
            result |= (base & 1);
            base >>= 1;
            result <<= 1;
        }
        result |= (base & 1);
        return result;
    }
    int upper2n(int x) {
        int result = 1;
        while (result < x) result *= 2;
        return result;
    }
    int_t readInt() {
        int_t result = 0;
        char chr = getchar();
        while (isdigit(chr) == false) chr = getchar();
        while (isdigit(chr)) {
            result = (result * 10 + chr - '0') % mod;
            chr = getchar();
        }
        return result;
    }
    int main() {
        for (int i = 0; (1 << i) <= LARGE; i++) {
            for (int j = 0; j < LARGE; j++) {
                revs[i][j] = bitRev(j, i);
            }
        }
        static int A[LARGE + 1], B[LARGE + 1];
        int_t n, k;
        scanf("%lld%lld", &n, &k);
        for (int i = 0; i < n; i++) scanf("%d", &A[i]);
    
        int_t x0pos = -1;
        for (int_t i = 0; i < n; i++) {
            if (A[i]) {
                x0pos = i;
                break;
            }
        }
        if (x0pos == -1) {
            for (int_t i = 0; i < n; i++) printf("0 ");
            return 0;
        }
        // x0pos += 1;
        for (int_t i = 0; i < n; i++) A[i] = A[i + x0pos];
        int_t x0 = A[0];
        const int_t inv = power(x0, -1);
        for (int_t i = 0; i < n; i++) A[i] = A[i] * inv % mod;
        poly_log(A, n, B);
        for (int_t i = 0; i < n; i++) B[i] = B[i] * k % mod;
        memset(A, 0, sizeof(A));
        poly_exp(B, n, A);
        for (int_t i = 0; i < n; i++) A[i] = (int_t)A[i] * power(x0, k) % mod;
        int_t time = x0pos * k;
        if (time <= n - 1) {
            for (int_t i = n - 1; i >= time; i--) A[i] = A[i - time];
            for (int_t i = time - 1; i >= 0; i--) A[i] = 0;
            for (int_t i = 0; i < n; i++) printf("%d ", A[i]);
        } else {
            for (int_t i = 0; i < n; i++) printf("0 ");
        }
        return 0;
    }
    //计算A(x)^index mod x^n
    void poly_power(const int* A, int n, int index, int* result) {
        int size2 = upper2n(2 * n);
        static int base[LARGE + 1];
        for (int i = 0; i < size2; i++) {
            if (i < n)
                base[i] = A[i];
            else
                base[i] = 0;
            result[i] = 0;
        }
        result[0] = 1;
        while (index) {
            transform(base, size2, 1);
            if (index & 1) {
                transform(result, size2, 1);
                for (int i = 0; i < size2; i++) {
                    result[i] = (int_t)result[i] * base[i] % mod;
                }
                transform(result, size2, -1);
                for (int i = 0; i < size2; i++) {
                    if (i < n)
                        result[i] = (int_t)result[i] * power(size2, -1) % mod;
                    else
                        result[i] = 0;
                }
            }
            for (int i = 0; i < size2; i++)
                base[i] = (int_t)base[i] * base[i] % mod;
            transform(base, size2, -1);
            for (int i = 0; i < size2; i++) {
                if (i < n)
                    base[i] = (int_t)base[i] * power(size2, -1) % mod;
                else
                    base[i] = 0;
            }
            index >>= 1;
        }
    }
    int power(int base, int index) {
        const int phi = mod - 1;
        index = (index % phi + phi) % phi;
        base = (base % mod + mod) % mod;
        int result = 1;
        while (index) {
            if (index & 1) result = (int_t)result * base % mod;
            base = (int_t)base * base % mod;
            index >>= 1;
        }
        return result;
    }
    void transform(int* A, int len, int arg) {
        int size2 = log2(len);
        for (int i = 0; i < len; i++) {
            int x = revs[size2][i];
            if (x > i) std::swap(A[i], A[x]);
        }
        for (int i = 2; i <= len; i *= 2) {
            int mr = power(g, (int_t)arg * (mod - 1) / i);
            for (int j = 0; j < len; j += i) {
                int curr = 1;
                for (int k = 0; k < i / 2; k++) {
                    int u = A[j + k];
                    int t = (int_t)A[j + k + i / 2] * curr % mod;
                    A[j + k] = (u + t) % mod;
                    A[j + k + i / 2] = (u - t + mod) % mod;
                    curr = (int_t)curr * mr % mod;
                }
            }
        }
    }
    //计算多项式A在模x^n下的逆元
    // C(x)<-2B(x)-A(x)B^2(x)
    void poly_inv(const int* A, int n, int* result) {
        if (n == 1) {
            result[0] = power(A[0], -1);
            return;
        }
        poly_inv(A, n / 2 + n % 2, result);
        static int temp[LARGE + 1];
        int size2 = upper2n(3 * n + 1);
        for (int i = 0; i < size2; i++) {
            if (i < n)
                temp[i] = A[i];
            else
                temp[i] = 0;
        }
        transform(temp, size2, 1);
        transform(result, size2, 1);
        for (int i = 0; i < size2; i++) {
            result[i] =
                ((int_t)2 * result[i] % mod -
                 (int_t)temp[i] * result[i] % mod * result[i] % mod + 2 * mod) %
                mod;
        }
        transform(result, size2, -1);
        for (int i = 0; i < size2; i++) {
            if (i < n) {
                result[i] = (int_t)result[i] * power(size2, -1) % mod;
            } else {
                result[i] = 0;
            }
        }
    }
    int modSqrt(int a) {
        const int_t b = 3;
        // mod-1=2^s*t
        int s = 23, t = 119;
        int x = power(a, (t + 1) / 2);
        int e = power(a, t);
        int k = 1;
        while (k < s) {
            if (power(e, 1 << (s - k - 1)) != 1) {
                x = (int_t)x * power(b, (1 << (k - 1)) * t) % mod;
            }
            e = (int_t)power(a, -1) * x % mod * x % mod;
            k++;
        }
        return x;
    }
    //计算多项式开根
    void poly_sqrt(const int* A, int n, int* result) {
        if (n == 1) {
            int p = modSqrt(A[0]);
            result[0] = std::min(p, mod - p);
            return;
        }
        poly_sqrt(A, n / 2 + n % 2, result);
        int size2 = upper2n(3 * n);
        static int Ax[LARGE + 1], pInv[LARGE];
        for (int i = 0; i < size2; i++) {
            if (i < n) {
                Ax[i] = A[i];
            } else {
                Ax[i] = 0;
            }
            pInv[i] = 0;
        }
        poly_inv(result, n, pInv);
        transform(Ax, size2, 1);
        transform(result, size2, 1);
        transform(pInv, size2, 1);
        const int inv2 = power(2, -1);
        for (int i = 0; i < size2; i++) {
            result[i] = (result[i] -
                         ((int_t)result[i] * result[i] % mod - Ax[i] + mod) % mod *
                             inv2 % mod * pInv[i] % mod +
                         mod) %
                        mod;
        }
        transform(result, size2, -1);
        for (int i = 0; i < size2; i++) {
            if (i < n)
                result[i] = (int_t)result[i] * power(size2, -1) % mod;
            else
                result[i] = 0;
        }
    }
    void poly_log(const int* A, int n, int* result) {
        int size2 = upper2n(2 * n);
        static int Ad[LARGE + 1];
        for (int i = 0; i < size2; i++) {
            if (i < n) {
                Ad[i] = (int_t)(i + 1) * A[i + 1] % mod;
    
            } else {
                Ad[i] = 0;
            }
            result[i] = 0;
        }
        transform(Ad, size2, 1);
        poly_inv(A, n, result);
        transform(result, size2, 1);
        for (int i = 0; i < size2; i++) {
            result[i] = (int_t)result[i] * Ad[i] % mod;
        }
        transform(result, size2, -1);
        for (int i = 0; i < size2; i++)
            if (i < n)
                result[i] = (int_t)result[i] * power(size2, -1) % mod;
            else
                result[i] = 0;
        for (int i = n - 1; i >= 1; i--) {
            result[i] = (int_t)result[i - 1] * power(i, -1) % mod;
        }
        result[0] = 0;
    }
    void poly_exp(const int* A, int n, int* result) {
        if (n == 1) {
            assert(A[0] == 0);
            result[0] = 1;
            return;
        }
        poly_exp(A, n / 2 + n % 2, result);
        static int G0[LARGE + 1], Ax[LARGE + 1];
        int size2 = upper2n(2 * n);
        poly_log(result, n, G0);
        for (int_t i = 0; i < size2; i++) {
            if (i < n)
                Ax[i] = A[i];
            else
                Ax[i] = 0;
        }
        transform(Ax, size2, 1);
        transform(G0, size2, 1);
        transform(result, size2, 1);
        for (int i = 0; i < size2; i++)
            result[i] =
                (result[i] - (int_t)result[i] * (G0[i] - Ax[i] + mod) % mod + mod) %
                mod;
        transform(result, size2, -1);
        for (int i = 0; i < size2; i++)
            if (i < n)
                result[i] = (int_t)result[i] * power(size2, -1) % mod;
            else
                result[i] = 0;
    }

     

  • 洛谷4841 城市规划

    带标号无向连通图计数,使用多项式ln。

    #include <iostream>
    #include <algorithm>
    #include <cmath>
    #include <assert.h>
    #include <sstream>
    #include <cstdio>
    typedef long long int int_t;
    using std::cin;
    using std::cout;
    using std::endl;
    //1004535809
    //const int_t mod = 998244353;
    const int_t mod = 1004535809;
    const int_t g = 3;
    int_t power(int_t base, int_t index);
    int_t upper2n(int_t x);
    int_t bitReverse(int_t bits, int_t size2);
    template <int_t arg>
    void transform(int_t *A, int_t size);
    void poly_inv(int_t *A, int_t *inv, int_t n);
    void poly_log(int_t *A, int_t n);
    const int_t LARGE = 1 << 22;
    int_t C2(int_t x){
        return x * (x - 1) % mod * power(2,-1) % mod;
    }
    //4 38
     
    int main()
    {
    //  freopen("count.in","r",stdin);
    //  freopen("count.out","w",stdout);
        static int_t A[LARGE];
        static int_t fact[LARGE + 1];
        fact[0] = 1;
        int_t n;
        cin >> n;
        for(int_t i = 1;i <= n ;i++){
            fact[i] = fact[i - 1] * i % mod;
        }
        for(int_t i = 0;i <= n ;i++){
            A[i] = power(2 , i*(i-1)/2) * power(fact[i] , -1)%mod;
        }
        poly_log(A,n+1);
        for(int_t i = 0 ;i <= n  ;i ++){
            A[i] = A[i] * fact[i] % mod;
        }
        cout << A[n] << endl;
        return 0;
    }
    void poly_log(int_t *A, int_t n)
    {
        int_t size = upper2n(2 * n);
        static int_t inv[LARGE];
        std::fill(inv, inv + size, 0);
        poly_inv(A, inv, n);
        for (int_t i = 0; i < n; i++)
        {
            A[i] = (A[i + 1] * (i + 1)) % mod;
        }
        transform<1>(A, size);
        transform<1>(inv, size);
        for (int_t i = 0; i < size; i++)
        {
            A[i] = (A[i] * inv[i]) % mod;
        }
        transform<-1>(A, size);
        for (int_t i = n - 1; i >= 1; i--)
        {
            A[i] = A[i - 1] * power(i, -1) % mod * power(size, -1) % mod;
        }
        A[0] = 0;
    }
    void poly_inv(int_t *A, int_t *inv, int_t n)
    {
        static int_t Ax[LARGE];
        if (n == 1)
        {
            inv[0] = power(A[0], -1);
            return;
        }
        poly_inv(A, inv, n / 2 + n % 2);
        int_t size = upper2n(n * 3 - 1);
        for (int_t i = 0; i < size; i++)
        {
     
            if (i < n)
            {
                Ax[i] = A[i];
            }
            else
            {
                Ax[i] = 0;
            }
        }
        transform<1>(inv, size);
        transform<1>(Ax, size);
        for (int_t i = 0; i < size; i++)
        {
            inv[i] = (2 * inv[i] - inv[i] * inv[i] % mod * Ax[i] % mod + mod) % mod;
        }
        transform<-1>(inv, size);
        for (int_t i = 0; i < n; i++)
        {
            inv[i] = (inv[i] * power(size, -1)) % mod;
        }
        for (int_t i = n; i < size; i++)
        {
            inv[i] = 0;
        }
    }
    template <int_t arg >
    void transform(int_t *A, int_t size)
    {
        int_t size2 = log2(size);
        for (int_t i = 0; i < size; i++)
        {
            int_t x = bitReverse(i, size2);
            if (x > i)
            {
                std::swap(A[i], A[x]);
            }
        }
        for (int_t i = 2; i <= size; i *= 2)
        {
            int_t mr = power(g, arg * (mod - 1) / i);
            for (int_t j = 0; j < size; j += i)
            {
                int_t curr = 1;
                for (int_t k = 0; k < i / 2; k++)
                {
                    int_t u = A[j + k];
                    int_t t = A[j + k + i / 2] * curr % mod;
                    A[j + k] = (u + t) % mod;
                    A[j + k + i / 2] = (u - t + mod) % mod;
                    curr = (curr * mr) % mod;
                }
            }
        }
    }
     
    int_t bitReverse(int_t bits, int_t size2)
    {
        int_t result = 0;
        for (int_t i = 1; i < size2; i++)
        {
            result |= (bits & 1);
            result <<= 1;
            bits >>= 1;
        }
        return result | (bits & 1);
    }
    int_t upper2n(int_t x)
    {
        int_t result = 1;
        while (result < x)
            result *= 2;
        return result;
    }
    int_t power(int_t base, int_t index)
    {
        int_t result = 1;
        if (index < 0)
        {
            index *= -1;
            base = power(base, mod - 2);
        }
        while (index)
        {
            if (index & 1)
            {
                result = (result * base) % mod;
            }
            index >>= 1;
            base = (base * base) % mod;
        }
        return result;
    }
    

     

  • LOJ150 YT2sOJ21 挑战多项式

    毕竟是板子集合,不管常数了。

    #include <assert.h>
    #include <algorithm>
    #include <cmath>
    #include <cstring>
    #include <ctime>
    #include <iostream>
    using int_t = long long int;
    using std::cin;
    using std::cout;
    using std::endl;
    const int mod = 998244353;
    const int g = 3;
    const int LARGE = 1 << 19;
    int revs[20][LARGE + 1];
    
    int power(int base, int index);
    void transform(int* A, int len, int arg);
    void poly_inv(const int* A, int n, int* result);
    void poly_sqrt(const int* A, int n, int* result);
    void poly_log(const int* A, int n, int* result);
    void poly_exp(const int* A, int n, int* result);
    void poly_power(const int* A, int n, int index, int* result);
    int modSqrt(int a);
    int bitRev(int base, int size2) {
        int result = 0;
        for (int i = 1; i < size2; i++) {
            result |= (base & 1);
            base >>= 1;
            result <<= 1;
        }
        result |= (base & 1);
        return result;
    }
    int upper2n(int x) {
        int result = 1;
        while (result < x) result *= 2;
        return result;
    }
    int main() {
    #ifdef TIME
        auto begin = clock();
    #endif
        for (int i = 0; (1 << i) <= LARGE; i++) {
            for (int j = 0; j < LARGE; j++) {
                revs[i][j] = bitRev(j, i);
            }
        }
        static int F[LARGE + 1], A[LARGE + 1], B[LARGE + 1];
        int n, k;
    
        scanf("%d%d", &n, &k);
        for (int i = 0; i <= n; i++) {
            scanf("%d", &A[i]);
            F[i] = A[i];
        }
        poly_sqrt(A, n + 1, B);
    
    #ifdef DEBUG
        cout << "sqrt ";
        for (int_t i = 0; i <= n; i++) cout << B[i] << " ";
        cout << endl;
    #endif
        memset(A, 0, sizeof(A));
        poly_inv(B, n + 1, A);
    #ifdef DEBUG
        cout << "sqrt inv ";
        for (int_t i = 0; i <= n; i++) cout << A[i] << " ";
        cout << endl;
    #endif
        for (int i = n; i >= 1; i--) B[i] = (int_t)A[i - 1] * power(i, -1) % mod;
        B[0] = 0;
    #ifdef DEBUG
        cout << "sqrt inv integrate ";
        for (int_t i = 0; i <= n; i++) cout << B[i] << " ";
        cout << endl;
    #endif
        memset(A, 0, sizeof(A));
        poly_exp(B, n + 1, A);
    #ifdef DEBUG
        cout << "sqrt inv integrate exp ";
        for (int_t i = 0; i <= n; i++) cout << A[i] << " ";
        cout << endl;
    #endif
        for (int i = 0; i < n + 1; i++) {
            A[i] = (F[i] - A[i] + mod) % mod;
        }
        A[0] = (A[0] + 2 - F[0] + mod) % mod;
    #ifdef DEBUG
        cout << "sqrt inv integrate exp process ";
        for (int_t i = 0; i <= n; i++) cout << A[i] << " ";
        cout << endl;
    #endif
        // memset(B, 0, sizeof(B));
        poly_log(A, n + 1, B);
    
        B[0] = 1;
    #ifdef DEBUG
        cout << "sqrt inv integrate exp process log +1 ";
        for (int_t i = 0; i <= n; i++) cout << B[i] << " ";
        cout << endl;
    #endif
        // memset(A, 0, sizeof());
        poly_power(B, n + 1, k, A);
    #ifdef DEBUG
        cout << "sqrt inv integrate exp process power ";
        for (int_t i = 0; i <= n; i++) cout << A[i] << " ";
        cout << endl;
    #endif
        for (int i = 0; i < n; i++) {
            A[i] = (int_t)A[i + 1] * (i + 1) % mod;
            printf("%d ", A[i]);
        }
    #ifdef TIME
        auto end = clock();
        cout << 1.0 * (end - begin) / CLOCKS_PER_SEC << endl;
    #endif
        return 0;
    }
    //计算A(x)^index mod x^n
    void poly_power(const int* A, int n, int index, int* result) {
        int size2 = upper2n(2 * n);
        static int base[LARGE + 1];
        for (int i = 0; i < size2; i++) {
            if (i < n)
                base[i] = A[i];
            else
                base[i] = 0;
            result[i] = 0;
        }
        result[0] = 1;
        while (index) {
            transform(base, size2, 1);
            if (index & 1) {
                transform(result, size2, 1);
                for (int i = 0; i < size2; i++) {
                    result[i] = (int_t)result[i] * base[i] % mod;
                }
                transform(result, size2, -1);
                for (int i = 0; i < size2; i++) {
                    if (i < n)
                        result[i] = (int_t)result[i] * power(size2, -1) % mod;
                    else
                        result[i] = 0;
                }
            }
            for (int i = 0; i < size2; i++)
                base[i] = (int_t)base[i] * base[i] % mod;
            transform(base, size2, -1);
            for (int i = 0; i < size2; i++) {
                if (i < n)
                    base[i] = (int_t)base[i] * power(size2, -1) % mod;
                else
                    base[i] = 0;
            }
            index >>= 1;
        }
    }
    int power(int base, int index) {
        const int phi = mod - 1;
        index = (index % phi + phi) % phi;
        base = (base % mod + mod) % mod;
        int result = 1;
        while (index) {
            if (index & 1) result = (int_t)result * base % mod;
            base = (int_t)base * base % mod;
            index >>= 1;
        }
        return result;
    }
    void transform(int* A, int len, int arg) {
        int size2 = log2(len);
        for (int i = 0; i < len; i++) {
            int x = revs[size2][i];
            if (x > i) std::swap(A[i], A[x]);
        }
        for (int i = 2; i <= len; i *= 2) {
            int mr = power(g, (int_t)arg * (mod - 1) / i);
            for (int j = 0; j < len; j += i) {
                int curr = 1;
                for (int k = 0; k < i / 2; k++) {
                    int u = A[j + k];
                    int t = (int_t)A[j + k + i / 2] * curr % mod;
                    A[j + k] = (u + t) % mod;
                    A[j + k + i / 2] = (u - t + mod) % mod;
                    curr = (int_t)curr * mr % mod;
                }
            }
        }
    }
    //计算多项式A在模x^n下的逆元
    // C(x)<-2B(x)-A(x)B^2(x)
    void poly_inv(const int* A, int n, int* result) {
        if (n == 1) {
            result[0] = power(A[0], -1);
            return;
        }
        poly_inv(A, n / 2 + n % 2, result);
        static int temp[LARGE + 1];
        int size2 = upper2n(3 * n + 1);
        for (int i = 0; i < size2; i++) {
            if (i < n)
                temp[i] = A[i];
            else
                temp[i] = 0;
        }
        transform(temp, size2, 1);
        transform(result, size2, 1);
        for (int i = 0; i < size2; i++) {
            result[i] =
                ((int_t)2 * result[i] % mod -
                 (int_t)temp[i] * result[i] % mod * result[i] % mod + 2 * mod) %
                mod;
        }
        transform(result, size2, -1);
        for (int i = 0; i < size2; i++) {
            if (i < n) {
                result[i] = (int_t)result[i] * power(size2, -1) % mod;
            } else {
                result[i] = 0;
            }
        }
    }
    int modSqrt(int a) {
        const int_t b = 3;
        // mod-1=2^s*t
        int s = 23, t = 119;
        int x = power(a, (t + 1) / 2);
        int e = power(a, t);
        int k = 1;
        while (k < s) {
            if (power(e, 1 << (s - k - 1)) != 1) {
                x = (int_t)x * power(b, (1 << (k - 1)) * t) % mod;
            }
            e = (int_t)power(a, -1) * x % mod * x % mod;
            k++;
        }
        return x;
    }
    //计算多项式开根
    void poly_sqrt(const int* A, int n, int* result) {
        if (n == 1) {
            int p = modSqrt(A[0]);
            result[0] = std::min(p, mod - p);
            return;
        }
        poly_sqrt(A, n / 2 + n % 2, result);
        int size2 = upper2n(3 * n);
        static int Ax[LARGE + 1], pInv[LARGE];
        for (int i = 0; i < size2; i++) {
            if (i < n) {
                Ax[i] = A[i];
            } else {
                Ax[i] = 0;
            }
            pInv[i] = 0;
        }
        poly_inv(result, n, pInv);
        transform(Ax, size2, 1);
        transform(result, size2, 1);
        transform(pInv, size2, 1);
        const int inv2 = power(2, -1);
        for (int i = 0; i < size2; i++) {
            result[i] = (result[i] -
                         ((int_t)result[i] * result[i] % mod - Ax[i] + mod) % mod *
                             inv2 % mod * pInv[i] % mod +
                         mod) %
                        mod;
        }
        transform(result, size2, -1);
        for (int i = 0; i < size2; i++) {
            if (i < n)
                result[i] = (int_t)result[i] * power(size2, -1) % mod;
            else
                result[i] = 0;
        }
    }
    void poly_log(const int* A, int n, int* result) {
        int size2 = upper2n(2 * n);
        static int Ad[LARGE + 1];
        for (int i = 0; i < size2; i++) {
            if (i < n) {
                Ad[i] = (int_t)(i + 1) * A[i + 1] % mod;
    
            } else {
                Ad[i] = 0;
            }
            result[i] = 0;
        }
        transform(Ad, size2, 1);
        poly_inv(A, n, result);
        transform(result, size2, 1);
        for (int i = 0; i < size2; i++) {
            result[i] = (int_t)result[i] * Ad[i] % mod;
        }
        transform(result, size2, -1);
        for (int i = 0; i < size2; i++)
            if (i < n)
                result[i] = (int_t)result[i] * power(size2, -1) % mod;
            else
                result[i] = 0;
        for (int i = n - 1; i >= 1; i--) {
            result[i] = (int_t)result[i - 1] * power(i, -1) % mod;
        }
        result[0] = 0;
    }
    void poly_exp(const int* A, int n, int* result) {
        if (n == 1) {
            assert(A[0] == 0);
            result[0] = 1;
            return;
        }
        poly_exp(A, n / 2 + n % 2, result);
        static int G0[LARGE + 1], Ax[LARGE + 1];
        int size2 = upper2n(2 * n);
        poly_log(result, n, G0);
        for (int_t i = 0; i < size2; i++) {
            if (i < n)
                Ax[i] = A[i];
            else
                Ax[i] = 0;
        }
        transform(Ax, size2, 1);
        transform(G0, size2, 1);
        transform(result, size2, 1);
        for (int i = 0; i < size2; i++)
            result[i] =
                (result[i] - (int_t)result[i] * (G0[i] - Ax[i] + mod) % mod + mod) %
                mod;
        transform(result, size2, -1);
        for (int i = 0; i < size2; i++)
            if (i < n)
                result[i] = (int_t)result[i] * power(size2, -1) % mod;
            else
                result[i] = 0;
    }