Bluestein's Algorithm

Bluestein’s Algorithm is an algorithm that computes $f(z^0),f(z^1),\dots ,f(z^{m-1})$ for a polynomial $f(x)={\sum }_{i=0}^{n-1}{a}_i{x}^i$ in $\mathcal{O}((n+m)\log (n+m))$ time and $\mathcal{O}(n+m)$ space.

Algorithm

Lemma

Let

$$P(x)=\sum _{i=0}^{n-1}{z}^{-\frac{i^2}{2}}{a}_i{x}^{n-i-1},Q(x)=\sum _{i=0}^{n+m-1}{z}^{\frac{i^2}{2}}{x}^i$$

then

$$\forall i\in \{0,1,\dots ,m-1\},f(z^i)=z^{-\frac{i^2}{2}}[x^{n-1+i}](P(x)Q(x))$$

Proof

$$\begin{array}{rl}{z}^{-\frac{i^2}{2}}[x^{n-1+i}](P(x)Q(x))& =z^{-\frac{i^2}{2}}\sum _{j=0}^{n-1}[x^{n-j-1}]P(x)[x^{i+j}]Q(x)\\ & =\sum _{j=0}^{n-1}{a}_j{z}^{\frac{(i+j{)}^2}{2}-\frac{i^2}{2}-\frac{j^2}{2}}\\ & =\sum _{j=0}^{n-1}{a}_j{z}^{ij}\\ & =f(z^i)\end{array}$$

0. Let $N=2^{\lceil {\log }_2(2n+m-1)\rceil }$. Apply the Fast Fourier Transform to find $P({\omega }_N^0),P({\omega }_N^1),\dots ,P({\omega }_N^{N-1})$ and $Q({\omega }_N^0),Q({\omega }_N^1),\dots ,Q({\omega }_N^{N-1})$.
1. Use the results from 0 to find $P({\omega }_N^0)Q({\omega }_N^0),P({\omega }_N^1)Q({\omega }_N^1),\dots ,P({\omega }_N^{N-1})Q({\omega }_N^{N-1})$.
2. Apply the Inverse Fast Fourier Transform to find $P(x)Q(x)$.
3. For each $i$ in $\{0,1,\dots ,m-1\}$, apply the lemma to find $f(z^i)$.

This algorithm solves the problem in $\mathcal{O}((n+m)\log (n+m))$ time and $\mathcal{O}(n+m)$ space.

std::vector<std::complex<double>> bluestein(int n, const std::vector<std::complex<double>> &a, std::complex<double> z, int m) {
	auto f = std::ranges::to<std::vector>(std::views::iota(0, n) | std::views::transform([&](int i) -> std::complex<double> {
		return a[i] / std::pow(z, 1. * i * i / 2);
	}) | std::views::reverse), g = std::ranges::to<std::vector>(std::views::iota(0, n + m) | std::views::transform([&](int i) -> std::complex<double> {
		return std::pow(z, 1. * i * i / 2);
	}));
 
	int k = std::bit_ceil<u32>(2 * n + m - 1);
	f.resize(k, 0), g.resize(k, 0);
	fft(k, f), fft(k, g);
	for (int i = 0; i < k; i++) {
		f[i] *= g[i];
	}
	ifft(k, f);
 
	return std::ranges::to<std::vector>(std::views::iota(0, m) | std::views::transform([&](int i) -> std::complex<double> {
		return f[n - 1 + i] / std::pow(z, 1. * i * i / 2);
	}));
}