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]);
}