Computed unitary radix-two FFT with certified twiddle phases

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=2^m$ with $m\ge 0$, $x\in {\mathbb{C}}^N$ whose real and imaginary coordinates all belong to $\mathbb{F}$, a direction forward or inverse, an exact rational twiddle tolerance $0<\tau \le 1\mathrm{/}4$, integer exponent limits $e_{\min }\le e_{\max }$ for the dyadic phase format, 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 }$. Output a vector in ${\mathbb{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 ${\widehat{w}}_k$, require its real and imaginary coordinates to belong to $\mathbb{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\le k<N\mathrm{/}2$, call Certified dyadic enclosures for unit-circle phases on $(q,\tau ,e_{\min },e_{\max },p_0,p_{\max },J)$ with $q=-2k\mathrm{/}N$, obtaining $({\widehat{w}}_k,{ϵ}_k)$ approximating $w_k=e^{-2\pi ik\mathrm{/}N}$. Compute the six-real-operation complex product $t_k={\mathrm{fl}}_{\mathbb{C}}({\widehat{w}}_k{O}_k)$, componentwise rounded sums $s_k^{+}={\mathrm{fl}}_{\mathbb{C}}(E_k+t_k)$ and $s_k^{-}={\mathrm{fl}}_{\mathbb{C}}(E_k-t_k)$, and $\widehat{c}=\mathrm{fl}(\sqrt{1\mathrm{/}2})$. Set output entries $k$ and $k+N\mathrm{/}2$ to the componentwise rounded real scalings ${\mathrm{fl}}_{\mathbb{C}}(\widehat{c}{s}_k^{+})$ and ${\mathrm{fl}}_{\mathbb{C}}(\widehat{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.