Bluestein's Algorithm

用下降幂搞搞事情

\[ \frac{m(m-1) + k(k+1) -(m-k)(m-k-1)}{2} \\ =\frac{m^2 - m + k^2 + k - m^2 - k^2 + 2mk - k + m}{2} \\ = mk \\ y_{m} = \sum_{k=0}^n A_{k}\omega^{mk}\\=\omega^{\frac{m(m-1)}{2}}\sum_{k=0}^n A_k\omega^{\frac{k(k+1)}{2}}\omega^{\frac{-(m-k)(m-k-1)}{2}} \]

后面是一个卷积形式，可以快速实现求 \(f(x),f(x^1), f(x^2) .. f(x^n)\) 以及循环卷积意义下的多项式快速幂。

例题 UOJ500 任意基DFT

观察发现

\[q_m = q_0x^m+\sum_{i=0}^{m-1} xy^i =(\frac{y}{x-1}+q_0)x^m-\frac{y}{x-1}\\ \]

​ 设 \(c = \frac{y}{x-1}\) ，那么有

​

\[f(q_m) = \sum_{i=0}^n a_i((c+q_0)x^m-c)^i \\ =\sum_{i=0}^n a_i\sum_{j=0}^i{i\choose j}(c+q_0)^jx^{mj}(-c)^{i-j} \\ =\sum_{i=0}^n(c+q_0)^ix^{mi}\sum_{j=i}^n{i \choose j}(-c)^{i-j}a_i \]

令 \(b_j = (c+q_0)^i \sum_{j=i}^n{i \choose j}(-c)^{i-j}a_i\) ，这是一个卷积形式，可以用一遍 \(\text{FFT}\) 求出所有的 \(b\) ，那么

\[f(q_m) = \sum_{i=0}^nb_i x^{mi} \]

再跑一遍 bluestein's 对 \(1-m\) 求解即可，注意用线性递推求幂次来优化常数。

code

/*program by mangoyang*/
#pragma GCC optimize("Ofast", "inline")
#include<bits/stdc++.h>
#define inf (0x7f7f7f7f)
#define Max(a, b) ((a) > (b) ? (a) : (b))
#define Min(a, b) ((a) < (b) ? (a) : (b))
typedef long long ll;
using namespace std;
template <class T>
inline void read(T &x){
    int ch = 0, f = 0; x = 0;
    for(; !isdigit(ch); ch = getchar()) if(ch == '-') f = 1;
    for(; isdigit(ch); ch = getchar()) x = x * 10 + ch - 48;
    if(f) x = -x;
}
const int N = 1 << 21, mod = 998244353, G = 3;
int a[N], b[N], A[N], B[N], C[N], d[N], e[N], js[N], inv[N], n, q, q0, qx, qy;
inline void up(int &x, int y){
	x = x + y >= mod ? x + y - mod : x + y;
}
inline int Pow(int a, int b){
	int ans = 1;
	for(; b; b >>= 1, a = 1ll * a * a % mod)
		if(b & 1) ans = 1ll * ans * a % mod;
	return ans;
}
namespace poly{
	int rev[N], len, lg;
	inline void init(int lenth){
		for(len = 1, lg = 0; len <= lenth; len <<= 1, lg++);
		for(int i = 0; i < len; i++)
			rev[i] = (rev[i>>1] >> 1) | ((i & 1) << (lg - 1));
	}
	inline void dft(int *a, int sgn){
		for(int i = 0; i < len; i++)
			if(i < rev[i]) swap(a[i], a[rev[i]]);
		for(int k = 2; k <= len; k <<= 1){
			int w = Pow(G, (mod - 1) / k);
			if(sgn == -1) w = Pow(w, mod - 2);
			for(int i = 0, now; i < len; i += k){
				now = 1;
				for(int j = i; j < i + (k >> 1); j++){
					int x = a[j], y = 1ll * a[j+(k>>1)] * now % mod;
					a[j] = x + y >= mod ? x + y - mod : x + y;
					a[j+(k>>1)] = x - y < 0 ? x - y + mod : x - y;
					now = 1ll * now * w % mod;
				}
			}
		}
		if(sgn == -1){
			int x = Pow(len, mod - 2);
			for(int i = 0; i < len; i++) a[i] = 1ll * a[i] * x % mod;
		}
	}
}
int main(){
	read(n), read(q), js[0] = 1, inv[0] = 1;
	for(int i = 1; i <= n; i++){
		js[i] = 1ll * js[i-1] * i % mod;
		inv[i] = Pow(js[i], mod - 2);
	}
	for(int i = 0; i <= n; i++) read(a[i]);
	read(q0), read(qx), read(qy);
	int c = 1ll * qy * Pow(qx - 1, mod - 2) % mod;
	for(int i = 0, ss = 1; i <= n; i++){
		A[n-i] = 1ll * js[i] * a[i] % mod;
		B[i] = 1ll * inv[i] * ss % mod;
		ss = 1ll * ss * (mod - c) % mod;
	}
	poly::init(n + n + 1);
	poly::dft(A, 1), poly::dft(B, 1);
	for(int i = 0; i < poly::len; i++)
		C[i] = 1ll * A[i] * B[i] % mod;
	poly::dft(C, -1);
	for(int i = 0; i <= n; i++){
		b[i] = 1ll * C[n-i] * inv[i] % mod;
		b[i] = 1ll * b[i] * Pow(c + q0, i) % mod;
	}
	memset(A, 0, sizeof(A));
	memset(B, 0, sizeof(B));
	memset(C, 0, sizeof(C));
	d[0] = 1, e[0] = 1;
	int ix = Pow(qx, mod - 2);
	for(int i = 1, ss = 1, tt = 1; i <= n + q + 1; i++){
		ss = 1ll * qx * ss % mod;
		tt = 1ll * ix * tt % mod;
		d[i] = 1ll * d[i-1] * ss % mod;
		e[i] = 1ll * e[i-1] * tt % mod;
	}
	for(int i = 0; i <= n; i++)
		A[i] = 1ll * b[i] * d[i] % mod;
	for(int i = -n + 1; i <= q; i++)
		B[i+n] = i <= 0 ? e[abs(i)] : e[i-1];
	poly::init(n + n + q);
	poly::dft(A, 1), poly::dft(B, 1);
	for(int i = 0; i < poly::len; i++)
		C[i] = 1ll * A[i] * B[i] % mod;
	poly::dft(C, -1);
	int ans = 0;
	for(int i = 1; i <= q; i++)
		ans ^= 1ll * C[i+n] * d[i-1] % mod;
	cout << ans << endl;
	return 0;
}