Merge pull request #49 from falseywinchnet/tcq1n0-main/audit-accuracy-of-odd/prime-techniques

High-precision twiddles and pairwise summation for generalized Bruun prototype; add accuracy probe and report

Author: falseywinchnet

documentation/reports/odd_prime_accuracy.md

@@ -0,0 +1,132 @@
+# Odd and prime-powered generalized-Bruun accuracy audit
+
+## Scope and verdict
+
+This report audits the odd-prime and prime-power path in `scratch_genbruun_exact.py`, especially the condition-1 cascade frame introduced for composite odd Bruun nodes.  The result is mathematically correct as a DFT factorization, and the new cascade is the right structural fix for the previous Chebyshev inverse-coefficient blow-up.  It is not a cheat: it is the direct complex block DFT written as two real planes.
+
+The metrology result is now much closer to the target after pairwise reductions and high-precision-generated twiddles.  The current Python prototype is usually equal to or better than NumPy's mixed-radix/T-C style `rfft` on the sampled odd and prime-powered sizes, with normalized maximum-error ratios from **0.33x to 2.00x** NumPy's and most composite prime-power cases below NumPy.  It is not a formal correctly rounded FFT, and it is not an unconditional win for every prime leaf, but the correct-twiddle path extracts the main available accuracy benefit from the odd/prime design.
+
+## Mathematical structure
+
+The top-level real transform factors `z^N - 1` by repeatedly removing an odd prime `p` from the active length `D = pM`.  The DC branch is the block sum, while each conjugate branch is evaluated at `g = j / p` turns for `1 <= j <= (p - 1) / 2`.
+
+For a real block decomposition `R_i`, the condition-1 seed is
+
+```text
+seed_lo = sum_i R_i cos(2 pi i j / p)
+seed_hi = sum_i R_i sin(2 pi i j / p)
+X_seed = seed_lo - i seed_hi.
+```
+
+For a composite child with residual size `M`, the cascade keeps `X_seed` as two real arrays rather than converting back through Chebyshev inverse coefficients.  If `lo - i hi` is already the node value at parent angle `f`, then the child with index `t` and angle `phi = (f + t) / p` is
+
+```text
+child_lo = sum_i A_i cos(2 pi i phi) - B_i sin(2 pi i phi)
+child_hi = sum_i B_i cos(2 pi i phi) + A_i sin(2 pi i phi)
+```
+
+where `A_i` and `B_i` are the `p` block slices of `lo` and `hi`.  Equivalently,
+
+```text
+child_lo - i child_hi = sum_i (A_i - i B_i) exp(-2 pi i i phi).
+```
+
+That is exactly the length-`p` complex DFT across the block index, evaluated with unit-modulus twiddles.  The power-of-two tail then enters the normalized Bruun subtree without changing representation.
+
+## Correctness argument
+
+The implementation tracks every phase as a `fractions.Fraction` of a full turn.  Therefore every recursive child phase is an exact rational expression until the single local `sin` or `cos` evaluation.  Bin placement uses `round(N * f) mod N`, which is an integer for the rational phases generated by this factor tree.
+
+The cascade recurrence is algebraically exact over real arithmetic because it is just complex multiplication by `exp(-2 pi i i phi)` and addition.  The real and imaginary formulas above are the expanded real form of that complex expression.  The base case places `lo[0] - i hi[0]`; the power-of-two case passes the unchanged two-plane normalized seed into the existing normalized Bruun subtree.  Inductively, every branch evaluates the same DFT value that the natural odd-prime factorization would evaluate.
+
+## Error sources and drift mechanisms
+
+1. **Trigonometric evaluation.**  Each `_cos` and `_sin` call now evaluates the exact rational turn with mpmath at 160 decimal digits and stores the resulting binary64 value in a cache.  This models a C++ table-generation path where twiddles are generated with enough precision to be correctly rounded before the FFT runs.
+2. **Block summation.**  The prototype now uses pairwise array reductions in the odd-prime path.  There is still no compensated summation or Shewchuk expansion in the transform itself, so cancellation-heavy blocks can still lose low-order bits.
+3. **Recursive accumulation.**  Each odd stage adds `p` rounded vector terms per child.  Prime powers such as `3^6` revisit the same kind of stage many times, so first-order error grows with both depth and the number of nonzero block contributions.
+4. **Power-of-two tail.**  When `M` becomes a power of two, accuracy inherits the normalized Bruun subtree's rounding behavior.
+5. **Reference conversion.**  The metrology script evaluates the reference DFT with mpmath at 90 decimal digits and converts the final reference bins to binary64 complex values for error measurement.  RMS reductions use `math.fsum` to avoid report-side summation drift.
+
+The removed Chebyshev inverse path had a worse conditioning failure: coefficients proportional to `1 / sin(phi)` become large when a composite odd child is near a singular angle.  The cascade frame replaces those coefficients with unit-modulus twiddles, so its local coefficient condition is 1.  That fixes the main structural instability, but ordinary rounded summation remains.
+
+## Measured accuracy
+
+Command:
+
+```sh
+PYTHONPATH=. python scripts/odd_prime_accuracy_probe.py
+```
+
+Reference: direct mpmath DFT at 90 decimal digits.  Inputs: deterministic `standard_normal` vectors seeded by `N`.  Errors: absolute bin error divided by `max(max(abs(reference)), 1)`.
+
+| N | factors | max Bruun | max NumPy | RMS Bruun | RMS NumPy | max ratio |
+| ---: | --- | ---: | ---: | ---: | ---: | ---: |
+| 3 | 3 | 6.874e-18 | 2.062e-17 | 4.861e-18 | 1.458e-17 | 0.33 |
+| 5 | 5 | 7.630e-17 | 3.815e-17 | 5.395e-17 | 2.203e-17 | 2.00 |
+| 7 | 7 | 1.103e-16 | 9.865e-17 | 7.601e-17 | 6.287e-17 | 1.12 |
+| 9 | 3^2 | 1.518e-16 | 1.224e-16 | 9.148e-17 | 1.044e-16 | 1.24 |
+| 15 | 3 * 5 | 1.621e-16 | 2.017e-16 | 1.068e-16 | 1.305e-16 | 0.80 |
+| 21 | 3 * 7 | 1.091e-16 | 8.987e-17 | 7.221e-17 | 5.253e-17 | 1.21 |
+| 25 | 5^2 | 1.463e-16 | 1.573e-16 | 8.416e-17 | 8.778e-17 | 0.93 |
+| 27 | 3^3 | 1.361e-16 | 1.361e-16 | 8.610e-17 | 8.492e-17 | 1.00 |
+| 45 | 3^2 * 5 | 1.814e-16 | 2.250e-16 | 9.514e-17 | 1.250e-16 | 0.81 |
+| 75 | 3 * 5^2 | 1.621e-16 | 2.563e-16 | 7.171e-17 | 1.090e-16 | 0.63 |
+| 81 | 3^4 | 1.244e-16 | 1.760e-16 | 7.048e-17 | 8.732e-17 | 0.71 |
+| 125 | 5^3 | 1.524e-16 | 2.170e-16 | 7.643e-17 | 9.807e-17 | 0.70 |
+| 225 | 3^2 * 5^2 | 2.260e-16 | 2.452e-16 | 9.168e-17 | 1.036e-16 | 0.92 |
+| 243 | 3^5 | 2.267e-16 | 3.356e-16 | 9.013e-17 | 1.208e-16 | 0.68 |
+| 405 | 3^4 * 5 | 2.040e-16 | 2.325e-16 | 8.717e-17 | 8.971e-17 | 0.88 |
+| 729 | 3^6 | 3.046e-16 | 3.137e-16 | 8.514e-17 | 1.114e-16 | 0.97 |
+
+## Bounds and comparison with a typical T-C FFT
+
+A first-order floating-point model gives a local `p`-term child sum error on the order of `gamma_{O(p)} * sum_i |block_i|`, where `gamma_m = m eps / (1 - m eps)`.  Because the cascade twiddles have magnitude 1, the coefficient growth factor at the odd split is approximately 1 rather than `1 / |sin(phi)|`.  Across a path with odd radices `p_1, ..., p_s` and a power-of-two tail, the expected forward error is therefore proportional to
+
+```text
+O(eps * (sum_r p_r + C_bruun log2(M_tail))) * moderate data-growth factors.
+```
+
+That is the same qualitative stability class as a standard mixed-radix Cooley-Tukey FFT.  With correctly generated twiddles, the current prototype has a comparable or smaller constant than NumPy's implementation on most measured composite odd cases, while isolated prime leaves can still be worse in max-error ratio.
+
+## Conclusion
+
+The condition-1 cascade is mathematically correct and removes the composite odd-node Chebyshev singularity.  It should remain the candidate structure for a C++ port.  With pairwise reductions and correctly generated twiddles, the prototype now meets or exceeds NumPy/T-C accuracy on most sampled composite odd and prime-power cases, but it is still not a formal correctly rounded FFT and still has prime-leaf cases where NumPy wins by max error.  Before porting as an accuracy-superior implementation, the next work should evaluate compensated odd child reductions, FMA kernels, and a larger automated mpmath/Shewchuk-style regression suite.
+
+## Follow-up design notes: reducing accumulation
+
+The prototype now uses pairwise array summation for the odd-prime DC branch, condition-1 seed construction, and cascade child construction.  This changes the reduction depth from a left-associated `p - 1` additions to approximately `ceil(log2(p))` rounded levels.  It does not make the result correctly rounded, but it reduces the expected summation constant and better matches how a SIMD C++ implementation should reduce vector lanes and radix blocks.
+
+For the C++ port, block summation is likely manageable rather than a fundamental blocker if the implementation uses a fixed pairwise reduction tree:
+
+1. For small odd radices, keep the explicit radix kernel and sum in a balanced order.
+2. For larger prime radices, accumulate real and imaginary planes with a deterministic vector tree, then horizontally reduce lanes with pairwise or compensated lane reduction.
+3. Use FMA for `a * c - b * s` and `b * c + a * s` when the target permits it, because each child term then has one final rounding after the product-add expression instead of separate rounded multiply and add steps.
+4. Avoid recurrence-generated twiddles in the accuracy path.  Generate or table exact-rational phases so each twiddle is rounded independently from the intended phase.
+
+Correct rounding of the whole FFT is a stronger target than FFT-grade stability.  To claim correct rounding for every output bin, the implementation would need a reproducible error bound below one half ulp for each bin after all twiddle, product, and summation errors.  In practice that means either adaptive extended precision for difficult bins, expansion arithmetic/Shewchuk-style summation in the odd reductions, correctly rounded or tightly bounded sin/cos tables, and interval/error tracking through the Bruun subtree.  Pairwise SIMD reduction alone is not enough for a formal correct-rounding claim.
+
+## Pairwise-summation probe update
+
+A short probe run during the pairwise change showed the same pattern on the smaller sample:
+
+| N | factors | max Bruun | max NumPy | RMS Bruun | RMS NumPy | max ratio |
+| ---: | --- | ---: | ---: | ---: | ---: | ---: |
+| 3 | 3 | 1.555e-16 | 2.062e-17 | 1.101e-16 | 1.458e-17 | 7.54 |
+| 5 | 5 | 1.145e-16 | 3.815e-17 | 8.531e-17 | 2.203e-17 | 3.00 |
+| 9 | 3^2 | 3.643e-16 | 1.224e-16 | 1.816e-16 | 1.044e-16 | 2.98 |
+| 15 | 3 * 5 | 2.087e-16 | 2.017e-16 | 1.347e-16 | 1.305e-16 | 1.03 |
+| 27 | 3^3 | 6.805e-16 | 1.361e-16 | 3.153e-16 | 8.492e-17 | 5.00 |
+| 45 | 3^2 * 5 | 7.374e-16 | 2.250e-16 | 3.579e-16 | 1.250e-16 | 3.28 |
+| 75 | 3 * 5^2 | 5.730e-16 | 2.563e-16 | 2.518e-16 | 1.090e-16 | 2.24 |
+| 81 | 3^4 | 6.700e-16 | 1.760e-16 | 2.994e-16 | 8.732e-17 | 3.81 |
+| 125 | 5^3 | 4.514e-16 | 2.170e-16 | 2.019e-16 | 9.807e-17 | 2.08 |
+| 225 | 3^2 * 5^2 | 7.232e-16 | 2.452e-16 | 2.966e-16 | 1.036e-16 | 2.95 |
+
+The improvement is visible for radix-5-heavy cases, especially `N = 5`, `75`, `125`, and `225`.  Radix-3 powers are mostly unchanged because a three-term balanced sum still has two addition levels and the dominant error is often twiddle/product rounding or inherited subtree error.  The conclusion therefore remains conservative: pairwise summation is a necessary C++ design choice, but not sufficient by itself to beat NumPy/T-C across the odd-prime stack.
+
+
+## Correctly rounded twiddle update
+
+The final prototype round changes twiddle generation from `np.cos(2*pi*float(f))` and `np.sin(2*pi*float(f))` to a cached high-precision path.  Each exact rational turn is evaluated with mpmath at 160 decimal digits, converted once to binary64, and then reused from the cache.  This models the intended C++ assumption that plan-time or offline twiddle tables are generated correctly and stored as doubles before SIMD kernels consume them.
+
+This change extracts most of the remaining accuracy benefit.  The largest measured Bruun/NumPy max-error ratio dropped from 7.54 in the libm-twiddle pairwise prototype to 2.00 with high-precision twiddles, and the composite prime-power cases in the table are now generally at or below NumPy.  The remaining gap is no longer dominated by phase generation; it is mostly the final rounded product/add schedule, cancellation in small prime leaves, and the inherited power-of-two subtree.

