diff --git a/research/periodic-gaussian-direct-solve/README.md b/research/periodic-gaussian-direct-solve/README.md new file mode 100644 index 0000000000..b6ee76a87a --- /dev/null +++ b/research/periodic-gaussian-direct-solve/README.md @@ -0,0 +1,214 @@ +# A certified direct consumer for the periodic Gaussian reference + +**Constructive finite-source result with an all-input error proof; no new +kappa.** The complete inverse action for the original +`Reference(384,1,2,4,18,320)` can use a sparse direct factorization instead of +119 whole-period updates. It retains all 385 rows, the true periodic entries, +hole 192, arbitrary signed complex inputs, and the original 33-bit input and +output grid. The internal factor and solve grid has 48 fractional bits. + +The solve uses **1,398 integer products**, compared with **731,136** in the +weighted recurrence of [PR #5216](https://github.com/CrocSwap/integer-mult-bounds/pull/5216). +Preparation adds 319 factor products, 349 factor divisions and 668 exact +factor-certificate products. Both methods retain the original coefficient +producer. Its cost is substantial and is not a free input oracle. +The new solve performs 2,310 component divisions/truncations, including final +rounding back to the source grid. These product counts are not a comparison +of total bit operations. [TIMING.json](TIMING.json) records local sequential +wall-clock trials, including direct factor preparation and reuse separately. + +This is a host-array consumer with a verified finite source. It is not a +physical tape implementation, a production-parameter theorem, native supplier +admission or a new multiplication exponent. The selected kappa remains +`60403/20000000`. + +## Predecessors and changed contract + +Jose Manuel supplied the immutable rational Reference. JoyboyBrian's +[PR #5159](https://github.com/CrocSwap/integer-mult-bounds/pull/5159) supplies the +complete periodic coefficient table, all-output recurrence, source gap, +input ownership and analytic perturbation budget. His PR #5216 reduces the +recurrence from 275 to 119 updates. Both are retained and replayed here. + +Dugongue's [PR #4554](https://github.com/CrocSwap/integer-mult-bounds/pull/4554) +already gives a native band/Woodbury construction. Ryanson Jonathan's +[band audit](../pr4554-banded-inverse-check/README.md) checks that construction. +Classical Gaussian elimination, diagonal dominance and periodic band solving +are not new results of this package. The existing `band_solve` executable +requires `||A-I||_infinity<=1/4`; its proof starts with a small-norm hypothesis. +The actual source here has a funded gap about 0.06334, and its stored +off-diagonal row norm exceeds one half. That executable contract does not +admit this fixture. No check in that implementation is weakened. + +The new result is an actual source-compatible direct consumer with a checked +factor residual. The certificate replaces the small-norm premise at this +boundary. The sparse kernel retains cyclic fill and checks every original +stored coefficient. It does not discard corner entries to claim a band. +All predecessor files and source pins remain unchanged. + +## Certificate theorem + +Use the maximum over coordinates of `|Re z|+|Im z|`. Let the exact source +be `H=I+A`, where `A` is real, nonnegative, has zero diagonal, and its row +sums are at most `1-d`, with `d>0`. Downward coefficient rounding to grid +`2^-P` gives `Hhat=I+Ahat`, with the same row bound and + +```text +||H-Hhat|| <= delta = m 2^-P, +``` + +where `m` bounds the number of stored off-diagonal entries per row. The +original producer proves these relations for its supplied source. The +verifier checks its bytes, literal support and original rational entries. + +Let the stored real factors be unit lower triangular `L` and upper +triangular `U` with positive diagonal. No claim about their quality is +assumed from the elimination algorithm. Multiply their sparse rows exactly +and check the coefficient residual + +```text +eta = ||LU-Hhat|| < d. +``` + +The checker includes every nonzero product and every original input entry. +Let `l=||L||`, `v=||U||`, `w=max_i U_ii` and `h=2^-b`, where `b>=P` is the +internal factor grid. These constants are computed with exact integers and +rational arithmetic. In this package `P=33`, `b=48` and `m=8`. + +Row diagonal dominance gives `||H^-1||,||Hhat^-1|| <= 1/d`. For example, +choose a coordinate of maximum complex-l1 norm and apply the triangle +inequality to that row. The same argument, with row perturbation at most +`eta`, gives + +```text +||(LU)^-1|| <= 1/(d-eta). +``` + +Let `u` be an arbitrary input on the source grid with norm at most `Y`. +Forward and backward substitution accumulate each dot product exactly and +round the solved component once toward zero. Thus their literal outputs +`z,y` satisfy + +```text +Lz = u + e_f, ||e_f|| <= 2h, +Uy = z + e_b, ||e_b|| <= 2wh. +``` + +There are two real components per coordinate. The second residual includes +the diagonal multiplier. Consequently + +```text +LUy = u + e_f + L e_b, +rounding = 2h(1+lw), +||y|| <= B = (Y+rounding)/(d-eta). +``` + +Also `L^-1=U(LU)^-1`, so +`||z|| <= v(Y+2h)/(d-eta)`. These bounds apply to the full solve and all +stored coordinates. They also bound each partial dot product by the row +absolute sum times the coordinate bound. The kernel checks its integer +word limit and these coordinate limits during execution. Preparation also +checks the uniform accumulator bound +`2^(2b) max(Y+(l-1)Z, Z+vB, Y, Z, B, 1) < 2^width`, where +`Z=v(Y+2h)/(d-eta)`. Thus every allowed input fits the accepted solve word +limit; input-dependent overflow is not hidden in the theorem. + +Now compare the computed output with the original source solution: + +```text +||y-H^-1 u|| <= ((delta+eta)B+rounding)/d. +``` + +The final two components are truncated to the original `P`-bit grid, adding +at most `2*2^-P`. Preparation accepts only when the sum is at most +`2^(-t)/2`. This proves the requested source error for **every** admitted +input, not only the finite test vectors. If the inherited analytic matrix +has reference error `E= 0 else -1)*(abs(a)//b) + + +def bounded(value, width): + need(abs(value) < 1 << width, 'arithmetic word overflow') + return value + + +@dataclass(frozen=True) +class Prepared: + lower: tuple + upper: tuple + bits: int + source_bits: int + width: int + radius: int + target: int + certificate: tuple + setup: tuple + + +def factors(rows, source_bits, bits, width): + n = len(rows) + need(n > 0 and type(source_bits) is int and 1 <= source_bits <= bits < width, 'factor dimensions') + scale = 1 << bits + a = [] + for i,row in enumerate(rows): + need(len({j for j,v in row}) == len(row), 'duplicate coefficient') + need(all(type(j) is int and 0 <= j < n and j != i and type(v) is int and v >= 0 + for j,v in row), 'coefficient contract') + a.append({i:scale, **{j:bounded(v << (bits-source_bits),width) for j,v in row if v}}) + columns = [set() for _ in range(n)] + for i,row in enumerate(a): + for j in row: + if i > j: + columns[j].add(i) + lower = [{} for _ in range(n)] + products = divisions = 0 + active = peak = sum(map(len,a)) + for k in range(n): + pivot = a[k][k] + need(pivot > 0, 'nonpositive pivot') + right = sorted((j,v) for j,v in a[k].items() if j > k) + for i in sorted(columns[k]): + value = a[i].pop(k) + active -= 1 + quotient = trunc(bounded(value*scale,width),pivot) + divisions += 1 + if quotient: + lower[i][k] = bounded(quotient,width) + active += 1 + for j,value in right: + prior = j in a[i] + new = bounded(a[i].get(j,0)-trunc(bounded(quotient*value,width),scale),width) + products += 1 + if new: + a[i][j] = new + active += int(not prior) + if i > j: + columns[j].add(i) + else: + a[i].pop(j,None) + active -= int(prior) + if i > j: + columns[j].discard(i) + peak = max(peak,active) + columns[k].clear() + l = tuple(tuple(sorted(row.items())) for row in lower) + u = tuple(tuple(sorted(row.items())) for row in a) + bill = dict(products=products,divisions=divisions,coefficient_updates=products, + lower_records=sum(map(len,l)),upper_records=sum(map(len,u)), + active_factor_coefficient_records_peak=peak) + return l,u,bill + + +def certify(rows, source_bits, lower, upper, bits, gap, radius, target): + """Check every input coefficient in L U using exact integers.""" + n = len(rows) + need(len(lower) == len(upper) == n and n > 0, 'certificate dimensions') + need(type(radius) is int and radius >= 1 and type(target) is int and target >= 1, + 'radius/target contract') + need(0 < gap <= 1 and bits >= source_bits, 'certificate scale/gap') + scale = 1 << bits + l = [dict(row) for row in lower] + u = [dict(row) for row in upper] + for i in range(n): + need(len({j for j,v in rows[i]}) == len(rows[i]) and + all(type(j) is int and 0 <= j < n and j != i and type(v) is int and v >= 0 + for j,v in rows[i]), 'certificate input contract') + need(len(l[i]) == len(lower[i]) and len(u[i]) == len(upper[i]), 'duplicate factor entry') + need(all(type(j) is int and 0 <= j < i and type(v) is int for j,v in lower[i]), 'lower triangular contract') + need(all(type(j) is int and i <= j < n and type(v) is int for j,v in upper[i]), 'upper triangular contract') + need(u[i].get(i,0) > 0, 'positive factor diagonal') + need(sum(v for j,v in rows[i]) <= (1-gap)*(1<