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
;
cd omega
, omegaInv
;
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
, b
, c
, d
;

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
, b, c, d, p
, g
, f
, c2
;
ll inv
, fac
, facInv
;
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]);
}