Pipeline-generated

How statement and proof provenance work

The first chip identifies the source of the statement or construction; the second identifies the source of its local proof or verification.

Literature-sourced: the exact statement appears in a cited source; only wording and notation differ.
AI-adapted: a semantically identical restatement of literature-sourced material, modulo indexing, notation, and boundary cases adopted by the library.
AI-generated: a genuinely novel statement formulated by AI, with no source for the claim itself.

These labels describe origin, not correctness: citations and verification chips remain separate evidence.

Computed fast transforms and convolution

1 · Prerequisites

2 · Summary

This page extends the canonical finite Fourier transform and radix-two FFT with computed error bounds, arbitrary-length reductions, trigonometric transforms, and convolution. Its radix-two stability estimate assumes round-to-nearest, normal-range successful operations, and a certified error bound for every twiddle. The product-stage constant is explicit: $η = 10 u + ϵ_{t w}$ ; the $O (u \log N)$ shorthand requires $ϵ_{t w} = O (u)$ and $u \log N$ small. A phase-evaluation budget that cannot certify its requested enclosure returns inconclusive. The interval evaluator's arctangent definition is shown to use the same canonical $π$ as the Fourier and chirp formulas.

Mixed-radix Cooley–Tukey and Bluestein supply exact arbitrary-length reductions. Bluestein pads to a supported power of two $M \geq 2 N - 1$ and records chirp-phase error separately from the FFT-convolution error. Arithmetic counts exclude integer factor search, interval phase evaluation, memory traffic, and bit complexity; no general mixed-radix stability theorem is asserted here.

DCT-II uses the even extension and the stated endpoint normalization; DST-I uses the odd extension and its fixed endpoint zeros. Both inverses are proved directly, and the DCT-II recovers Chebyshev coefficients at first-kind roots. Linear convolution uses zero padding to prevent circular aliasing. Its error bound separates forward transforms, pointwise products, inverse transform, and final scaling, conditional on successful operations.

3 · Logical flowchart

4 · Definitions, theorems and proofs

DefinitionDefinition: Literature-sourcedProof: Not applicablejudge pass (gpt-6-luna)Open item page →

DCT-II and DST-I with fixed even and odd extension conventions

Definition

For $N \geq 1$ and $x = (x_{0}, \dots, x_{N - 1}) \in R^{N}$ , the DCT-II is the vector $X \in R^{N}$ with $X_{k} = \sum_{j = 0}^{N - 1} x_{j} \cos (π (j + 1 / 2) k / N)$ for $0 \leq k < N$ . Its even extension $E_{N} x \in R^{2 N}$ is $(E_{N} x)_{j} = x_{j}$ and $(E_{N} x)_{2 N - 1 - j} = x_{j}$ for $0 \leq j < N$ .

For $N \geq 1$ and $x = (x_{1}, \dots, x_{N}) \in R^{N}$ , the DST-I is the vector $Y \in R^{N}$ with $Y_{k} = \sum_{j = 1}^{N} x_{j} \sin (π j k / (N + 1))$ for $1 \leq k \leq N$ . Its odd extension $O_{N} x \in R^{2 (N + 1)}$ has $(O_{N} x)_{0} = (O_{N} x)_{N + 1} = 0$ , $(O_{N} x)_{j} = x_{j}$ , and $(O_{N} x)_{2 (N + 1) - j} = - x_{j}$ for $1 \leq j \leq N$ .

These definitions fix the transform type, index range, extension symmetry, and endpoint normalization. Their orthogonality, inverse formulas, and relation to the unitary DFT are proved in separate results.

LemmaStatement: Literature-sourcedProof: AI-generatedOpen item page →

Bluestein’s chirp identity embeds an arbitrary DFT in a padded cyclic convolution

Statement

Let $N \geq 1$ , $M \geq 2 N - 1$ , and $x \in C^{N}$ . Put $α_{j} = e^{- π i j^{2} / N}$ for $0 \leq j < N$ , $β_{t} = e^{π i t^{2} / N}$ for $- (N - 1) \leq t \leq N - 1$ , and $γ_{k} = e^{- π i k^{2} / N}$ for $0 \leq k < N$ . Embed $a_{j} = α_{j} x_{j}$ for $0 \leq j < N$ and $a_{j} = 0$ otherwise in $C^{M}$ . Embed $b_{[t]_{M}} = β_{t}$ for $- (N - 1) \leq t \leq N - 1$ and set its other entries to zero; this is well-defined because these $2 N - 1$ integers have distinct classes modulo $M$ . Let $c = a *_{M} b$ be the unnormalised cyclic convolution on $Z / M Z$ . Then, for $0 \leq k < N$ ,

(F_{N} x) (k) = N^{- 1 / 2} γ_{k} c_{k},

where $γ_{k} = e^{- π i k^{2} / N}$ . Moreover $c = \sqrt{M} F_{M}^{- 1} ((F_{M} a) ⊙ (F_{M} b))$ , where $⊙$ is coordinatewise multiplication. The identity is exact for every $N$ ; padding to a radix-two $M$ is an implementation choice, not a mathematical restriction.

Facts & Assumptions

Given: Integers $N \geq 1$ , $M \geq 2 N - 1$ , a vector $x \in C^{N}$ , and the definitions of $F_{N}$ , $F_{M}$ , and cyclic convolution.

[F1]

The unnormalised DFT is $X_{k} = \sum_{j = 0}^{N - 1} x_{j} e^{- 2 π i j k / N}$ and $F_{N} x = N^{- 1 / 2} X$ (The unnormalised engineering DFT and its conversion to the unitary transform, The unitary discrete Fourier transform on $Z / N Z$ ).

[F2]

Cyclic convolution is $(a *_{M} b)_{k} = \sum_{j = 0}^{M - 1} a_{j} b_{[k - j]_{M}}$ , and its unitary DFT satisfies $F_{M} (a *_{M} b) = \sqrt{M} (F_{M} a) ⊙ (F_{M} b)$ (The unnormalised cyclic convolution on $Z / N Z$ , The DFT turns cyclic convolution into a scaled pointwise product).

[F3]

The unitary DFT on $Z / M Z$ is invertible with inverse given by the positive-sign unitary DFT (Finite Fourier inversion for the unitary transform on $Z / N Z$ ).

[F4]

The map $t \mapsto [t]_{M}$ is injective on the integer interval $[- (N - 1), N - 1]$ when $M \geq 2 N - 1$ ; complex exponentials satisfy their addition law (The congruence class $[a]_{n}$ and the quotient set $Z / n$ , For every natural $n$ , $(Z / n, +)$ is an abelian group, multiplication is a commutative monoid operation, and both distributive laws hold, For $n \geq 1$ , every class in $Z / n$ has one representative $r$ with $0 \leq r < n$ , so $∣ Z / n ∣ = n$ ; while $Z / 0$ is in bijection with $Z$ , $\exp (z + w) = \exp z \exp w$ , and the complex exponential extends the real exponential, $C = R [x] / (x^{2} + 1)$ is a field, every element is uniquely $a + b i$ , and every nonzero element has inverse $(a - b i) / (a^{2} + b^{2})$ ).

[F5]

Positive rational powers satisfy $(r s)^{1 / 2} = r^{1 / 2} s^{1 / 2}$ for positive $r, s$ (Rational powers $a^{r}$ of a positive base, Laws of rational exponents).

Proof

1.1givenF4choose

If two integers $t, t^{'}$ in $[- (N - 1), N - 1]$ have the same class modulo $M$ , then $M$ divides $t - t^{'}$ . But $∣ t - t^{'} ∣ \leq 2 N - 2 < M$ ; the only multiple of $M$ in this range is zero, so $t = t^{'}$ . Thus the definition of $b$ assigns at most one value to each residue. For each $0 \leq k < N$ and $0 \leq j < N$ , the difference $k - j$ lies in $[- (N - 1), N - 1]$ , so the cyclic-convolution sum reduces exactly to $c_{k} = \sum_{j = 0}^{N - 1} α_{j} x_{j} β_{k - j}$ .

2.1step 1.1F1F4algebra

For every such $j, k$ , $- j^{2} + (k - j)^{2} - k^{2} = - 2 k j$ , hence $α_{j} β_{k - j} γ_{k} = e^{- 2 π i j k / N}$ . Summing this identity against $x_{j}$ and using [F1] gives $γ_{k} c_{k} = X_{k}$ and therefore $(F_{N} x) (k) = N^{- 1 / 2} γ_{k} c_{k}$ .

3.1F2F3F5algebra∎

Applying [F2] to $a, b$ gives $F_{M} c = \sqrt{M} ((F_{M} a) ⊙ (F_{M} b))$ . Applying the positive-sign inverse in [F3] and pulling the scalar $\sqrt{M}$ through its finite sum gives $c = \sqrt{M} F_{M}^{- 1} ((F_{M} a) ⊙ (F_{M} b))$ , as claimed.

LemmaStatement: Literature-sourcedProof: AI-generatedjudge pass (gpt-6-luna)Open item page →

Mixed-radix Cooley–Tukey factorization with fixed index maps

Statement

Let $r, s \geq 2$ , $N = r s$ , and $x = (x_{j})_{j = 0}^{N - 1} \in C^{N}$ . For $0 \leq u < r$ define $x_{v}^{(u)} = x_{u + r v}$ for $0 \leq v < s$ . For $0 \leq k_{2} < s$ set $y_{u}^{(k_{2})} = e^{- 2 π i u k_{2} / N} (F_{s} x^{(u)}) (k_{2})$ for $0 \leq u < r$ . Then for $0 \leq k_{1} < r$ and $0 \leq k_{2} < s$ ,

(F_{N} x) (k_{2} + s k_{1}) = (F_{r} y^{(k_{2})}) (k_{1}) .

Here every $F_{q}$ is the negative-sign unitary DFT with normalization $q^{- 1 / 2}$ from The unitary discrete Fourier transform on $Z / N Z$ . Thus the inner transforms are the $r$ length- $s$ transforms of strided input subsequences; the twiddle multiplication is explicit; and the outer transforms are the $s$ length- $r$ transforms. No coprimality of $r$ and $s$ is assumed.

Facts & Assumptions

Given: Integers $r, s \geq 2$ , $N = r s$ , a vector $x \in C^{N}$ , and indices $k_{1} \in {0, \dots, r - 1}$ and $k_{2} \in {0, \dots, s - 1}$ .

[F1]

The unitary DFT is $(F_{q} z) (k) = q^{- 1 / 2} \sum_{v = 0}^{q - 1} z_{v} e^{- 2 π i v k / q}$ for $q \geq 1$ (The unitary discrete Fourier transform on $Z / N Z$ ).

[F2]

Every integer $j$ with $0 \leq j < r s$ is uniquely $j = u + r v$ with $0 \leq u < r$ and $0 \leq v < s$ ; every $k$ with $0 \leq k < r s$ is uniquely $k = k_{2} + s k_{1}$ in the stated ranges (For every natural $n$ , $(Z / n, +)$ is an abelian group, multiplication is a commutative monoid operation, and both distributive laws hold, For $n \geq 1$ , every class in $Z / n$ has one representative $r$ with $0 \leq r < n$ , so $∣ Z / n ∣ = n$ ; while $Z / 0$ is in bijection with $Z$ ).

[F3]

