hdu 4656 Evaluation [任意模数fft trick]

hdu 4656 Evaluation

题意：给出\(n,b,c,d,f(x) = \sum_{i=1}^{n-1} a_ix^i\)，求\(f(b\cdot c^{2k}+d):0\le k < n\)

取模\(10^6+3\)

---

昨天刚看过《具体数学》上求和一章

代入\(b\cdot c^{2k}+d\)然后展开，交换求和顺序，得到
\[ f(k) = \sum_{j=0}^n \frac{b^j c^{2kj}}{j!} \sum_{i=j}^n a_i i! \frac{d^{i-j}}{(i-j)!} \]
后面的式子反转后变成卷积，fft预处理\(p_j = \sum_{i=j}^n a_i i! \frac{d^{i-j}}{(i-j)!}\)

那么我们要求
\[ f(k) = \sum_{j=0}^n \frac{b^j c^{2kj} p_j}{j!} \]
到这里就不会了，看了叉姐的课件，很神的一步来转化成卷积
\[ (k-j)^2 = k^2 + j^2 - 2kj \]
代入之后得到
\[ f_k = c^{k^2}\sum_{j=0}^n \frac{b^j c^{j^2}p_j}{j!c^{(k-j)^2}} \]
再次fft就行了

注意，指数带着平方可以为负，所以我们先乘上\(x^n\)

然后还要用拆系数fft.... 完虐ntt怒拿rank1

#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
#include <ctime>
using namespace std;
typedef long long ll;
const int N = (1<<18) + 5, mo = 1e6+3, P = 1e6+3;
const double PI = acos(-1.0);
inline int read() {
    char c=getchar(); int x=0,f=1;
    while(c<'0'||c>'9') {if(c=='-')f=-1;c=getchar();}
    while(c>='0'&&c<='9') {x=x*10+c-'0';c=getchar();}
    return x*f;
}

struct meow{
    double x, y;
    meow(double a=0, double b=0):x(a), y(b){}
};
meow operator +(meow a, meow b) {return meow(a.x+b.x, a.y+b.y);}
meow operator -(meow a, meow b) {return meow(a.x-b.x, a.y-b.y);}
meow operator *(meow a, meow b) {return meow(a.x*b.x-a.y*b.y, a.x*b.y+a.y*b.x);}
meow conj(meow a) {return meow(a.x, -a.y);}
typedef meow cd;

namespace fft {
    int maxlen, rev[N];
    cd omega[N], omegaInv[N];
    void init(int lim) {
        maxlen = 1; while(maxlen < lim) maxlen<<=1;
        for(int i=0; i<maxlen; i++) {
            omega[i] = cd(cos(2*PI/maxlen*i), sin(2*PI/maxlen*i));
            omegaInv[i] = conj(omega[i]);
        }
    }
    void dft(cd *a, int n, int flag) {
        cd *w = flag == 1 ? omega : omegaInv;
        int k = 0; while((1<<k) < n) k++;
        for(int i=0; i<n; i++) {
            rev[i] = (rev[i>>1]>>1) | ((i&1)<<(k-1));
            if(i < rev[i]) swap(a[i], a[rev[i]]);
        }
        for(int l=2; l<=n; l<<=1) {
            int m = l>>1;
            for(cd *p = a; p != a+n; p += l) 
                for(int k=0; k<m; k++) {
                    cd t = w[maxlen/l*k] * p[k+m];
                    p[k+m] = p[k] - t;
                    p[k] = p[k] + t;
                }
        }
        if(flag == -1) for(int i=0; i<n; i++) a[i].x /= n;
    }

    cd a[N], b[N], c[N], d[N];
    
    void mul_any(int *x, int *y, int lim) {
        int n = maxlen;
        for(int i=0; i<lim; i++) {
            a[i] = cd(x[i]>>15), b[i] = cd(x[i]&32767);
            c[i] = cd(y[i]>>15), d[i] = cd(y[i]&32767);
        }
        for(int i=lim; i<n; i++) a[i] = b[i] = c[i] = d[i] = cd();
        dft(a, n, 1); dft(b, n, 1); dft(c, n, 1); dft(d, n, 1);
        for(int i=0; i<n; i++) {
            cd _a = a[i], _b = b[i], _c = c[i], _d = d[i];
            a[i] = _a * _c;
            b[i] = _a * _d + _b * _c;
            c[i] = _b * _d;
        }
        dft(a, n, -1); dft(b, n, -1); dft(c, n, -1);
        for(int i=0; i<lim; i++) x[i] = ( (ll(a[i].x + 0.5) %mo <<30) + (ll(b[i].x + 0.5) %mo <<15) + ll(c[i].x + 0.5)%mo) %mo;
    }
}

inline ll Pow(ll a, int b) {
    ll ans = 1;
    for(; b; b>>=1, a=a*a%P)
        if(b&1) ans=ans*a%P;
    return ans;
}
int n, a[N], b, c, d, p[N], g[N], f[N], c2[N];
ll inv[N], fac[N], facInv[N];
int main() {
    freopen("in", "r", stdin);
    n=read(); b=read(); c=read(); d=read();
    for(int i=0; i<n; i++) a[i] = read();
    fft::init(n+n+1);
    
    inv[1] = fac[0] = facInv[0] = 1;
    for(int i=1; i<=n; i++) {
        if(i != 1) inv[i] = (P-P/i) * inv[P%i] %P;
        fac[i] = fac[i-1] * i %P;
        facInv[i] = facInv[i-1] * inv[i] %P;
    }

    ll mi = 1;
    for(int i=0; i<=n; i++, mi = mi * d %P) p[i] = a[n-i] * fac[n-i] %P, g[i] = mi * facInv[i] %P;
    fft::mul_any(p, g, n+1);
    for(int i=0; i<=n>>1; i++) swap(p[i], p[n-i]);

    mi = 1;
    memset(g, 0, sizeof(g));
    for(int i=0; i<=n; i++, mi = mi * b %P) {
        int t = Pow(c, (ll) i * i %(P-1));
        f[i] = mi * t %P * p[i] %P * facInv[i] %P;
        if(i != n) g[i + n] = g[n - i] = Pow(t, P-2);
        //if(i != n) g[i + n-1] = g[n-1 - i] = Pow(t, P-2);
        c2[i] = t; 
    }
    fft::mul_any(f, g, n+n+1);
    for(int i=0; i<n; i++) f[i+n] = (ll) f[i+n] * c2[i] %P, printf("%d\n", f[i+n]);
    //for(int i=0; i<n; i++) f[i+n-1] = (ll) f[i+n-1] * c2[i] %P, printf("%d\n", f[i+n-1]);
}