scratch_genbruun_exact.py

@@ -5,6 +5,7 @@
 accumulation); leaf bins are exact integers round(N*f).
 """
 import numpy as np
+import mpmath as mp
 from fractions import Fraction as Fr
 from scratch_normbruun import rfft_normbruun
 
@@ -18,11 +19,57 @@ def _odd_prime(n):
         d += 2
     return n
 
-TWO_PI = 2*np.pi
+TWIDDLE_DPS = 160
+_TWIDDLE_CACHE = {}
+
+def _twiddle_key(kind, f):
+    phase = f % 1
+    return kind, phase.numerator, phase.denominator
+
+def _phase_to_mpf(f):
+    phase = f % 1
+    return mp.mpf(phase.numerator) / mp.mpf(phase.denominator)
+
+def _rounded_trig(kind, f):
+    """Return a high-precision-generated binary64 twiddle component.
+
+    This models the C++ design target: twiddle tables are generated offline or at
+    plan construction with enough precision that each stored double is the
+    correctly rounded value of the exact rational turn.
+    """
+    key = _twiddle_key(kind, f)
+    if key in _TWIDDLE_CACHE:
+        return _TWIDDLE_CACHE[key]
+
+    with mp.workdps(TWIDDLE_DPS):
+        angle = 2 * mp.pi * _phase_to_mpf(f)
+        if kind == "cos":
+            value = float(mp.cos(angle))
+        else:
+            value = float(mp.sin(angle))
+
+    _TWIDDLE_CACHE[key] = value
+    return value
+
 def _cos(f):  # f: Fraction turns
-    return np.cos(TWO_PI*float(f % 1))
+    return _rounded_trig("cos", f)
 def _sin(f):
-    return np.sin(TWO_PI*float(f % 1))
+    return _rounded_trig("sin", f)
+
+def _pairwise_sum(terms):
+    """Pairwise sum array terms to reduce source-order accumulation."""
+    items = list(terms)
+    if not items:
+        return 0.0
+    while len(items) > 1:
+        next_items = []
+        limit = len(items) - 1
+        for i in range(0, limit, 2):
+            next_items.append(items[i] + items[i + 1])
+        if len(items) % 2 == 1:
+            next_items.append(items[-1])
+        items = next_items
+    return items[0]
 
 def _cheb_exact(g, n):
     """Chebyshev reduction coeffs a_i=-U_{i-2}(c), b_i=U_{i-1}(c), c=cos(2pi g).
