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

多项式三角函数

根据欧拉公式，$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;
}