跑的比较快的多项式板子

#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;
}

 

评论

发表回复

您的邮箱地址不会被公开。 必填项已用 * 标注

这个站点使用 Akismet 来减少垃圾评论。了解你的评论数据如何被处理