@@ -87,8 +134,8 @@ def reduce_bruun_cascade(lo, hi, M, f):
         B = [hi[i*Mp:(i+1)*Mp] for i in range(p)]
         for t in range(p):
             phi = (f + t) / p             # (theta + 2pi t)/p in turns
-            child_lo = sum(A[i]*_cos(i*phi) - B[i]*_sin(i*phi) for i in range(p))
-            child_hi = sum(B[i]*_cos(i*phi) + A[i]*_sin(i*phi) for i in range(p))
+            child_lo = _pairwise_sum(A[i]*_cos(i*phi) - B[i]*_sin(i*phi) for i in range(p))
+            child_hi = _pairwise_sum(B[i]*_cos(i*phi) + A[i]*_sin(i*phi) for i in range(p))
             reduce_bruun_cascade(child_lo, child_hi, Mp, phi)
 
     def reduce_bruun(r, twoM, f):       # factor angle = 2*pi*f
@@ -104,13 +151,13 @@ def reduce_minus(r, D, sigma):
             return
         p = _odd_prime(D); M = D // p
         R = [r[i*M:(i+1)*M] for i in range(p)]
-        reduce_minus(sum(R), M, sigma*p)
+        reduce_minus(_pairwise_sum(R), M, sigma*p)
         for j in range(1, (p-1)//2 + 1):
             g = Fr(j, p)
             # CONDITION-1 cascade seed: B*(natural) = (sum_i R_i cos(i*phi),
             # sum_i R_i sin(i*phi)), bounded coeffs (no 1/sin(phi) growth). Exact phase.
-            seed_lo = sum(R[i]*_cos(Fr(i*j, p)) for i in range(p))
-            seed_hi = sum(R[i]*_sin(Fr(i*j, p)) for i in range(p))
+            seed_lo = _pairwise_sum(R[i]*_cos(Fr(i*j, p)) for i in range(p))
+            seed_hi = _pairwise_sum(R[i]*_sin(Fr(i*j, p)) for i in range(p))
             if M == 1:                                   # prime leaf: X = seed_lo - i seed_hi
                 place(int(round(N*float(g % 1))) % N, seed_lo[0] - 1j*seed_hi[0])
             else:                                        # composite/pow2 M: stay in cascade frame

scripts/odd_prime_accuracy_probe.py

@@ -0,0 +1,122 @@
+#!/usr/bin/env python3
+"""Odd/composite generalized-Bruun accuracy probe.
+
+This is an intentionally slow metrology script for scratch_genbruun_exact.py.  It
+compares the generalized Bruun prototype and NumPy rfft against an mpmath DFT
+reference evaluated at high precision.  The final reductions use math.fsum so
+that the reported RMS values are not dominated by Python summation order.
+"""
+
+from __future__ import annotations
+
+import argparse
+import math
+from pathlib import Path
+import sys
+
+import mpmath as mp
+import numpy as np
+
+REPO_ROOT = Path(__file__).resolve().parents[1]
+if str(REPO_ROOT) not in sys.path:
+    sys.path.insert(0, str(REPO_ROOT))
+
+from scratch_genbruun_exact import rfft_gen_exact  # noqa: E402
+
+
+def factor_string(n: int) -> str:
+    factors: list[str] = []
+    d = 2
+    value = n
+    while d * d <= value:
+        count = 0
+        while value % d == 0:
+            value //= d
+            count += 1
+        if count == 1:
+            factors.append(str(d))
+        if count > 1:
+            factors.append(f"{d}^{count}")
+        d += 1
+    if value > 1:
+        factors.append(str(value))
+    return " * ".join(factors)
+
+
+def mpmath_rfft(x: np.ndarray, dps: int) -> np.ndarray:
+    mp.mp.dps = dps
+    n = int(x.shape[0])
+    result: list[complex] = []
+    for k in range(n // 2 + 1):
+        real_part = mp.mpf("0")
+        imag_part = mp.mpf("0")
+        for sample_index, sample_value in enumerate(x):
+            angle = -2 * mp.pi * k * sample_index / n
+            exact_sample = mp.mpf(str(float(sample_value)))
+            real_part += exact_sample * mp.cos(angle)
+            imag_part += exact_sample * mp.sin(angle)
+        result.append(complex(float(real_part), float(imag_part)))
+    return np.asarray(result, dtype=np.complex128)
+
+
+def row_for_size(n: int, dps: int) -> tuple[str, ...]:
+    x = np.random.default_rng(n).standard_normal(n)
+    reference = mpmath_rfft(x, dps)
+    generalized = rfft_gen_exact(x)
+    numpy_rfft = np.fft.rfft(x)
+
+    scale = max(float(np.max(np.abs(reference))), 1.0)
+    generalized_error = np.abs(generalized - reference) / scale
+    numpy_error = np.abs(numpy_rfft - reference) / scale
+
+    generalized_rms = math.sqrt(
+        math.fsum(float(v * v) for v in generalized_error) / len(generalized_error)
+    )
+    numpy_rms = math.sqrt(
+        math.fsum(float(v * v) for v in numpy_error) / len(numpy_error)
+    )
+
+    if float(np.max(numpy_error)) == 0.0:
+        ratio = math.inf
+    else:
+        ratio = float(np.max(generalized_error)) / float(np.max(numpy_error))
+
+    return (
+        str(n),
+        factor_string(n),
+        f"{float(np.max(generalized_error)):.3e}",
+        f"{float(np.max(numpy_error)):.3e}",
+        f"{generalized_rms:.3e}",
+        f"{numpy_rms:.3e}",
+        f"{ratio:.2f}",
+    )
+
+
+def main() -> int:
+    parser = argparse.ArgumentParser()
+    parser.add_argument("--dps", type=int, default=90)
+    parser.add_argument(
+        "sizes",
+        nargs="*",
+        type=int,
+        default=[3, 5, 7, 9, 15, 21, 25, 27, 45, 75, 81, 125, 225, 243, 405, 729],
+    )
+    args = parser.parse_args()
+
+    header = (
+        "N",
+        "factors",
+        "max_bruun",
+        "max_numpy",
+        "rms_bruun",
+        "rms_numpy",
+        "max_ratio",
+    )
+    print("\t".join(header))
+    for n in args.sizes:
+        print("\t".join(row_for_size(n, args.dps)), flush=True)
+    return 0
+
+
+if __name__ == "__main__":
+    raise SystemExit(main())