Finite sums over these index sets can be reindexed by the bijections in [F2], split into finite sums, and reordered (A finite sum in a commutative monoid indexed by an arbitrary finite set, Finite commutative-monoid sums are invariant under bijective reindexing, split over disjoint unions, and satisfy the finite Fubini rule).

[F4]

For complex exponentials, $e^{z + w} = e^{z} e^{w}$ and $e^{z} = 1$ exactly when $z \in 2 π i Z$ , so $e^{2 π i n} = 1$ for every integer $n$ ; positive rational powers obey their exponent laws ( $\exp (z + w) = \exp z \exp w$ , and the complex exponential extends the real exponential, $\ker (\exp) = 2 π i Z$ , and $\exp z = \exp w$ exactly when $z - w \in 2 π i Z$ , $C = R [x] / (x^{2} + 1)$ is a field, every element is uniquely $a + b i$ , and every nonzero element has inverse $(a - b i) / (a^{2} + b^{2})$ , Rational powers $a^{r}$ of a positive base, Laws of rational exponents).

Proof

1.1givenF1F2F3

By [F2], the pair $(u, v)$ with $j = u + r v$ runs exactly once through the $N$ input indices, and $(k_{1}, k_{2})$ with $k = k_{2} + s k_{1}$ runs exactly once through the $N$ output indices. Therefore the definition [F1] gives $(F_{N} x) (k_{2} + s k_{1}) = (r s)^{- 1 / 2} \sum_{u = 0}^{r - 1} \sum_{v = 0}^{s - 1} x_{u + r v} e^{- 2 π i (u + r v) (k_{2} + s k_{1}) / (r s)}$ .

2.1step 1.1F3F4algebra

Expanding the exponent gives $(u + r v) (k_{2} + s k_{1}) / (r s) = u k_{2} / (r s) + u k_{1} / r + v k_{2} / s + v k_{1}$ . Since $v k_{1}$ is an integer, its exponential factor is $1$ by [F4]. Also $(r s)^{- 1 / 2} = r^{- 1 / 2} s^{- 1 / 2}$ by [F4]. Regrouping the finite sums in [F3] now yields $(r s)^{- 1 / 2} \sum_{u = 0}^{r - 1} e^{- 2 π i u k_{2} / N} e^{- 2 π i u k_{1} / r} \sum_{v = 0}^{s - 1} x_{u + r v} e^{- 2 π i v k_{2} / s}$ .

3.1step 2.1F1given∎

The inner sum in step 2.1 is $s^{1 / 2} (F_{s} x^{(u)}) (k_{2})$ by [F1]. The remaining coefficient and sum are exactly $r^{- 1 / 2} \sum_{u = 0}^{r - 1} y_{u}^{(k_{2})} e^{- 2 π i u k_{1} / r} = (F_{r} y^{(k_{2})}) (k_{1})$ , again by [F1] and the definition of $y^{(k_{2})}$ . This is the claimed identity.

LemmaStatement: Literature-sourcedProof: AI-generatedOpen item page →

Zero padding embeds linear convolution in cyclic convolution without aliasing

Statement

Let $a_{0}, \dots, a_{n - 1}$ and $b_{0}, \dots, b_{m - 1}$ be finite complex sequences, and let $M \geq n + m - 1$ . Extend both sequences by zero to indices $0, \dots, M - 1$ . For every $0 \leq k < n + m - 1$ , their length- $M$ cyclic convolution equals the ordinary linear convolution $c_{k} = \sum_{j = 0}^{n - 1} a_{j} b_{k - j}$ , with $b_{q} = 0$ outside $0 \leq q < m$ . For every remaining cyclic index $n + m - 1 \leq k < M$ , the cyclic convolution is zero. Thus the first $n + m - 1$ cyclic-convolution outputs are exactly the full linear convolution, with no wraparound alias.

Facts & Assumptions

Given: Positive sequence lengths $n, m$ , an integer $M \geq n + m - 1$ , and the zero-padded vectors $a, b \in C^{M}$ .

[F1]

The length- $M$ cyclic convolution is $(a *_{M} b)_{k} = \sum_{j = 0}^{M - 1} a_{j} b_{[k - j]_{M}}$ (The unnormalised cyclic convolution on $Z / N Z$ ).

[F2]

By the given zero extension, the supports of the padded vectors are contained in ${0, \dots, n - 1}$ and ${0, \dots, m - 1}$ , respectively. Each residue modulo $M$ has exactly one representative in ${0, \dots, M - 1}$ (For $n \geq 1$ , every class in $Z / n$ has one representative $r$ with $0 \leq r < n$ , so $∣ Z / n ∣ = n$ ; while $Z / 0$ is in bijection with $Z$ ).

Proof

1.1givenF1F2choose

A nonzero summand in [F1] has $0 \leq j < n$ and its second index is a residue in ${0, \dots, m - 1}$ . If $k < n + m - 1$ and $j > k$ , then $k - j < 0$ and $M + k - j \geq M - (n - 1) \geq m$ , so $[k - j]_{M}$ lies outside the support of $b$ and the summand is zero. If $j \leq k$ , then $k - j \geq 0$ ; whenever $k - j \geq m$ the summand is again zero, and otherwise the cyclic representative equals the ordinary index $k - j$ . Therefore the nonzero terms of [F1] are exactly the linear-convolution terms.

2.1step 1.1F1F2algebra∎

Every index sum $j + q$ with $a_{j} b_{q} \neq 0$ satisfies $0 \leq j + q \leq n + m - 2 < M$ . Hence no such sum wraps modulo $M$ . For $k \geq n + m - 1$ , no supported pair has sum congruent to $k$ modulo $M$ : a sum is at most $n + m - 2$ , and its difference from $k$ has magnitude less than $M$ for the output indices $0 \leq k < M$ . Thus all cyclic terms vanish there. This proves the asserted zero tail and completes the claim.

LemmaStatement: Literature-sourcedProof: AI-generatedOpen item page →

Six-operation complex multiplication has a normwise roundoff bound

Statement

Let $u$ be the unit roundoff of the operationwise model in Operationwise real and complex floating point arithmetic with rounded square roots and fused multiply-add, with $0 < u \leq 1 / 16$ . For $z, w \in C$ , compute ${fl}_{C} (z w)$ from the four real products, one real subtraction, and one real addition specified in Operationwise real and complex floating point arithmetic with rounded square roots and fused multiply-add. If every nonzero intermediate is normal and no operation overflows, then

∣ {fl}_{C} (z w) - z w ∣ \leq 3 u ∣ z ∣ ∣ w ∣ .

Facts & Assumptions

Given: $z = z_{r} + i z_{i}$ , $w = w_{r} + i w_{i}$ , the operationwise model with unit roundoff $0 < u \leq 1 / 16$ , and success of all six real operations.

[F1]

A complex product computes four separately rounded real products and then rounded real/imaginary sums, for six real operations; successful operations have exact-zero behavior or relative error at most $u$ (Operationwise real and complex floating point arithmetic with rounded square roots and fused multiply-add, The standard relative floating-point model and unit roundoff).

[F2]

$∣ z ∣ = (z_{r}^{2} + z_{i}^{2})^{1 / 2}$ and likewise for $w$ (Real and imaginary parts, complex conjugation, and modulus); Cauchy–Schwarz in $R^{2}$ gives $∣ z_{r} w_{r} ∣ + ∣ z_{i} w_{i} ∣ \leq ∣ z ∣ ∣ w ∣$ and $∣ z_{r} w_{i} ∣ + ∣ z_{i} w_{r} ∣ \leq ∣ z ∣ ∣ w ∣$ (Absolute value in an ordered field, Cauchy–Schwarz: $∣ ⟨ u, v ⟩ ∣ \leq ∥ u ∥ ∥ v ∥$ , with equality exactly for linearly dependent vectors).

Proof

1.1givenF1algebra

Write $p_{1} = fl (z_{r} w_{r})$ and $p_{2} = fl (z_{i} w_{i})$ . By [F1], $∣ p_{1} - z_{r} w_{r} ∣ \leq u ∣ z_{r} w_{r} ∣$ and $∣ p_{2} - z_{i} w_{i} ∣ \leq u ∣ z_{i} w_{i} ∣$ . Rounding their difference adds error at most $u ∣ p_{1} - p_{2} ∣ \leq u (1 + u) (∣ z_{r} w_{r} ∣ + ∣ z_{i} w_{i} ∣)$ . Therefore the real-part error is at most $(2 u + u^{2}) (∣ z_{r} w_{r} ∣ + ∣ z_{i} w_{i} ∣)$ . The same calculation for the imaginary part gives $(2 u + u^{2}) (∣ z_{r} w_{i} ∣ + ∣ z_{i} w_{r} ∣)$ .

2.1step 1.1F2algebra∎

Applying [F2] to both component sums shows that each component error is at most $(2 u + u^{2}) ∣ z ∣ ∣ w ∣$ . The complex modulus of the two-component error is consequently at most $\sqrt{2} (2 u + u^{2}) ∣ z ∣ ∣ w ∣$ . Since $u \leq 1 / 16$ , $\sqrt{2} (2 + u) u \leq \sqrt{2} (33 / 16) u < 3 u$ , proving the stated bound.

TheoremStatement: Literature-sourcedProof: AI-generatedjudge pass (gpt-6-luna)Open item page →

DCT-II is orthogonal with the stated endpoint normalization and even extension

Statement

Let $C_{k j} = \cos (π (j + 1 / 2) k / N)$ for $0 \leq j, k < N$ . Then $C C^{T} = diag (N, N / 2, \dots, N / 2)$ , so the inverse DCT-II is $x_{j} = N^{- 1} (X_{0} + 2 \sum_{k = 1}^{N - 1} X_{k} \cos (π (j + 1 / 2) k / N))$ . For the even extension $E_{N} x$ of DCT-II and DST-I with fixed even and odd extension conventions, the unitary DFT also satisfies $(F_{2 N} E_{N} x) (k) = \sqrt{2 / N} e^{π i k / (2 N)} X_{k}$ for $0 \leq k < N$ .

Facts & Assumptions

Given: An integer $N \geq 1$ , a vector $x \in R^{N}$ , the DCT-II matrix $C$ , and the even extension $E_{N} x$ .

[F1]

$X_{k} = \sum_{j = 0}^{N - 1} x_{j} \cos (π (j + 1 / 2) k / N)$ and $(E_{N} x)_{j} = (E_{N} x)_{2 N - 1 - j} = x_{j}$ (DCT-II and DST-I with fixed even and odd extension conventions).

[F2]

$(F_{2 N} z) (k) = (2 N)^{- 1 / 2} \sum_{m = 0}^{2 N - 1} z_{m} e^{- 2 π i m k / (2 N)}$ (The unitary discrete Fourier transform on $Z / N Z$ ).

[F3]

$e^{i θ} = \cos θ + i \sin θ$ , $e^{u + v} = e^{u} e^{v}$ , and the complex field laws hold ( $\exp (x + i y) = e^{x} (\cos y + i \sin y)$ , $∣ \exp (x + i y) ∣ = e^{x}$ , and $e^{i π} + 1 = 0$ , $\exp (z + w) = \exp z \exp w$ , and the complex exponential extends the real exponential, $C = R [x] / (x^{2} + 1)$ is a field, every element is uniquely $a + b i$ , and every nonzero element has inverse $(a - b i) / (a^{2} + b^{2})$ ).

[F4]

