hdu4656Evaluation 生成函数 FFT jk项的处理

hdu4656Evaluation

题目传送门

分析

沉迷生成函数，无法自拔~
先瞎推一波式子。
$$\sum _{i=0}^{n-1}{a}_i(b{c}^{2k}+d{)}^i=\sum _{i=0}^{n-1}\sum _{j=0}^i{a}_i{C}_i^j{b}^j{c}^{2kj}{d}^{i-j}$$
$$=\sum _{j=0}^{n-1}\frac{b^j{c}^{2kj}}{j!}\sum _{i=j}^{n-1}{a}_{i}i!\frac{d^{i-j}}{(i-j)!}$$
令
$$H(j)=\frac{b^j}{j!}\sum _{i=j}^{n-1}{a}_{i}i!\frac{d^{i-j}}{(i-j)!}$$
这个瞎卷一波。
然后就不会了>_<
要求
$$Ans(k)=\sum _{j=0}^{n-1}{c}^{2kj}H(j)$$
一个神奇的操作是
$$-2kj=(k-j{)}^2-k^2-j^2$$
$$Ans(k)=\sum _{j=0}^{n-1}{c}^{k^2+j^2-(k-j{)}^2}H(j)$$
$$=c^{k^2}\sum _{j=0}^{n-1}{c}^{j^2}H(j)c^{-(k-j{)}^2}$$
仍然瞎卷一波就可以了。
注意一下$k-j$可能是负数，所以要倍长一下下。。

代码

$MTT$被卡精度了。。。调了半天

#include<bits/stdc++.h>
const int N = 524288, P = 1e6 + 3;
const double pi = acos(-1.0);
int ri() {
    char c = getchar(); int x = 0, f = 1; for(;c < '0' || c > '9'; c = getchar()) if(c == '-') f = -1;
    for(;c >= '0' && c <= '9'; c = getchar()) x = (x << 1) + (x << 3) - '0' + c; return x * f;
}
int mul(int a, int b) {return 1LL * a * b % P;}
int add(int a, int b) {return (a + b) % P;}
int fix(int x) {return (x >> 31 & P) + x;}
int Pow(int x, long long k) {
    int r = 1; k %= (P - 1);
    if(k < 0) k += P - 1;
    for(;k; x = mul(x, x), k >>= 1)
        if(k & 1)
            r = mul(r, x);
    return r;
}
struct cp {
    double r, i;
    cp(double _r = 0, double _i = 0) : r(_r), i(_i) {}
    cp operator + (const cp &a) {return cp(r + a.r, i + a.i);}
    cp operator - (const cp &a) {return cp(r - a.r, i - a.i);}
    cp operator * (const cp &a) {return cp(r * a.r - i * a.i, i *a.r + a.i * r);}
}Da[N], Db[N], Dc[N], Dd[N], w[N];
int L, n, b, c, d, A[N], B[N], R[N], iv[N], fac[N], ivf[N];
void Pre(int m) {
    L = 1; int x = 0;
    for(;(L <<= 1) < m; ++x) ;
    for(int i = 1;i < L; ++i)
        R[i] = R[i >> 1] >> 1 | (i & 1) << x;
    for(int i = 0;i < L; ++i)     
        w[i] = cp(cos(2 * pi * i / L), sin(2 * pi * i / L));
}
void DFT(cp *F) {
    for(int i = 1;i < L; ++i)
        if(i < R[i])
            std::swap(F[i], F[R[i]]);
    for(int i = 1, d = L >> 1;i < L; i <<= 1, d >>= 1)
        for(int j = 0;j < L; j += i << 1) {
            cp *l = F + j, *r = F + j + i, *p = w, tp;
            for(int k = i; k--; ++l, ++r, p += d)
                tp = *r * *p, *r = *l - tp, *l = *l + tp;
        }
}
void Mul(int *a, int *b, int *c, int n, int m) {
    Pre(n + m);
    for(int i = 0;i < n; ++i) {
        Da[i] = cp(a[i] >> 15, 0);
        Db[i] = cp(a[i] & 32767, 0);
    }
    for(int i = n;i < L; ++i) 
        Da[i] = Db[i] = cp(0, 0);
    for(int i = 0;i < m; ++i) {
        Dc[i] = cp(b[i] >> 15, 0);
        Dd[i] = cp(b[i] & 32767, 0);
    }
    for(int i = m;i < L; ++i) 
        Dc[i] = Dd[i] = cp(0, 0);
    DFT(Da); DFT(Db); DFT(Dc); DFT(Dd);
    for(int i = 0;i < L; ++i) {
        cp a = Da[i], b = Db[i], c = Dc[i], d = Dd[i];
        Da[i] = a * c;
        Db[i] = b * c + a * d;
        Dc[i] = b * d;
    }
    DFT(Da); DFT(Db); DFT(Dc);
    for(int i = 0;i < n + m; ++i) {
        int j = L - i & L - 1;
        c[i] = (long long)(Da[j].r / L + 0.5) % P;
        c[i] = ((long long)c[i] << 15) % P;
        c[i] = add(c[i], (long long)(Db[j].r / L + 0.5) % P);
        c[i] = ((long long)c[i] << 15) % P;
        c[i] = add(c[i], (long long)(Dc[j].r / L + 0.5) % P);
    }
}
int main() {
    n = ri(); b = ri(); c = ri(); d = ri();
    fac[0] = ivf[0] = iv[1] = 1;
    for(int i = 1;i < n; ++i) {
        if(i != 1) iv[i] = mul(fix(-iv[P % i]), P / i);
        fac[i] = mul(fac[i - 1], i);
        ivf[i] = mul(ivf[i - 1], iv[i]);
    }
    for(int i = 0;i < n; ++i) {
        A[i] = mul(ri(), fac[i]);
        B[i] = mul(Pow(d, n - i - 1), ivf[n - i - 1]);
    } 
    Mul(A, B, A, n, n);
    for(int i = 0;i < n; ++i) {
        A[i] = mul(mul(A[i + n - 1], Pow(b, i)), mul(ivf[i], Pow(c, (long long)i * i))); 
        A[i + n - 1] = 0;
    }
    for(int i = 0;i < (n << 1); ++i)
        B[i] = Pow(c, -(long long)(n - i) * (n - i));
    Mul(A, B, A, n, n << 1);
    for(int i = 0;i < n; ++i)
        printf("%d\n", mul(A[i + n], Pow(c, (long long)i * i)));
    return 0;
}