Arbitrary-length mixed-radix and Bluestein DFT

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 $\mathbb{F}$. Input $N\ge 1$, $x\in {\mathbb{C}}^N$ whose real and imaginary coordinates all belong to $\mathbb{F}$, an exact rational phase tolerance $0<\tau \le 1\mathrm{/}4$, integer dyadic exponent limits $e_{\min }\le e_{\max }$, an initial phase precision $p_0\ge 2$, a finite precision cap $p_{\max }\ge p_0$, a finite elementary-series cap $J\ge 1$, and a maximum radix-two recursion depth $d_{\max }\ge 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 $\mathbb{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\pi ijk\mathrm{/}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.

1. For $N=1$, return $x_0$.
2. If $N=2^m$ with $m\ge 1$, call Computed unitary radix-two FFT with certified twiddle phases in forward mode with controls $(\tau ,e_{\min },e_{\max },p_0,p_{\max },J,d_{\max })$ and multiply each output by the computed scalar $\mathrm{fl}(\sqrt{N})$. A failed subcall or scale operation returns its status.
3. 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\mathrm{/}r$. For $u=0,\dots ,r-1$, recursively call this procedure with the same controls on the strided vector $(x_{u+rv}{)}_{v=0}^{s-1}$ to obtain an approximate unnormalised length-$s$ DFT $U_u$. For each $0\le k_2<s$ and $0\le u<r$, call Certified dyadic enclosures for unit-circle phases on $(q,\tau ,e_{\min },e_{\max },p_0,p_{\max },J)$ with $q=-2u{k}_2\mathrm{/}N$ to obtain the twiddle enclosure, then multiply its returned phase by $U_u(k_2)$ using ${\mathrm{fl}}_{\mathbb{C}}$. For each $0\le 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\le k_1<r$ store its output $k_1$ at index $k_2+s{k}_1$. Propagate any failed phase or arithmetic status.
4. If $N$ is prime, set $M$ to the least power of two with $M\ge 2N-1$. For each $0\le j<N$, call Certified dyadic enclosures for unit-circle phases on $(q,\tau ,e_{\min },e_{\max },p_0,p_{\max },J)$ with $q=-j^2\mathrm{/}N$ to obtain ${\widehat{\alpha }}_j$, and set $a_j={\mathrm{fl}}_{\mathbb{C}}({\widehat{\alpha }}_j{x}_j)$; set the remaining entries of $a\in {\mathbb{C}}^M$ to zero. For each $-(N-1)\le t\le N-1$, obtain ${\widehat{\beta }}_t$ by the same call with $q=t^2\mathrm{/}N$ and set $b_{[t{]}_M}={\widehat{\beta }}_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 $(\tau ,e_{\min },e_{\max },p_0,p_{\max },J,d_{\max })$. For each $0\le k<N$, obtain ${\widehat{\gamma }}_k$ by the same phase call with $q=-k^2\mathrm{/}N$ and set ${\widehat{X}}_k={\mathrm{fl}}_{\mathbb{C}}({\widehat{\gamma }}_k{\widehat{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.