For real $u, v$ , $\cos u \cos v = \frac{1}{2} (\cos (u - v) + \cos (u + v))$ (Product-to-sum and sum-to-product identities).

[F5]

$e^{z} = 1$ exactly when $z \in 2 π i Z$ , and $\sin t = 0$ exactly when $t \in π Z$ ( $\ker (\exp) = 2 π i Z$ , and $\exp z = \exp w$ exactly when $z - w \in 2 π i Z$ , The zero sets of sine and cosine and the least positive common period 2 pi).

[F6]

Rank–nullity holds for a linear map with finite-dimensional domain (Rank-nullity: $\dim_{F} V = nullity T + rank T$ ).

[F7]

A square matrix is invertible exactly when its multiplication map is a linear isomorphism; equivalently, unique solvability for every right-hand side implies invertibility (Invertible matrix theorem: invertibility, full pivot rank, RREF $I$ , trivial nullspace and unique solvability are equivalent).

Proof

1.1givenF3F5algebra

For $1 \leq r < 2 N$ , put $θ = π r / N$ and $S_{r} = \sum_{j = 0}^{N - 1} \cos (θ (j + 1 / 2))$ . Since $0 < θ < 2 π$ , [F5] gives $e^{i θ} \neq 1$ ; the finite geometric identity with $z = e^{i θ}$ therefore gives the complex sum $e^{i θ / 2} (1 - e^{i N θ}) / (1 - e^{i θ})$ . If $r$ is even, its numerator is zero. If $r$ is odd, $e^{i N θ} = - 1$ and $1 - e^{i θ} = - 2 i e^{i θ / 2} \sin (θ / 2)$ by [F3]. Here $0 < θ / 2 < π$ , so [F5] also implies $\sin (θ / 2) \neq 0$ ; the quotient is purely imaginary. In either case its real part is zero, hence $S_{r} = 0$ ; also $S_{0} = N$ .

2.1F1F4step 1.1algebra

Applying [F4] to each product $C_{k j} C_{ℓ j}$ gives $(C C^{T})_{k ℓ} = \frac{1}{2} (S_{∣ k - ℓ ∣} + S_{k + ℓ})$ . If $k = ℓ = 0$ , this is $N$ . If $k = ℓ > 0$ , then $S_{0} = N$ and $0 < 2 k < 2 N$ , so it is $N / 2$ . If $k \neq ℓ$ , then $0 < ∣ k - ℓ ∣ < N$ and $0 < k + ℓ < 2 N$ , so both sums vanish by step 1.1. Thus $C C^{T} = diag (N, N / 2, \dots, N / 2)$ .

3.1step 2.1F6F7algebra

Let $D = diag (N, N / 2, \dots, N / 2)$ . Since $C C^{T} = D$ and $D$ is invertible, $R = C^{T} D^{- 1}$ satisfies $C R = I$ . Hence every $y \in R^{N}$ has a preimage $R y$ under $C$ , so the linear map defined by $C$ is surjective. Its rank is $N$ , and [F6] gives nullity zero; therefore it is bijective and [F7] makes $C$ invertible. The identity $C R = I$ then gives $R = C^{- 1}$ , whose $(j, k)$ entry yields the stated inverse formula.

4.1F1F2F3algebra∎

For each $j$ , the two terms of the DFT of $E_{N} x$ at positions $j$ and $2 N - 1 - j$ are $x_{j} e^{- π i j k / N}$ and $x_{j} e^{π i (j + 1) k / N}$ , since $e^{- 2 π i k} = 1$ . Their sum is $2 x_{j} e^{π i k / (2 N)} \cos (π (j + 1 / 2) k / N)$ by [F3]. Summing over $j$ and applying the normalization in [F2] gives $(F_{2 N} E_{N} x) (k) = 2 (2 N)^{- 1 / 2} e^{π i k / (2 N)} X_{k} = \sqrt{2 / N} e^{π i k / (2 N)} X_{k}$ .

TheoremStatement: Literature-sourcedProof: AI-generatedjudge pass (gpt-6-luna)Open item page →

DST-I is orthogonal with the stated endpoint normalization and odd extension

Statement

Let $S_{k j} = \sin (π j k / (N + 1))$ for $1 \leq j, k \leq N$ . Then $S^{T} S = ((N + 1) / 2) I_{N}$ , so the inverse DST-I is $x_{j} = 2 (N + 1)^{- 1} \sum_{k = 1}^{N} Y_{k} \sin (π j k / (N + 1))$ . For the odd extension $O_{N} x$ of DCT-II and DST-I with fixed even and odd extension conventions, $(F_{2 (N + 1)} O_{N} x) (k) = - 2 i (2 (N + 1))^{- 1 / 2} Y_{k}$ for $1 \leq k \leq N$ .

Facts & Assumptions

Given: An integer $N \geq 1$ , a vector $x \in R^{N}$ , the DST-I matrix $S$ , and the odd extension $O_{N} x$ .

[F1]

$Y_{k} = \sum_{j = 1}^{N} x_{j} \sin (π j k / (N + 1))$ ; $(O_{N} x)_{0} = (O_{N} x)_{N + 1} = 0$ , $(O_{N} x)_{j} = x_{j}$ , and $(O_{N} x)_{2 (N + 1) - j} = - x_{j}$ (DCT-II and DST-I with fixed even and odd extension conventions).

[F2]

$(F_{2 (N + 1)} z) (k) = (2 (N + 1))^{- 1 / 2} \sum_{m = 0}^{2 (N + 1) - 1} z_{m} e^{- 2 π i m k / (2 (N + 1))}$ (The unitary discrete Fourier transform on $Z / N Z$ ).

[F3]

Euler's identity, the complex exponential addition law, and field laws hold ( $\exp (x + i y) = e^{x} (\cos y + i \sin y)$ , $∣ \exp (x + i y) ∣ = e^{x}$ , and $e^{i π} + 1 = 0$ , $\exp (z + w) = \exp z \exp w$ , and the complex exponential extends the real exponential, $C = R [x] / (x^{2} + 1)$ is a field, every element is uniquely $a + b i$ , and every nonzero element has inverse $(a - b i) / (a^{2} + b^{2})$ ).

[F4]

For real $u, v$ , $\sin u \sin v = \frac{1}{2} (\cos (u - v) - \cos (u + v))$ (Product-to-sum and sum-to-product identities).

[F5]

$e^{z} = 1$ exactly when $z \in 2 π i Z$ ( $\ker (\exp) = 2 π i Z$ , and $\exp z = \exp w$ exactly when $z - w \in 2 π i Z$ ).

Proof

1.1givenF3F5algebra

For $1 \leq r \leq 2 N$ , let $T_{r} = \sum_{j = 1}^{N} \cos (π r j / (N + 1))$ . Put $z = e^{π i r / (N + 1)}$ . Since $0 < π r / (N + 1) < 2 π$ , [F5] gives $z \neq 1$ ; also $z^{2 (N + 1)} = 1$ , so the geometric sum over $j = 0, \dots, 2 (N + 1) - 1$ is zero. Pair the terms with indices $j$ and $2 (N + 1) - j$ for $1 \leq j \leq N$ ; they are conjugates, while the unpaired terms are $1$ and $(- 1)^{r}$ . Taking real parts gives $1 + (- 1)^{r} + 2 T_{r} = 0$ . Hence $T_{r} = - 1$ for even $r$ and $T_{r} = 0$ for odd $r$ ; also define $T_{0} = N$ .

2.1F1F4step 1.1algebra

Applying [F4] to each product $S_{k j} S_{ℓ j}$ gives $(S^{T} S)_{k ℓ} = \frac{1}{2} (T_{∣ k - ℓ ∣} - T_{k + ℓ})$ . If $k = ℓ$ , then the first term is $T_{0} = N$ and the second is $T_{2 k} = - 1$ , giving $(N + 1) / 2$ . If $k \neq ℓ$ , then $∣ k - ℓ ∣$ and $k + ℓ$ have the same parity and both are between $1$ and $2 N$ ; their $T$ values from step 1.1 agree, giving zero. Therefore $S^{T} S = ((N + 1) / 2) I_{N}$ .

3.1F1F2F3step 2.1algebra∎

Since $S$ is symmetric and $(N + 1) / 2 > 0$ , step 2.1 gives $S^{- 1} = 2 (N + 1)^{- 1} S$ , proving the inverse formula. For the odd extension, pair its positions $j$ and $2 (N + 1) - j$ : their DFT terms are $x_{j} e^{- π i j k / (N + 1)} - x_{j} e^{π i j k / (N + 1)} = - 2 i x_{j} \sin (π j k / (N + 1))$ by [F3]. Summing and applying [F2] gives the claimed DFT relation.

TheoremStatement: Literature-sourcedProof: AI-generatedjudge pass (gpt-6-luna)Open item page →

DCT-II samples extract Chebyshev coefficients at first-kind roots

Statement

Let $p (t) = \sum_{k = 0}^{N - 1} a_{k} T_{k} (t)$ be a polynomial of degree less than $N$ , where $T_{k}$ is the Chebyshev polynomial of the first kind, and sample $x_{j} = p (\cos (π (j + 1 / 2) / N))$ for $0 \leq j < N$ . If $X$ is the DCT-II of $x$ , then $a_{0} = X_{0} / N$ and $a_{k} = 2 X_{k} / N$ for $1 \leq k < N$ .

Facts & Assumptions

Given: $N \geq 1$ , coefficients $a_{0}, \dots, a_{N - 1}$ , Chebyshev polynomials $T_{k}$ , and the DCT-II convention of DCT-II and DST-I with fixed even and odd extension conventions.

[F1]

$T_{k} (\cos θ) = \cos (k θ)$ for every integer $k \geq 0$ , by the recurrence definition together with the addition formulas (Chebyshev polynomials of the first and second kinds by their three-term recurrences, $T_{n} (\cos θ) = \cos (n θ)$ and $U_{n} (\cos θ) \sin θ = \sin ((n + 1) θ)$ for every $n \in N$ ).

[F2]

The DCT-II matrix is $C_{k j} = \cos (π (j + 1 / 2) k / N)$ , and its row-orthogonality factor is $C C^{T} = diag (N, N / 2, \dots, N / 2)$ (DCT-II and DST-I with fixed even and odd extension conventions, DCT-II is orthogonal with the stated endpoint normalization and even extension).

Proof

1.1givenF1F2

Put $θ_{j} = π (j + 1 / 2) / N$ . By [F1], $x_{j} = \sum_{k = 0}^{N - 1} a_{k} \cos (k θ_{j}) = \sum_{k = 0}^{N - 1} C_{k j} a_{k}$ , so in vector form $x = C^{T} a$ .

2.1step 1.1F2algebra∎

The DCT output is $X = C x = C C^{T} a$ . By [F2], $C C^{T}$ is diagonal with entries $N$ at index zero and $N / 2$ at positive indices. Thus $X_{0} = N a_{0}$ and $X_{k} = (N / 2) a_{k}$ for $k \geq 1$ . Dividing by these positive constants gives the claimed coefficients.

TheoremStatement: Literature-sourcedProof: AI-adaptedOpen item page →

The interval evaluator's arctangent definition uses the canonical Fourier pi

Statement

The value $4 \arctan (1)$ stipulated for $π$ by Real elementary functions and their domains equals the canonical circle constant $π$ used in the finite Fourier conventions, including Pi as twice the smallest positive zero of cosine. Thus the Machin enclosure of $π$ used by the interval evaluator encloses the same constant as the DFT and Bluestein formulas.

Facts & Assumptions

Given: The interval-evaluator definition of $\arctan$ and $π$ , and the published principal inverse tangent.

[F1]

The interval-evaluator definition sets $\arctan x = \int_{0}^{x} (1 + t^{2})^{- 1} d t$ and $π = 4 \arctan (1)$ (Real elementary functions and their domains).

[F2]

The principal arctangent equals the same integral for every real input; for $∣ x ∣ < 1$ it equals its alternating power series, and the endpoint alternating series has sum $π / 4$ (Principal arctangent: derivative, integral, power series, and the Gregory–Leibniz series).

[F3]

If an ordinary real series converges to $s$ , its power series tends to $s$ as its argument increases to $1$ (Abel's limit theorem: if a real series converges to $s$ , then its power series tends to $s$ as $x ↑ 1$ ).

[F4]

A function differentiable at a limit point is continuous there (A function differentiable at $c$ is continuous at $c$ ).

[F5]

The canonical Fourier constant is the positive constant $π = 2 γ$ where $γ$ is the smallest positive zero of cosine (Pi as twice the smallest positive zero of cosine).

[F6]

For the integral-defined arctangent $A$ and $π_{A} = 4 A (1)$ , $16 A (1 / 5) - 4 A (1 / 239) = π_{A}$ (Trigonometric identities, Machin's formula and reduced Taylor tails, Proof step 3.1).

Proof

1.1givenF1F2

By [F1] and the integral identity in [F2], the interval-evaluator arctangent and the principal arctangent have the same value at $1$ .

1.2F2F3F5algebra

Let $a_{n} = (- 1)^{n} / (2 n + 1)$ . The endpoint statement in [F2] says $\sum_{n \geq 0} a_{n} = π / 4$ ; [F3] gives $\lim_{x ↑ 1} \sum_{n \geq 0} a_{n} (x^{2})^{n} = π / 4$ . Multiplication by $x$ and the interior power-series identity in [F2] show $\lim_{x ↑ 1} \arctan x = π / 4$ , where [F5] fixes this as the canonical Fourier constant.

2.1F1F2F4F6step 1.1step 1.2algebra∎

The principal arctangent is differentiable at $1$ by [F2], hence continuous there by [F4]. Therefore step 1.2 gives $\arctan (1) = π / 4$ . Step 1.1 and [F1] now give $4 \arctan (1) = π$ , proving that the interval evaluator and the Fourier formulas use the same constant. Finally [F6] identifies $16 \arctan (1 / 5) - 4 \arctan (1 / 239)$ with $4 \arctan (1)$ , so its Machin enclosure contains this canonical $π$ .

AlgorithmAlgorithm: AI-adaptedProof: Not applicableOpen item page →

Certified dyadic enclosures for unit-circle phases

Definition

Input an exact rational $q$ , an exact rational target $0 < τ \leq 1 / 4$ , integer dyadic exponent limits $e_{\min} \leq e_{\max}$ , an initial precision $p_{0} \geq 2$ , a finite precision cap $p_{\max} \geq p_{0}$ , and a finite elementary-series cap $J \geq 1$ . The target value is $w = e^{π i q} = \cos (π q) + i \sin (π q)$ by $\exp (x + i y) = e^{x} (\cos y + i \sin y)$ , $∣ \exp (x + i y) ∣ = e^{x}$ , and $e^{i π} + 1 = 0$ , for the canonical Fourier $π$ identified with the interval evaluator's constant by The interval evaluator's arctangent definition uses the canonical Fourier pi. For each precision $p = p_{0}, p_{0} + 1, \dots, p_{\max}$ , and $t \in {1 / 5, 1 / 239}$ , form the exact rational interval $I_{J} (t) = [A_{J} (t) - R_{J} (t), A_{J} (t) + R_{J} (t)]$ , where $A_{J} (t) = \sum_{j = 0}^{J} \frac{(- 1)^{j} t^{2 j + 1}}{2 j + 1}, R_{J} (t) = \frac{∣ t ∣^{2 J + 3}}{(2 J + 3) (1 - t^{2})} .$ The Machin identity gives the interval $P = 16 I_{J} (1 / 5) - 4 I_{J} (1 / 239)$ containing $π$ . Form the exact rational interval $q P$ by taking the minimum and maximum of its two endpoint products, then round its endpoints down and up to obtain a dyadic input interval $Θ$ in format $(p, e_{\min}, e_{\max})$ . If this conversion overflows, treat that precision as a non-finite evaluation. Otherwise build two finite real expression DAGs with input interval $Θ$ and requested outputs $\cos (Θ)$ and $\sin (Θ)$ , and call Certified range evaluation for elementary real and complex functions using the same dyadic format and truncation index $J$ ; its fixed finite-cell rules govern range reduction. If either result is non-finite, continue at the next precision unless the precision cap is reached, in which case return inconclusive with the supplier status. If the results are finite intervals $C = [c_{-}, c_{+}]$ and $S = [s_{-}, s_{+}]$ , set $\hat{w} = c_{-} + i s_{-}$ and compute an upward-rounded dyadic $ϵ \geq (c_{+} - c_{-}) + (s_{+} - s_{-})$ . If $ϵ \leq τ$ , return certified(\widehat w,\epsilon); otherwise continue or return inconclusive at the cap. A non-finite endpoint or overflow returns its explicit failure status. The routine makes no finite-error claim on a failure branch.

Remarks

The lower-left rectangle corner is already a representable dyadic pair. The radius uses the sum of coordinate widths, which bounds complex modulus error without a square-root computation. Precision refinement is finite and may return inconclusive; only the certified branch is used by dependent algorithms.

AlgorithmAlgorithm: AI-adaptedProof: Not applicableOpen item page →

Computed unitary radix-two FFT with certified twiddle phases

Definition

Fix an instance of the operationwise floating-point model in Operationwise real and complex floating point arithmetic with rounded square roots and fused multiply-add, with representable real set $F$ . Input $N = 2^{m}$ with $m \geq 0$ , $x \in C^{N}$ whose real and imaginary coordinates all belong to $F$ , a direction forward or inverse, an exact rational twiddle tolerance $0 < τ \leq 1 / 4$ , integer exponent limits $e_{\min} \leq e_{\max}$ for the dyadic phase format, an initial phase precision $p_{0} \geq 2$ , a finite precision cap $p_{\max} \geq p_{0}$ , a finite elementary-series cap $J \geq 1$ , and a maximum radix-two recursion depth $d_{\max} \geq 0$ . At precision $p$ , phase endpoints use the dyadic format with precision $p$ and exponent limits $e_{\min}, e_{\max}$ . Output a vector in $C^{N}$ and the certificate for each twiddle used, or the first range-failure/inconclusive status and its stage, index, or phase. Before multiplying by a certified phase ${\hat{w}}_{k}$ , require its real and imaginary coordinates to belong to $F$ ; otherwise return range-failure at that twiddle index before performing the product. Use the same floating-point model instance throughout the recursion.

For direction forward, at $N = 1$ return $x_{0}$ . For $N > 1$ , if $d_{\max} = 0$ return inconclusive with the subdivision-cap-exhausted status. Otherwise split $x$ into even and odd subsequences and recursively compute their unitary transforms $E, O$ by calling this routine with the same phase parameters and depth cap $d_{\max} - 1$ . For each $0 \leq k < N / 2$ , call Certified dyadic enclosures for unit-circle phases on $(q, τ, e_{\min}, e_{\max}, p_{0}, p_{\max}, J)$ with $q = - 2 k / N$ , obtaining $({\hat{w}}_{k}, ϵ_{k})$ approximating $w_{k} = e^{- 2 π i k / N}$ . Compute the six-real-operation complex product $t_{k} = {fl}_{C} ({\hat{w}}_{k} O_{k})$ , componentwise rounded sums $s_{k}^{+} = {fl}_{C} (E_{k} + t_{k})$ and $s_{k}^{-} = {fl}_{C} (E_{k} - t_{k})$ , and $\hat{c} = fl (\sqrt{1 / 2})$ . Set output entries $k$ and $k + N / 2$ to the componentwise rounded real scalings ${fl}_{C} (\hat{c} s_{k}^{+})$ and ${fl}_{C} (\hat{c} s_{k}^{-})$ , respectively. For direction inverse, conjugate the input, run the forward routine with the same phase parameters and $d_{\max}$ , and conjugate its output; conjugation is exact in the operationwise model. Propagate any failed phase certificate; on arithmetic failure return the operationwise range status and index. Accuracy is asserted only by the separate theorem and only on a successful no-range-failure execution.

TheoremStatement: AI-adaptedProof: AI-generatedOpen item page →

Certified phase output encloses the exact unit-circle value

Statement

Whenever Certified dyadic enclosures for unit-circle phases returns certified(\widehat w,\epsilon) on rational input $q$ , its output satisfies $∣ \hat{w} - e^{π i q} ∣ \leq ϵ \leq τ$ .

Facts & Assumptions

Given: Rational $q$ , target $0 < τ \leq 1 / 4$ , finite precision and evaluation budgets, and a successful output certified(\widehat w,\epsilon).

[F1]

On a finite return, the real interval $C$ contains $\cos (π q)$ and $S$ contains $\sin (π q)$ , and Euler's identity identifies their pair with $e^{π i q}$ ; the interval and Fourier definitions use the same $π$ (Certified real elementary-function evaluations enclose their ranges, Real intervals, complex rectangles and ball enclosures, $\exp (x + i y) = e^{x} (\cos y + i \sin y)$ , $∣ \exp (x + i y) ∣ = e^{x}$ , and $e^{i π} + 1 = 0$ , The interval evaluator's arctangent definition uses the canonical Fourier pi).

[F2]

The returned lower endpoints are the lower endpoints of $C, S$ , and $ϵ$ is an upward-rounded dyadic at least $(c_{+} - c_{-}) + (s_{+} - s_{-})$ ; directed interval arithmetic encloses these exact widths (Certified dyadic enclosures for unit-circle phases, Outward-rounded operations are inclusion-isotone and enclose DAG evaluations).

[F3]

For real $a, b$ , $∣ a + i b ∣ \leq ∣ a ∣ + ∣ b ∣$ by the complex modulus triangle inequality (Conjugation is an involutive real-field automorphism, $z \overline{z} = ∣ z ∣^{2}$ , and modulus is definite, multiplicative, and subadditive).

Proof

1.1givenF1F2

By [F1], the real coordinate of $w = e^{π i q}$ lies in $[c_{-}, c_{+}]$ and its imaginary coordinate lies in $[s_{-}, s_{+}]$ . Since $\hat{w} = c_{-} + i s_{-}$ by [F2], the coordinate differences satisfy $∣ ℜ (w - \hat{w}) ∣ \leq c_{+} - c_{-}$ and $∣ ℑ (w - \hat{w}) ∣ \leq s_{+} - s_{-}$ .

2.1step 1.1F2F3algebra∎

Applying [F3] to these two coordinate differences gives $∣ w - \hat{w} ∣ \leq (c_{+} - c_{-}) + (s_{+} - s_{-}) \leq ϵ$ . The algorithm returns only when $ϵ \leq τ$ , proving the claim.

AlgorithmAlgorithm: Literature-sourcedProof: Not applicableOpen item page →

Computed FFT cyclic convolution and zero-padded linear convolution

Definition

Fix an instance of the operationwise floating-point model in Operationwise real and complex floating point arithmetic with rounded square roots and fused multiply-add, with representable real set $F$ . Input two finite complex sequences $a, b$ whose real and imaginary coordinates all belong to $F$ , a mode cyclic or linear, and phase/FFT controls: an exact rational twiddle tolerance $0 < τ \leq 1 / 4$ , integer dyadic exponent limits $e_{\min} \leq e_{\max}$ , an initial phase precision $p_{0} \geq 2$ , a finite precision cap $p_{\max} \geq p_{0}$ , a finite elementary-series cap $J \geq 1$ , and a maximum radix-two recursion depth $d_{\max} \geq 0$ . At precision $p$ , phase endpoints use the dyadic format with precision $p$ and exponent limits $e_{\min}, e_{\max}$ . In cyclic mode the inputs must have the same length $M = 2^{q}$ with integer $q \geq 0$ ; the requested exact operation is cyclic convolution on $Z / M Z$ . In linear mode, with lengths $n, m \geq 1$ , set $L = n + m - 1$ , choose the least power of two $M \geq L$ , and append zeros to both inputs to length $M$ ; the requested output is the first $L$ entries of their cyclic convolution. Validate the mode, positive dimensions, input representability, exact rational tolerance, integer control parameters and their displayed inequalities, and the cyclic power-of-two condition before any transform call. If any such input constraint fails, return inconclusive with diagnostic invalid-input and the offending parameter or coordinate, without a vector output. An arithmetic range failure returns range-failure with its operation and index, without a vector output.

For either mode, call Computed unitary radix-two FFT with certified twiddle phases on each length- $M$ input with direction forward and the same $(τ, e_{\min}, e_{\max}, p_{0}, p_{\max}, J, d_{\max})$ controls. Multiply corresponding transform entries as ${fl}_{C} ({\hat{A}}_{k} {\hat{B}}_{k})$ using the six-real-operation complex product of Operationwise real and complex floating point arithmetic with rounded square roots and fused multiply-add. Call the same FFT routine on this product vector with direction inverse and the same controls, and multiply each output component by the computed real scalar $fl (\sqrt{M})$ using componentwise rounded real scaling. Return the length- $M$ cyclic result, or its first $L$ entries in linear mode, only if every phase, recursion-budget, and arithmetic operation succeeds; otherwise return the first inconclusive or range-failure status and index. The exact convolution identity and numerical-error claim are separate results.

LemmaStatement: AI-adaptedProof: AI-generatedOpen item page →

A computed unitary radix-two butterfly has error at most (10u + εtw) times its input norm

Statement

Let $u \leq 1 / 16$ , let $∣ w ∣ = 1$ , and let $∣ \hat{w} - w ∣ \leq ϵ \leq 1 / 4$ . For $a, b \in C$ , perform the successful butterfly of Computed unitary radix-two FFT with certified twiddle phases using $\hat{w}$ . If each exact nonzero real intermediate is normal and every operation succeeds under the operationwise model, then its computed pair $({\hat{y}}_{+}, {\hat{y}}_{-})$ satisfies

{∥ ({\hat{y}}_{+}, {\hat{y}}_{-}) - 2^{- 1 / 2} (a + w b, a - w b) ∥}_{2} \leq (10 u + ϵ) {∥ (a, b) ∥}_{2} .

Facts & Assumptions

Given: $a, b, w, \hat{w}, u, ϵ$ with the stated bounds, and a successful rounded butterfly.

[F5]

For a pair $(z_{1}, z_{2}) \in C^{2}$ , $∥ (z_{1}, z_{2}) ∥_{2} = (∣ z_{1} ∣^{2} + ∣ z_{2} ∣^{2})^{1 / 2}$ (The counting inner product on $C^{Z / N Z}$ ).

[F1]

Six-operation complex multiplication obeys $∣ {fl}_{C} (z w) - z w ∣ \leq 3 u ∣ z ∣ ∣ w ∣$ when its six operations succeed (Six-operation complex multiplication has a normwise roundoff bound).

[F2]

$∣ w ∣ = 1$ and $∣ \hat{w} - w ∣ \leq ϵ$ ; hence $∣ \hat{w} ∣ \leq 1 + ϵ$ ( $\exp (x + i y) = e^{x} (\cos y + i \sin y)$ , $∣ \exp (x + i y) ∣ = e^{x}$ , and $e^{i π} + 1 = 0$ , Certified phase output encloses the exact unit-circle value, Conjugation is an involutive real-field automorphism, $z \overline{z} = ∣ z ∣^{2}$ , and modulus is definite, multiplicative, and subadditive).

[F3]

Each real rounded addition, multiplication, or positive in-range square root has relative error at most $u$ ; complex addition and real scaling round their two components separately (The standard relative floating-point model and unit roundoff, Operationwise real and complex floating point arithmetic with rounded square roots and fused multiply-add).

[F4]

The complex modulus is subadditive; Cauchy–Schwarz bounds the Euclidean norm of two coordinatewise error bounds (Conjugation is an involutive real-field automorphism, $z \overline{z} = ∣ z ∣^{2}$ , and modulus is definite, multiplicative, and subadditive, Cauchy–Schwarz: $∣ ⟨ u, v ⟩ ∣ \leq ∥ u ∥ ∥ v ∥$ , with equality exactly for linearly dependent vectors).

Proof

1.1givenF1F2algebra

Let $t = {fl}_{C} (\hat{w} b)$ . By [F1] and [F2], $∣ t - w b ∣ \leq 3 u ∣ \hat{w} ∣ ∣ b ∣ + ϵ ∣ b ∣ \leq (15 u / 4 + ϵ) ∣ b ∣$ . Thus $∣ t ∣ \leq (1 + 15 u / 4 + ϵ) ∣ b ∣ \leq 95 ∣ b ∣ / 64$ .

1.2givenF3F4algebra

For $q_{\pm} = {fl}_{C} (a \pm t)$ , [F3] gives $∣ q_{\pm} - (a \pm t) ∣ \leq u (∣ a ∣ + ∣ t ∣)$ by [F4]. If $c = 1 / \sqrt{2}$ and $\hat{c} = fl (\sqrt{1 / 2})$ , then correct rounding gives $∣ \hat{c} - c ∣ \leq u c$ . The componentwise rounded scaling therefore satisfies $∣ {fl}_{C} (\hat{c} q_{\pm}) - c q_{\pm} ∣ \leq (2 u + u^{2}) c ∣ q_{\pm} ∣$ , and $∣ q_{\pm} ∣ \leq (1 + u) (∣ a ∣ + ∣ t ∣)$ .

2.1step 1.1step 1.2algebra

Combining steps 1.1–1.2 and using $\sqrt{2} c = 1$ , the two-output error norm is at most $(15 u / 4 + ϵ) ∣ b ∣ + α (∣ a ∣ + 95 ∣ b ∣ / 64)$ , where $α = u + (2 u + u^{2}) (1 + u) = 3 u + 3 u^{2} + u^{3} \leq 817 u / 256$ . Since $15 / 4 + (817 / 256) (95 / 64) < 9$ , this is at most $(817 u / 256) ∣ a ∣ + (9 u + ϵ) ∣ b ∣$ .

3.1step 2.1F4F5algebra∎

By [F4], $\sqrt{(817 / 256)^{2} + 9^{2}} < 10$ , so $(817 u / 256) ∣ a ∣ + 9 u ∣ b ∣ \leq 10 u \sqrt{∣ a ∣^{2} + ∣ b ∣^{2}}$ ; [F5] identifies the square-root factor as $∥ (a, b) ∥_{2}$ . Also $ϵ ∣ b ∣ \leq ϵ \sqrt{∣ a ∣^{2} + ∣ b ∣^{2}}$ . Adding the bounds proves the claim.

AlgorithmAlgorithm: Literature-sourcedProof: Not applicableOpen item page →

Arbitrary-length mixed-radix and Bluestein DFT

Definition

Fix an instance of the operationwise floating-point model in Operationwise real and complex floating point arithmetic with rounded square roots and fused multiply-add, with representable real set $F$ . Input $N \geq 1$ , $x \in C^{N}$ whose real and imaginary coordinates all belong to $F$ , an exact rational phase tolerance $0 < τ \leq 1 / 4$ , integer dyadic exponent limits $e_{\min} \leq e_{\max}$ , an initial phase precision $p_{0} \geq 2$ , a finite precision cap $p_{\max} \geq p_{0}$ , a finite elementary-series cap $J \geq 1$ , and a maximum radix-two recursion depth $d_{\max} \geq 0$ . At precision $p$ , phase endpoints use the dyadic format with precision $p$ and exponent limits $e_{\min}, e_{\max}$ . Before any certified phase is used as an arithmetic operand or stored in an input vector to a computed subroutine, require its real and imaginary coordinates to belong to $F$ ; otherwise return range-failure with that phase index before the operation or subcall. Use the same floating-point model instance throughout all recursive calls. The output convention is the unnormalised DFT $X_{k} = \sum_{j = 0}^{N - 1} x_{j} e^{- 2 π i j k / N}$ (The unnormalised engineering DFT and its conversion to the unitary transform); the procedure returns an approximate vector or the first inconclusive/range-failure status.

For $N = 1$ , return $x_{0}$ .
If $N = 2^{m}$ with $m \geq 1$ , call Computed unitary radix-two FFT with certified twiddle phases in forward mode with controls $(τ, e_{\min}, e_{\max}, p_{0}, p_{\max}, J, d_{\max})$ and multiply each output by the computed scalar $fl (\sqrt{N})$ . A failed subcall or scale operation returns its status.
If $N$ is composite and not a power of two, find its least divisor $r$ with $1 < r < N$ by trial division and put $s = N / r$ . For $u = 0, \dots, r - 1$ , recursively call this procedure with the same controls on the strided vector $(x_{u + r v})_{v = 0}^{s - 1}$ to obtain an approximate unnormalised length- $s$ DFT $U_{u}$ . For each $0 \leq k_{2} < s$ and $0 \leq u < r$ , call Certified dyadic enclosures for unit-circle phases on $(q, τ, e_{\min}, e_{\max}, p_{0}, p_{\max}, J)$ with $q = - 2 u k_{2} / N$ to obtain the twiddle enclosure, then multiply its returned phase by $U_{u} (k_{2})$ using ${fl}_{C}$ . For each $0 \leq k_{2} < s$ , apply this procedure recursively, with the same controls, to the length- $r$ vector of these products indexed by $u$ , and for every $0 \leq k_{1} < r$ store its output $k_{1}$ at index $k_{2} + s k_{1}$ . Propagate any failed phase or arithmetic status.
If $N$ is prime, set $M$ to the least power of two with $M \geq 2 N - 1$ . For each $0 \leq j < N$ , call Certified dyadic enclosures for unit-circle phases on $(q, τ, e_{\min}, e_{\max}, p_{0}, p_{\max}, J)$ with $q = - j^{2} / N$ to obtain ${\hat{α}}_{j}$ , and set $a_{j} = {fl}_{C} ({\hat{α}}_{j} x_{j})$ ; set the remaining entries of $a \in C^{M}$ to zero. For each $- (N - 1) \leq t \leq N - 1$ , obtain ${\hat{β}}_{t}$ by the same call with $q = t^{2} / N$ and set $b_{[t]_{M}} = {\hat{β}}_{t}$ , setting other entries of $b$ to zero. Compute the cyclic convolution with Computed FFT cyclic convolution and zero-padded linear convolution in cyclic mode using controls $(τ, e_{\min}, e_{\max}, p_{0}, p_{\max}, J, d_{\max})$ . For each $0 \leq k < N$ , obtain ${\hat{γ}}_{k}$ by the same phase call with $q = - k^{2} / N$ and set ${\hat{X}}_{k} = {fl}_{C} ({\hat{γ}}_{k} {\hat{c}}_{k})$ . Return the first $N$ entries or the first failure status.

The divisor search and phase-enclosure work have separate integer/interval budgets and are not counted as complex arithmetic operations in the operation-count theorem. The exact factorization lemmas and the chirp-phase error theorem state what is guaranteed; this procedure by itself makes no general mixed-radix stability claim.

TheoremStatement: Literature-sourcedProof: AI-generatedOpen item page →

Computed unitary radix-two FFT has the multiplicative stage error bound

Statement

Let $N = 2^{m}$ , $m \geq 0$ , and suppose a successful execution of Computed unitary radix-two FFT with certified twiddle phases has no range failures. If every twiddle certificate satisfies $ϵ_{ℓ} \leq ϵ_{t w} \leq 1 / 4$ and $u \leq 1 / 16$ , put $η = 10 u + ϵ_{t w}$ . Then

∥ {\hat{F}}_{N} x - F_{N} x ∥_{2} \leq ((1 + η)^{m} - 1) ∥ x ∥_{2} .

The conjugated inverse routine has the same bound. If $m η \leq 1$ , the right side is at most $2 m η ∥ x ∥_{2}$ . Thus an $O (u \log N)$ shorthand requires both $ϵ_{t w} = O (u)$ and $u \log N$ small; the theorem makes no such shorthand when twiddle errors are uncontrolled.

Facts & Assumptions

Given: $N = 2^{m}$ , $x \in C^{N}$ , a successful execution, $u \leq 1 / 16$ , and $ϵ_{ℓ} \leq ϵ_{t w} \leq 1 / 4$ for every twiddle certificate.

[F1]

For $N = 2 M$ , split $x$ into even and odd subsequences $e, o$ ; the unnormalised DFT obeys $X_{k} (x) = X_{k} (e) + e^{- 2 π i k / N} X_{k} (o)$ and $X_{k + M} (x) = X_{k} (e) - e^{- 2 π i k / N} X_{k} (o)$ , while $F_{N} = N^{- 1 / 2} X_{N}$ (The radix-two even/odd factorisation of the DFT, The unnormalised engineering DFT and its conversion to the unitary transform, The unitary discrete Fourier transform on $Z / N Z$ ). The recursive FFT computes $X_{N}$ exactly (The recursive radix-two fast Fourier transform, Correctness of the recursive radix-two FFT).

[F2]

Every computed normalized butterfly differs from its exact block by at most $η$ times the pair input norm (A computed unitary radix-two butterfly has error at most (10u + εtw) times its input norm).

[F3]

The unitary DFT preserves the squared Euclidean norm; in particular the exact transforms of the even and odd subsequences have squared norms summing to $∥ x ∥_{2}^{2}$ (Finite Parseval and Plancherel identity for the unitary DFT).

[F4]

Complex conjugation preserves modulus and the Euclidean norm, and the inverse unitary DFT is obtained by conjugating input, applying the forward DFT, and conjugating output (Conjugation is an involutive real-field automorphism, $z \overline{z} = ∣ z ∣^{2}$ , and modulus is definite, multiplicative, and subadditive, Finite Fourier inversion for the unitary transform on $Z / N Z$ ).

[F5]

The real exponential is its power series (Real elementary functions and their domains).

[F6]

For $z \in C^{N}$ , $∥ z ∥_{2}^{2} = \sum_{j} ∣ z_{j} ∣^{2}$ (The counting inner product on $C^{Z / N Z}$ ).

[F7]

Every exact twiddle $w = e^{- 2 π i k / N}$ has $∣ w ∣ = 1$ , and $\sqrt{N} = \sqrt{2} \sqrt{N / 2}$ for $N \geq 2$ ( $\exp (x + i y) = e^{x} (\cos y + i \sin y)$ , $∣ \exp (x + i y) ∣ = e^{x}$ , and $e^{i π} + 1 = 0$ , Rational powers $a^{r}$ of a positive base, Laws of rational exponents).

Proof

1.1givenF1F7algebrabase

For $m = 0$ , $N = 1$ and the algorithm returns the input exactly, so the error is $0 = (1 + η)^{0} - 1$ . For $m \geq 1$ , [F1] gives the exact unnormalised top outputs $X_{k} (e) + w X_{k} (o)$ and $X_{k} (e) - w X_{k} (o)$ , with $w = e^{- 2 π i k / N}$ . The child unitary outputs are $E_{k} = X_{k} (e) / \sqrt{N / 2}$ and $O_{k} = X_{k} (o) / \sqrt{N / 2}$ ; dividing the parent unnormalised outputs by $\sqrt{N} = \sqrt{2} \sqrt{N / 2}$ gives exactly $U_{w} (E_{k}, O_{k}) = 2^{- 1 / 2} (E_{k} + w O_{k}, E_{k} - w O_{k})$ . By [F7], $∣ w ∣ = 1$ and the scaling identity holds, so the columns of $U_{w}$ have squared norm $1$ and inner product $2^{- 1} (w - w) = 0$ ; hence $U_{w}$ preserves the pair norm.

2.1givenF2F3F6step 1.1algebrainductionih

Suppose the recursive transform of length $N / 2$ has error bound $ρ_{m - 1} = (1 + η)^{m - 1} - 1$ . By [F3] and [F6], stacking the exact even and odd transforms has norm $∥ x ∥_{2}$ ; stacking their computed errors has norm at most $ρ_{m - 1} ∥ x ∥_{2}$ , since the even and odd input coordinates partition $x$ . Hence the stacked computed suboutputs have norm at most $(1 + ρ_{m - 1}) ∥ x ∥_{2}$ . For each $k$ , [F2] bounds the top butterfly error by $η$ times the norm of its computed input pair. Squaring these bounds and summing over $k$ shows that the total top-stage error norm is at most $η (1 + ρ_{m - 1}) ∥ x ∥_{2}$ . The prior subtransform errors contribute at most $ρ_{m - 1} ∥ x ∥_{2}$ through the direct sum of exact norm-preserving blocks in step 1.1. Therefore the total error is at most $[ρ_{m - 1} + η (1 + ρ_{m - 1})] ∥ x ∥_{2} = ((1 + η)^{m} - 1) ∥ x ∥_{2}$ . Induction proves the forward bound for every $m$ .

3.1F4step 2.1

By [F4], the inverse implementation is conjugation before and after the forward routine and conjugation preserves every coordinate modulus. Applying step 2.1 to $\overline{x}$ gives the same error norm for the inverse direction.

4.1givenF5step 2.1algebradischarge-induction: step 2.1∎

If $m η \leq 1$ , then $1 + η \leq e^{η}$ by the power series [F5], so $(1 + η)^{m} \leq e^{m η}$ . For $0 < t \leq 1$ , $e^{t} - 1 = \sum_{j \geq 1} t^{j} / j! \leq t \sum_{j \geq 1} 1 / j! < 2 t$ ; at $t = 0$ , $e^{t} - 1 = 2 t = 0$ . Thus $((1 + η)^{m} - 1) \leq 2 m η$ . If $ϵ_{t w} = O (u)$ , this is $O (m u) = O (u \log_{2} N)$ , with the smallness condition already imposed.

TheoremStatement: Literature-sourcedProof: AI-generatedOpen item page →

Mixed-radix and Bluestein reductions compute every finite DFT in O(N log N) arithmetic

Statement

For this theorem, define the formal exact-constant model as the complete finite branch schedule of Arbitrary-length mixed-radix and Bluestein DFT, with its radix-two subdivision cap sufficient to finish every call, every certified-phase output (including nested calls) replaced by its exact target $e^{π i q}$ , and every arithmetic instruction interpreted exactly. Phase-evaluation and arithmetic-failure statuses are suppressed in this idealized model; it is not an exactness claim about the computed routine, whose phase values are enclosures. For every integer $N \geq 1$ and $x \in C^{N}$ , the formal model returns the unnormalised DFT $X_{k} = \sum_{j = 0}^{N - 1} x_{j} e^{- 2 π i j k / N}$ . Count complex additions, multiplications, and scalar operations at each butterfly, but exclude integer factor search, exact phase generation, interval evaluation, storage, and bit complexity. If $M$ is the least power of two at least $2 N - 1$ in a Bluestein prime branch, that branch costs $O (M \log M)$ ; for $N \geq 2$ , the full arbitrary-length recursion costs $O (N \log N)$ , while $N = 1$ takes constant work.

Facts & Assumptions

Given: $N \geq 1$ , $x \in C^{N}$ , and the formal exact-constant model defined in the Statement, with exact arithmetic instructions.

[F1]

The formal exact-constant model follows the branch schedule, unnormalised output convention, and exact phase targets of the algorithm; powers of two use the unitary radix-two transform scaled by $\sqrt{N}$ , composites use the fixed mixed-radix schedule, and primes use a length- $M$ chirp convolution (Arbitrary-length mixed-radix and Bluestein DFT, The unnormalised engineering DFT and its conversion to the unitary transform).

[F2]

The exact mixed-radix branch equals the DFT by Mixed-radix Cooley–Tukey factorization with fixed index maps, and the exact prime-length branch equals it by Bluestein’s chirp identity embeds an arbitrary DFT in a padded cyclic convolution.

[F3]

The radix-two arithmetic count is at most $2 n \log_{2} n$ for power-of-two length $n$ , with twiddle generation excluded (The radix-two FFT uses $O (N \log_{2} N)$ complex arithmetic operations).

[F4]

In exact arithmetic, the cyclic convolution algorithm takes two unitary DFTs, pointwise multiplication, the inverse unitary DFT, and the $\sqrt{M}$ scale; the published DFT convolution identity and inversion theorem identify this with $a *_{M} b$ (Computed FFT cyclic convolution and zero-padded linear convolution, The DFT turns cyclic convolution into a scaled pointwise product, Finite Fourier inversion for the unitary transform on $Z / N Z$ , The unnormalised cyclic convolution on $Z / N Z$ ).

[F5]

For $N = 2^{m}$ , the exact radix-two recursion computes $X_{N} = \sqrt{N} F_{N} x$ (Correctness of the recursive radix-two FFT, The unitary discrete Fourier transform on $Z / N Z$ , The unnormalised engineering DFT and its conversion to the unitary transform).

[F6]

For every positive $N$ , $X_{N} = \sqrt{N} F_{N} x$ by the two transform normalizations; $\sqrt{r s} = \sqrt{r} \sqrt{s}$ for positive integers (The unitary discrete Fourier transform on $Z / N Z$ , The unnormalised engineering DFT and its conversion to the unitary transform, Rational powers $a^{r}$ of a positive base, Laws of rational exponents).

[F7]

The natural logarithm is strictly increasing, satisfies $\log (x y) = \log x + \log y$ for positive $x, y$ , and $\log_{2} x = \log x / \log 2$ ; since $2 > 1$ , strict increase and $\log 1 = 0$ give $\log 2 > 0$ , so $\log_{2}$ is strictly increasing. Repeated application of the product law gives $\log_{2} (2^{L}) = L$ and $\log_{2} (\prod_{i} n_{i}) = \sum_{i} \log_{2} n_{i}$ for each finite list of positive integers (The natural logarithm as the inverse of the exponential function, Order, continuity, range, and the product, quotient, and reciprocal laws for the natural logarithm, The logarithm to a positive base other than one, The principle of mathematical induction).

[F8]

Strong induction on positive integer lengths follows by choosing the least counterexample using the well-ordering principle; ordinary induction on a finite binary-tree height is the induction principle (The well-ordering principle, The principle of mathematical induction).

Proof

1.1givenF1F5F8algebrabase

We prove correctness by strong induction on $N$ , justified by [F8]. At $N = 1$ , step 1 of the formal model returns $x_{0}$ , the one-term unnormalised DFT. If $N$ is a power of two, its radix-two branch calls the exact unitary transform and scales by $\sqrt{N}$ ; [F5] identifies that result with the unnormalised DFT.

1.2ih

For the induction step, assume the formal model returns the unnormalised DFT for every positive length smaller than $N$ .

1.3F1F3F7algebra

For a power-of-two size $M \geq 2$ , [F3] gives $O (M \log_{2} M)$ ; the output normalization uses one square root and $O (M)$ component scalings, which preserve this order. In a Bluestein branch, the least power of two $M \geq 2 N - 1$ satisfies $M < 4 N$ , and three radix-two transforms plus $O (M)$ pointwise products/scalings cost $O (M \log_{2} M)$ . Since $\log_{2} M \leq \log_{2} (4 N) = 2 + \log_{2} N \leq 3 \log_{2} N$ for $N \geq 2$ , this is $O (N \log_{2} N)$ for a prime input $N$ .

2.1givenstep 1.2F1F2F6inductionalgebra

If $N$ is composite and not a power of two, the least-divisor branch splits $N = r s$ with $r, s < N$ . The hypothesis in step 1.2 makes the recursive unnormalised length- $s$ and length- $r$ transforms exact. By [F6], the length- $s$ output on each strided vector is $\sqrt{s} F_{s} x^{(u)}$ , and a length- $r$ unnormalised transform is $\sqrt{r}$ times the unitary length- $r$ DFT. The fixed strided inputs, exact phase factors, and output indexing are those in [F2], so the algorithm combines the recursive outputs as $X_{k_{2} + s k_{1}} = \sum_{u = 0}^{r - 1} U_{u} (k_{2}) e^{- 2 π i u k_{2} / N} e^{- 2 π i u k_{1} / r}$ . This is the unnormalised mixed-radix formula; hence every output is exact.

2.2givenF1F2F4F6step 1.1algebra

If $N$ is prime, [F2] gives $F_{N} x = N^{- 1 / 2} γ ⊙ c$ for the exact cyclic convolution $c$ in the Bluestein branch. The formal model computes $c$ exactly by [F4]. Multiplying the unitary output identity by $\sqrt{N}$ gives $X_{N} = γ ⊙ c$ , which is exactly the unnormalised output of the branch. The base, power-of-two, composite and prime cases exhaust all positive $N$ .

3.1givenF1F7F8step 1.1step 1.3algebrainductiondischarge-induction: step 1.2∎

At a composite non-power-of-two node $N = r s$ , the algorithm's actual call schedule invokes $r$ transforms of size $s$ and $s$ transforms of size $r$ , plus at most $C N$ arithmetic operations for a fixed $C$ . Thus $T (N) \leq r T (s) + s T (r) + C N$ , so the normalized count $S (N) : = T (N) / N$ satisfies $S (N) \leq S (r) + S (s) + C$ . For the bound, build a separate binary comparison tree by expanding this normalized recurrence at each composite non-power-of-two size; stop at prime or power-of-two sizes. This tree is finite because every child size is smaller than its parent, and its terminal sizes $n_{1}, \dots, n_{L}$ multiply to $N$ because each split replaces $n$ by factors $r, s$ with $n = r s$ . Every terminal size is at least $2$ , so $2^{L} \leq N$ and [F7] gives $L \leq \log_{2} N$ . A full binary tree with $L$ leaves has $L - 1$ internal nodes: a leaf has zero internal nodes, while a root with subtrees having $L_{1}, L_{2}$ leaves has $1 + (L_{1} - 1) + (L_{2} - 1) = L_{1} + L_{2} - 1$ internal nodes, by induction on height [F8]. The product law in [F7] gives $\sum_{i = 1}^{L} \log_{2} n_{i} = \log_{2} N$ . By steps 1.1 and 1.3, each terminal size $n_{i}$ , whether prime or a power of two, contributes $O (\log_{2} n_{i})$ to $S$ , with one uniform constant. Summing the terminal contributions and the at most $C$ contribution at each internal node yields $S (N) = O (\log_{2} N)$ and hence $T (N) = O (N \log_{2} N)$ . The comparison tree bounds the normalized recurrence; it is not the algorithm's actual call tree.

TheoremStatement: Literature-sourcedProof: AI-generatedOpen item page →

Computed FFT convolution error separates transform, product, inverse, and final-scale terms

Statement

Let $M = 2^{q}$ with integer $q \geq 0$ and let $a, b \in C^{M}$ be the arrays actually transformed by Computed FFT cyclic convolution and zero-padded linear convolution (zero-padded in linear mode). Assume all operations succeed, $u \leq 1 / 16$ , the twiddle certificates satisfy $ϵ_{ℓ} \leq ϵ_{t w} \leq 1 / 4$ , and put $η = 10 u + ϵ_{t w}$ , $ρ = (1 + η)^{q} - 1$ , $ϵ_{p} = 3 u$ , $ϵ_{s} = 2 u + u^{2}$ , and $λ = 2 ρ + ρ^{2} + ϵ_{p} (1 + ρ)^{2}$ . If $c = a *_{M} b$ is the exact unnormalised cyclic convolution, then the computed length- $M$ result $\hat{c}$ satisfies

∥ \hat{c} - c ∥_{2} \leq \sqrt{M} [λ + ρ (1 + λ) + ϵ_{s} (1 + ρ) (1 + λ)] ∥ a ∥_{2} ∥ b ∥_{2} .

In linear mode, the same right side bounds the error of the first $n + m - 1$ outputs relative to the exact linear convolution, because zero padding removes cyclic aliasing. This is a conditional error estimate, not an unconditional success certificate.

Facts & Assumptions

Given: The successful computation, length- $M$ input arrays, unitary transforms $A = F_{M} a$ and $B = F_{M} b$ , $u, ϵ_{t w}$ , and the constants in the Statement.

[F1]

The transform and inverse-transform errors satisfy $∥ \hat{A} - A ∥_{2} \leq ρ ∥ a ∥_{2}$ , $∥ \hat{B} - B ∥_{2} \leq ρ ∥ b ∥_{2}$ , and $∥ \hat{z} - F_{M}^{- 1} \hat{p} ∥_{2} \leq ρ ∥ \hat{p} ∥_{2}$ (Computed unitary radix-two FFT has the multiplicative stage error bound, Finite Fourier inversion for the unitary transform on $Z / N Z$ ).

[F2]

For vectors $u, v \in C^{M}$ , $∥ u ⊙ v ∥_{2} \leq ∥ u ∥_{2} ∥ v ∥_{2}$ , since $\sum_{k} ∣ u_{k} ∣^{2} ∣ v_{k} ∣^{2} \leq (\sum_{k} ∣ u_{k} ∣^{2}) (\sum_{k} ∣ v_{k} ∣^{2})$ (The counting inner product on $C^{Z / N Z}$ , Real and imaginary parts, complex conjugation, and modulus).

[F3]

Finite Parseval gives $∥ A ∥_{2} = ∥ a ∥_{2}$ , $∥ B ∥_{2} = ∥ b ∥_{2}$ , and $∥ F_{M}^{- 1} z ∥_{2} = ∥ z ∥_{2}$ ; the cyclic convolution identity is $a *_{M} b = \sqrt{M} F_{M}^{- 1} (A ⊙ B)$ (Finite Parseval and Plancherel identity for the unitary DFT, Finite Fourier inversion for the unitary transform on $Z / N Z$ , The DFT turns cyclic convolution into a scaled pointwise product).

[F4]

Each pointwise complex product error is at most $3 u ∣ {\hat{A}}_{k} ∣ ∣ {\hat{B}}_{k} ∣$ (Six-operation complex multiplication has a normwise roundoff bound).

[F5]

Put $s = fl (\sqrt{M})$ . The model gives $∣ s - \sqrt{M} ∣ \leq u \sqrt{M}$ and each componentwise rounded product with $s$ has error at most $u ∣ s z_{j} ∣$ . Thus the combined scale error is at most $[u (1 + u) + u] \sqrt{M} ∥ z ∥_{2} = (2 u + u^{2}) \sqrt{M} ∥ z ∥_{2} = ϵ_{s} \sqrt{M} ∥ z ∥_{2}$ on a successful normal-range result (Operationwise real and complex floating point arithmetic with rounded square roots and fused multiply-add, The standard relative floating-point model and unit roundoff).

[F6]

The Euclidean norm on $C^{M}$ is subadditive, by expansion of its square and Cauchy–Schwarz (The counting inner product on $C^{Z / N Z}$ , Cauchy–Schwarz: $∣ ⟨ u, v ⟩ ∣ \leq ∥ u ∥ ∥ v ∥$ , with equality exactly for linearly dependent vectors).

Proof

1.1givenF1F2F3F4algebra

Write $\hat{A} = A + d_{A}$ and $\hat{B} = B + d_{B}$ . By [F1], $∥ d_{A} ∥_{2} \leq ρ ∥ a ∥_{2}$ and $∥ d_{B} ∥_{2} \leq ρ ∥ b ∥_{2}$ , while [F3] gives $∥ A ∥_{2} = ∥ a ∥_{2}$ and $∥ B ∥_{2} = ∥ b ∥_{2}$ . Expand $\hat{A} ⊙ \hat{B} - A ⊙ B = d_{A} ⊙ B + A ⊙ d_{B} + d_{A} ⊙ d_{B}$ . Applying [F2] termwise bounds this transform perturbation by $(2 ρ + ρ^{2}) ∥ a ∥_{2} ∥ b ∥_{2}$ . For the rounded products, [F4] gives coordinate errors $∣ e_{p, k} ∣ \leq ϵ_{p} ∣ {\hat{A}}_{k} ∣ ∣ {\hat{B}}_{k} ∣$ ; squaring and using [F2] yields $∥ e_{p} ∥_{2} \leq ϵ_{p} ∥ \hat{A} ∥_{2} ∥ \hat{B} ∥_{2} \leq ϵ_{p} (1 + ρ)^{2} ∥ a ∥_{2} ∥ b ∥_{2}$ . Thus $∥ \hat{p} - A ⊙ B ∥_{2} \leq λ ∥ a ∥_{2} ∥ b ∥_{2}$ and, by [F2], $∥ \hat{p} ∥_{2} \leq (1 + λ) ∥ a ∥_{2} ∥ b ∥_{2}$ .

2.1step 1.1F1F3F6algebra

Let $z = F_{M}^{- 1} (A ⊙ B) = c / \sqrt{M}$ and let $\hat{z}$ be the computed inverse result before final scaling. By Parseval [F3], $∥ F_{M}^{- 1} (\hat{p} - A ⊙ B) ∥_{2} = ∥ \hat{p} - A ⊙ B ∥_{2}$ ; by the inverse part of [F1], $∥ \hat{z} - F_{M}^{- 1} \hat{p} ∥_{2} \leq ρ ∥ \hat{p} ∥_{2}$ . Subadditivity [F6] and step 1.1 give $∥ \hat{z} - z ∥_{2} \leq [λ + ρ (1 + λ)] ∥ a ∥_{2} ∥ b ∥_{2}$ . Also $∥ \hat{z} ∥_{2} \leq (1 + ρ) ∥ \hat{p} ∥_{2} \leq (1 + ρ) (1 + λ) ∥ a ∥_{2} ∥ b ∥_{2}$ .

3.1step 2.1F5F6algebra∎

The exact scaled target is $\sqrt{M} z = c$ . By [F5], scaling $\hat{z}$ incurs at most $ϵ_{s} \sqrt{M} ∥ \hat{z} ∥_{2}$ additional error. Subadditivity [F6] and step 2.1 therefore give $∥ \hat{c} - c ∥_{2} \leq \sqrt{M} [λ + ρ (1 + λ) + ϵ_{s} (1 + ρ) (1 + λ)] ∥ a ∥_{2} ∥ b ∥_{2}$ , as claimed. In linear mode, Zero padding embeds linear convolution in cyclic convolution without aliasing identifies the exact first $n + m - 1$ cyclic entries with the linear convolution, and coordinate restriction cannot increase the Euclidean norm.

TheoremStatement: Literature-sourcedProof: AI-generatedOpen item page →

Bluestein chirp-phase error is bounded separately from FFT-convolution error

Statement

Let $N \geq 2$ be prime, $M \geq 2 N - 1$ a power of two, and $K = 2 N - 1$ . Use the Bluestein representation of Bluestein’s chirp identity embeds an arbitrary DFT in a padded cyclic convolution. Suppose every certified chirp phase has absolute error at most $0 < τ \leq 1 / 4$ , $u \leq 1 / 16$ , $E \geq 0$ , the modulation and demodulation products succeed under the six-operation complex model, and the computed cyclic convolution satisfies $∥ \hat{c} - \hat{a} *_{M} \hat{b} ∥_{2} \leq \sqrt{M} E ∥ \hat{a} ∥_{2} ∥ \hat{b} ∥_{2}$ . Put $δ = τ + 3 u (1 + τ)$ and $H = δ (2 + δ) + E (1 + δ)^{2}$ . Then the computed unnormalised length- $N$ DFT output $\hat{X}$ satisfies

∥ \hat{X} - X ∥_{2} \leq \sqrt{M K} [δ + (1 + δ) H] ∥ x ∥_{2} .

The term $δ$ records chirp enclosure and chirp multiplication roundoff; $E$ records the separate FFT, pointwise-product, inverse, and cyclic-convolution scaling error. No mixed-radix stability claim is inferred.

Facts & Assumptions

Given: A prime $N \geq 2$ , $M \geq 2 N - 1$ a power of two, $x \in C^{N}$ , a successful Bluestein computation, and the bounds $τ, u, E$ in the Statement.

[F1]

The exact unnormalised DFT is $X_{k} = γ_{k} (a *_{M} b)_{k}$ , where $a_{j} = α_{j} x_{j}$ , $α_{j} = e^{- π i j^{2} / N}$ , $b_{[t]_{M}} = β_{t} = e^{π i t^{2} / N}$ , and $γ_{k} = e^{- π i k^{2} / N}$ (Bluestein’s chirp identity embeds an arbitrary DFT in a padded cyclic convolution).

[F2]

The certified phase algorithm gives $∣ {\hat{α}}_{j} - α_{j} ∣ \leq τ$ , $∣ {\hat{β}}_{t} - β_{t} ∣ \leq τ$ , and $∣ {\hat{γ}}_{k} - γ_{k} ∣ \leq τ$ ; all exact chirp moduli equal $1$ (Certified phase output encloses the exact unit-circle value, $\exp (x + i y) = e^{x} (\cos y + i \sin y)$ , $∣ \exp (x + i y) ∣ = e^{x}$ , and $e^{i π} + 1 = 0$ ).

[F3]

A successful complex multiplication contributes at most $3 u ∣ z ∣ ∣ w ∣$ absolute error (Six-operation complex multiplication has a normwise roundoff bound).

[F4]

For length- $M$ vectors, $∥ a *_{M} b ∥_{2} \leq \sqrt{M} ∥ a ∥_{2} ∥ b ∥_{2}$ , and Parseval preserves norms under $F_{M}$ and its inverse (The DFT turns cyclic convolution into a scaled pointwise product, Finite Parseval and Plancherel identity for the unitary DFT, The counting inner product on $C^{Z / N Z}$ ).

[F5]

The computed cyclic convolution satisfies $∥ \hat{c} - \hat{a} *_{M} \hat{b} ∥_{2} \leq \sqrt{M} E ∥ \hat{a} ∥_{2} ∥ \hat{b} ∥_{2}$ (Computed FFT convolution error separates transform, product, inverse, and final-scale terms).

[F6]

For length- $M$ vectors, $∥ u ⊙ v ∥_{2} \leq ∥ u ∥_{2} ∥ v ∥_{2}$ : after squaring, the product $(\sum_{k} ∣ u_{k} ∣^{2}) (\sum_{k} ∣ v_{k} ∣^{2})$ contains the terms $∣ u_{k} ∣^{2} ∣ v_{k} ∣^{2}$ and only additional nonnegative cross terms.

[F7]

The Euclidean norm on $C^{M}$ is subadditive (Cauchy–Schwarz: $∣ ⟨ u, v ⟩ ∣ \leq ∥ u ∥ ∥ v ∥$ , with equality exactly for linearly dependent vectors, The counting inner product on $C^{Z / N Z}$ ).

Proof

1.1givenF2F3algebra

Let ${\hat{a}}_{j} = {fl}_{C} ({\hat{α}}_{j} x_{j})$ on $0 \leq j < N$ and zero elsewhere. By [F2]–[F3], $∣ {\hat{a}}_{j} - a_{j} ∣ \leq (τ + 3 u (1 + τ)) ∣ x_{j} ∣ = δ ∣ x_{j} ∣$ , so $∥ \hat{a} - a ∥_{2} \leq δ ∥ x ∥_{2}$ and $∥ \hat{a} ∥_{2} \leq (1 + δ) ∥ x ∥_{2}$ . The vector $b$ has exactly $K$ unit-modulus entries; [F2] gives $∥ \hat{b} - b ∥_{2} \leq τ \sqrt{K} \leq δ \sqrt{K}$ and $∥ \hat{b} ∥_{2} \leq (1 + δ) \sqrt{K}$ .

2.1F1F2F4F5F6F7step 1.1algebra

By bilinearity, $\hat{a} *_{M} \hat{b} - a *_{M} b = (\hat{a} - a) *_{M} \hat{b} + a *_{M} (\hat{b} - b)$ . Applying the exact convolution formula [F4], Parseval, and the pointwise product bound [F6] gives $∥ \hat{a} *_{M} \hat{b} - a *_{M} b ∥_{2} \leq \sqrt{M} (∥ \hat{a} - a ∥_{2} ∥ \hat{b} ∥_{2} + ∥ a ∥_{2} ∥ \hat{b} - b ∥_{2}) \leq \sqrt{M K} δ (2 + δ) ∥ x ∥_{2}$ by step 1.1. The computed-convolution error [F5] is at most $\sqrt{M K} E (1 + δ)^{2} ∥ x ∥_{2}$ . Thus, with $c = a *_{M} b$ and $\hat{c}$ the computed convolution, $∥ \hat{c} - c ∥_{2} \leq \sqrt{M K} H ∥ x ∥_{2}$ ; [F4] and step 1.1 give $∥ c ∥_{2} \leq \sqrt{M K} ∥ x ∥_{2}$ , so $∥ \hat{c} ∥_{2} \leq \sqrt{M K} (1 + H) ∥ x ∥_{2}$ by [F7].

3.1F1F2F3F7step 2.1algebra∎

The output is the componentwise rounded product ${\hat{X}}_{k} = {fl}_{C} ({\hat{γ}}_{k} {\hat{c}}_{k})$ . By [F2]–[F3], $∣ {\hat{X}}_{k} - γ_{k} {\hat{c}}_{k} ∣ \leq δ ∣ {\hat{c}}_{k} ∣$ ; since $∣ γ_{k} ∣ = 1$ , step 2.1 and subadditivity [F7] give $∥ \hat{X} - X ∥_{2} \leq δ ∥ \hat{c} ∥_{2} + ∥ \hat{c} - c ∥_{2} \leq \sqrt{M K} [δ (1 + H) + H] ∥ x ∥_{2} = \sqrt{M K} [δ + (1 + δ) H] ∥ x ∥_{2}$ . This is the claimed separated chirp and convolution bound.

5 · Examples, counterexamples and false statements

None yet.

Computed fast transforms and convolution

1 · Prerequisites

2 · Summary

3 · Logical flowchart

4 · Definitions, theorems and proofs

DCT-II and DST-I with fixed even and odd extension conventions

Definition

Bluestein’s chirp identity embeds an arbitrary DFT in a padded cyclic convolution

Statement

Facts & Assumptions

Proof

Mixed-radix Cooley–Tukey factorization with fixed index maps

Statement

Facts & Assumptions

Proof

Zero padding embeds linear convolution in cyclic convolution without aliasing

Statement

Facts & Assumptions

Proof

Six-operation complex multiplication has a normwise roundoff bound

Statement

Facts & Assumptions

Proof

DCT-II is orthogonal with the stated endpoint normalization and even extension

Statement

Facts & Assumptions

Proof

DST-I is orthogonal with the stated endpoint normalization and odd extension

Statement

Facts & Assumptions

Proof

DCT-II samples extract Chebyshev coefficients at first-kind roots

Statement

Facts & Assumptions

Proof

The interval evaluator's arctangent definition uses the canonical Fourier pi

Statement

Facts & Assumptions

Proof

Certified dyadic enclosures for unit-circle phases

Definition

Remarks

Computed unitary radix-two FFT with certified twiddle phases

Definition

Certified phase output encloses the exact unit-circle value

Statement

Facts & Assumptions

Proof

Computed FFT cyclic convolution and zero-padded linear convolution

Definition

A computed unitary radix-two butterfly has error at most (10u + εtw) times its input norm

Statement

Facts & Assumptions

Proof

Arbitrary-length mixed-radix and Bluestein DFT

Definition

Computed unitary radix-two FFT has the multiplicative stage error bound

Statement

Facts & Assumptions

Proof

Mixed-radix and Bluestein reductions compute every finite DFT in O(N log N) arithmetic

Statement

Facts & Assumptions

Proof

Computed FFT convolution error separates transform, product, inverse, and final-scale terms

Statement

Facts & Assumptions

Proof

Bluestein chirp-phase error is bounded separately from FFT-convolution error

Statement

Facts & Assumptions

Proof

5 · Examples, counterexamples and false statements

Sources