diff --git a/research/rational-fft-fixed-h/DENSE-PROOF.md b/research/rational-fft-fixed-h/DENSE-PROOF.md new file mode 100644 index 0000000000..3ee98912ec --- /dev/null +++ b/research/rational-fft-fixed-h/DENSE-PROOF.md @@ -0,0 +1,336 @@ +# Conditioned Toeplitz programs on fixed Hadamard stages + +I extend the completed FFT connection to the actual stored full Toeplitz +program and its stored inverse at PR4205 head +78c72e563173f3822fcb4efed9ceeccf05630258. I retain its source and proof +bytes unchanged in [full_source](vendor/full_source/PROVENANCE.json). +I use Vlado Vojdanovski's existing zero-safe words and routing algorithm. +I add the nonunitary amplification argument and bind its parameters to +the admitted root, coefficient preparation and executed stages. + +## Generic domain + +Take a field of $N=2^m$ coordinates and an ordered program of $L$ +matching layers. In each layer, allow disjoint one-wire nonzero scalars, +upper or lower shears, invertible diagonal pairs, swaps and the rational +unitary butterflies or stored adjoints from [PROOF.md](PROOF.md). +The caller supplies exact rational bounds $\Pi,E_0,R\ge1$ such that: + +- Every original logical prefix and its inverse has norm at most $\Pi$. +- The exact stored endpoint $U$ has norm at most $E_0$. +- The incoming vector satisfies $\|x\|_2\le R$ and has complex components + on the $2^{-F_0}$ grid, for integer $F_0\ge0$. + +I use Euclidean vector norms and induced operator norms. These premises +concern linear maps with the stored coefficients and unrounded data. +The rounded execution itself is nonlinear. + +For $N\ge2$, I retain +\[ +J=L(16m+5),\qquad M=L(8m+3),\qquad h=L(8m+2), +\] +for the primitive, diagonal and whole-field Hadamard stage counts. +For $N=1$, I execute each scalar or idle layer as one diagonal: +$J=M=L$ and $h=0$. I make no bit-zero Hadamard call. +Put +\[ +\Gamma=2^{\lceil\log_2(\Pi^2)\rceil},\qquad +C=2^{14}\Gamma\Pi. +\tag{1} +\] +At a logical layer boundary, $A_j=P_jP_{j-1}^{-1}$ gives +$\|A_j\|,\|A_j^{-1}\|\le\Pi^2\le\Gamma$. Since the layer is a +direct sum, each scalar and its reciprocal has modulus at most $\Gamma$, +and each shear parameter has modulus at most $\Gamma$. + +For target integer $t\ge1$ and $L>0$, I choose +\[ +F=\max\left\{F_0,\ +t+\left\lceil\log_2(16E_0MC^2(R+1))\right\rceil,\ +t+\left\lceil\log_2(16E_0CJN)\right\rceil\right\}. +\tag{2} +\] +I prepare the exact reference words, replace each $H$ by $H_0=H/2$, +double its following diagonal, then round the diagonal components to +the $F$ grid. I round the data after each emitted stage, including +routing and idle words. I retain distinct adjacent diagonal stages. +I prove the following error against the stored logical endpoint: +\[ +\|\widehat y-Ux\|_2 +\le 4E_0MC^2\,2^{-F}R+4E_0CJN\,2^{-F} +<\tfrac12\,2^{-t}. +\tag{3} +\] +For $L=0$, I return the input, choose $F=\max(F_0,1)$, and incur zero +error. I give the connection to the exact Toeplitz root below. + +## Zero-safe shear and scalar words + +For an upper shear parameter $v$, set $c=2\Gamma+1$ and $z=v+c$. +The reverse triangle inequality gives +\[ +\Gamma+1\le |z|\le3\Gamma+1, +\] +so the retained construction can call its nonzero helper twice, even +at $v=0$. In execution order it uses the seven diagonals +\[ +(-5/(4c),1),\quad(-3,1),\quad(2,1),\quad +(c/(8z),-1/6),\quad(-3,1),\quad(2,1),\quad(-z/10,-1/6) +\tag{4} +\] +with a raw $H$ between neighbors. The two helper endpoints compose +as $U(z)U(-c)=U(v)$. Reversing the symmetric factors gives the +ordinary transpose $U(v)^T$, hence the lower shear. I do not conjugate +the complex parameter. + +I bound contiguous segments by products of factor norms, replacing +factors smaller than one by one. In (4), the middle diagonal has +norm at most $1/4$ and inverse norm at most 12. The last diagonal +has norm at most $(3\Gamma+1)/10\le\Gamma$, and inverse norm at most +six. The first inverse diagonal has norm at most $3\Gamma$. +Together with $\|H\|=\sqrt2$ and $\|H^{-1}\|=1/\sqrt2$, I obtain +raw segment bounds $288\Gamma$ and $216\Gamma$ for the direct and +inverse maps. The same bounds hold after word reversal, because +ordinary transposition preserves the spectral norm. + +After the $H_0$ compensation, a contiguous segment differs from its +raw counterpart by a scalar in $\{1/2,1,2\}$. Thus $1024\Gamma$ +bounds the native shear segments and their inverses. For a diagonal +pair, the source appends three complete two-H identity words; +the raw segment and inverse bounds are $8\Gamma$, and the native +bounds are $16\Gamma$. Swaps and idle pairs have smaller bounds. +For rational butterflies and their stored adjoints, I use the 4096 +bound from PROOF.md. Therefore +\[ +4096\Gamma +\tag{5} +\] +bounds each native central word segment and its inverse. + +The compensated diagonal entries have modulus at most $6\Gamma$ +and reciprocal modulus at most $200\Gamma$. I complete each mixed +layer's unused wires in pairs, using its actual singleton scalars +and identity values on idle wires. I preserve the oriented order +of each existing two-wire gate. All-one diagonal padding would +multiply the inactive pair by eight, so I use complete identity words. + +The source routes input $i$ to output $p[i]$. Its $2m-1$ switch +layers have unitary completed endpoints and native prefix/inverse +norms at most $\sqrt2$, as proved in PROOF.md. I use the reversed +switch schedule for the return route. Within a logical layer, +the central word and the two routes contribute less than +$2^{14}\Gamma$ to a partial prefix or inverse prefix. The completed +earlier logical program contributes at most $\Pi$. Equation (1) +therefore bounds exact compiled prefixes and their inverses by $C$. +I do not multiply the fragment allowance once per logical layer. + +## Nonunitary endpoint and rounding + +Let $P_j$ be the exact compiled prefix. An exact suffix is +$UP_j^{-1}$, so I must retain its bound $E_0C$. +For a rounded diagonal perturbation $\Delta_j$, put $u=2^{-F}$ +and $\varepsilon=2u$. Component rounding implies +$\|\Delta_j\|\le\varepsilon$, regardless of the field size. +In exact-prefix coordinates, +\[ +T_j=P_j^{-1}\Delta_jP_{j-1},\qquad +\|T_j\|\le C^2\varepsilon,\qquad +v:=MC^2\varepsilon<1/2. +\] +For any consecutive subproduct of the factors $I+T_j$, its norm +is at most $e^v<2$ and its distance from identity is at most $2v$. +Multiplying back by the endpoints gives: +\[ +\|\widehat U-U\|\le2E_0MC^2\varepsilon,\quad +\|\widehat P_j\|\le2C,\quad +\|\widehat S_j\|\le2E_0C,\quad +\|\widehat P_{b:a}\|\le2C^2. +\tag{6} +\] +Here the hatted maps use rounded coefficients and unrounded data. +The segment assertion follows by writing the segment in the same +exact-prefix coordinates; it does not require inverting a rounded +prefix. + +Each data rounding changes at most $N$ complex entries, with +Euclidean error at most $2Nu$. I propagate those residuals through +the perturbed suffixes in (6). I obtain (3), including the errors +from the completed identity words. In (2), the first precision +term pays strictly less than $2^{-t}/4$, because $R/(R+1)<1$; +the second pays at most $2^{-t}/4$. + +I retain a discriminating counterexample to omission of $E_0$. +Take $q=2^{16}$ and the three two-wire logical layers +\[ +\operatorname{diag}(q^{-1},q),\quad U(0),\quad +\operatorname{diag}(q^2,q^{-2}). +\] +The logical prefix and inverse-prefix bound is $\Pi=q$ and the +endpoint is $U=\operatorname{diag}(q,q^{-1})$, with $E_0=q$. +Use $\Gamma=q^2$, $c=2q^2+1$, and the upper word for the zero shear. +After its first central diagonal, the exact prefix is +\[ +P=\operatorname{diag}(-5/(4cq),q). +\] +Consequently +\[ +\|UP^{-1}\|=\frac{4cq^2}{5},\qquad +\frac{\|UP^{-1}\|}{C} +=\frac{32}{5}+\frac{1}{5\cdot2^{28}}>1. +\tag{7} +\] +The corrected bound $E_0C$ covers this suffix. I check the actual +retained word and routing prefix at this cut. This example refutes +the attempted unit-endpoint extension; it does not contradict the +completed FFT result. + +## Connection to the retained full program + +I admit an exact rational Toeplitz matrix $T$ of dyadic size $N$. +The retained admission routine checks its Toeplitz shape, inverts +each leading principal section and obtains component row/column +envelopes $\lambda,Q\ge1$. These imply +$\|T\|\le\lambda$ and inverse norm at most $Q$ for each leading +section. An invertible full matrix with a singular leading section +does not meet this admission. + +I invoke the unchanged full allocation and actual recursive plan +at target $t$ and radius $R$. I bind their matrix, budget, wire +order and both stored plans. Let $\chi=\lambda Q$, +$\mathcal B=256N^3\chi^{12}$ and $a=3N-2$. +The retained [FULL-PROOF.md](vendor/full_source/FULL-PROOF.md), +equations (12) and (14), gives +\[ +\|\widetilde T-T\|,\|\widetilde T^{-1}-T^{-1}\| +\le e=2a\mathcal B^3\eta,\qquad +\Pi=2\max(1,\lambda,Q)\mathcal B\mathcal H. +\tag{8} +\] +Here $\mathcal H$ is the retained generator/scalar fragment bound. +The source field execution_bound supplies $\Pi$ for stored logical +prefixes, suffixes and their inverses, in either direction. +I retain that author-proved dependency rather than interpret the +source's rounded-data amplitude field as an operator norm. + +I use $E_0=\lambda+e$ for the forward program and $E_0=Q+e$ for +its stored inverse. I compile that inverse's actual reversed gate +schedule. Nonunitary scalar leaves require reciprocals, so a global +adjoint would give a different endpoint. + +The retained choice +$\eta=2^{-s}$, with +$s=t+\lceil\log_2(4a\mathcal B^3(R+1))\rceil$, implies +\[ +eR\le\frac{2^{-t}R}{2(R+1)}<\tfrac12\,2^{-t}. +\] +I add this operator defect times $R$ to (3), obtaining +\[ +\boxed{\|\widehat y-T_I^{\pm1}x\|_2 +\le eR+4E_0MC^2uR+4E_0CJNu<2^{-t}.} +\tag{9} +\] +For the declared ordered physical wires $I$, I write $T_I$ for +the conjugate of $T$ by that coordinate permutation. The forward +direction uses $T_I$ and the reverse uses $T_I^{-1}$. This identification +preserves the norms and applies the exact root to the intended wire order. +I discard the predecessor's data-rounding allowance for this purpose; +the new stage schedule incurs the rounding term in (9). +I authenticate the source bytes and preserve source assertions under +optimized Python. Runtime seals reject changes to the admitted root, +request bounds, logical plans, prepared pages and native precision contract. + +## Widths, setup and remaining physical obligations + +Using the segment bound in (6), I bound actual states by +\[ +2CR+4C^2JNu +\le 2CR+\frac{C2^{-t}}{4E_0} +", "exec", optimize=0), module.__dict__) + loaded[name] = module + finally: + for name, old in previous.items(): + if old is None: + sys.modules.pop(name, None) + else: + sys.modules[name] = old + self.full = loaded["full"] + self.raw = MappingProxyType(dict(self.raw)) + + +class NativeProgram(FFT): + """A prefix-bound compiler; callers must establish the declared logical bound.""" + CONTRACT_FIELDS = ("size", "radius", "target", "incoming", "prefix", "endpoint", "upstream_error", + "logical_layers", "stage_count", "diagonal_count", "h_count", "gamma_bits", "gamma", + "compiled_bound", "bits", "coefficient_bits", "coefficient_error", "rounding_error", + "native_error", "total_error", "amplitude", "temporary_amplitude", "integer_bits", + "temporary_bits") + + def __init__(self, source, size, plan, prefix, endpoint, radius, target, incoming=0, *, + upstream_error=0, context=None, data_bits=None): + need(type(size) is int and size > 0 and size & (size - 1) == 0, "dyadic native field") + need(all(type(x) in (int, Q) and x >= 1 for x in (prefix, endpoint, radius)), + "exact positive circuit envelopes") + need(type(target) is int and target >= 1 and type(incoming) is int and incoming >= 0, + "native target and incoming grid") + need(type(upstream_error) in (int, Q) and upstream_error >= 0, "upstream operator defect") + self.source, self.e, self.context = source, source.e, context + self.size, self.radius, self.target, self.incoming = size, Q(radius), target, incoming + self.prefix, self.endpoint, self.upstream_error = Q(prefix), Q(endpoint), Q(upstream_error) + self.plan = [list(layer) for layer in plan] + self.reference_hash = plan_hash(source, self.plan) + self.logical_layers = len(self.plan) + m = size.bit_length() - 1 + self.stage_count = self.logical_layers * (16 * m + 5) if m else self.logical_layers + self.diagonal_count = self.logical_layers * (8 * m + 3) if m else self.logical_layers + self.h_count = self.logical_layers * (8 * m + 2) if m else 0 + self.gamma_bits = max(0, self.e.ceil_log2(self.prefix ** 2)) + self.gamma = 1 << self.gamma_bits + self.compiled_bound = (1 << 14) * self.gamma * self.prefix + c, e0, r = self.compiled_bound, self.endpoint, self.radius + minimum = max(incoming, target + self.e.ceil_log2(16 * e0 * self.diagonal_count * c ** 2 * (r + 1)), + target + self.e.ceil_log2(16 * e0 * c * self.stage_count * size)) if self.logical_layers else max(incoming, 1) + self.bits = minimum if data_bits is None else data_bits + need(type(self.bits) is int and self.bits >= minimum, "funded nonunitary native grid") + self.raw_stages, self.layers = [], [] + self.kinds = Counter() + self.input_height = 1 + self.word_count = 0 + for layer in self.plan: + used, pairs, words, singles, kinds = set(), [], [], {}, [] + for ids, gate in layer: + need(len(ids) == len(gate) in (1, 2) and all(len(row) == len(ids) for row in gate), + "one- or two-wire source gate") + need(len(set(ids)) == len(ids) and not used.intersection(ids) + and all(type(i) is int and 0 <= i < size for i in ids), "disjoint oriented source ownership") + used.update(ids) + self.input_height = max(self.input_height, *(self.e.qheight(p) + for row in gate for z in row for p in (z.re, z.im))) + if len(ids) == 1: + self.scalar(gate[0][0]) + singles[ids[0]] = gate[0][0] + else: + kind, word = self.gate_word(gate) + pairs.append(tuple(ids)) + words.append(word) + kinds.append(kind) + if size == 1: + value = singles.get(0, self.e.C(1)) + self.raw_stages.append(("D", self.e.Vec([value]))) + self.layers.append(dict(begin=len(self.raw_stages) - 1, end=len(self.raw_stages), + pairs=[], kinds=["scalar" if singles else "idle"])) + self.kinds.update(self.layers[-1]["kinds"]) + continue + paired = {i for pair in pairs for i in pair} + free = [i for i in range(size) if i not in paired] + for j in range(0, len(free), 2): + pair = tuple(free[j:j + 2]) + values = tuple(singles.get(i, self.e.C(1)) for i in pair) + for value in values: + self.scalar(value) + word = source.ns["gateword"]("diag", values) + need(source.word_endpoint(word) == self.e.Mat([[values[0], 0], [0, values[1]]]), + "paired scalar endpoint") + pairs.append(pair) + words.append(word) + kinds.append("scalar-pair" if any(i in singles for i in pair) else "idle") + need(len(pairs) == len(words) == size // 2, "complete mixed matching") + permutation = source.ns["matching_permutation"](size, pairs) + self.check_network(permutation, source.ns["benes"](permutation)) + stages = source.ns["local_layer"](size, pairs, words) + need(len(stages) == 16 * m + 5, "source matching-layer stage count") + begin = len(self.raw_stages) + self.raw_stages.extend(stages) + self.layers.append(dict(begin=begin, end=len(self.raw_stages), pairs=pairs, + kinds=kinds, permutation=permutation)) + self.kinds.update(kinds) + self.word_count += len(words) + self.reference_stages = [] + self.setup_height = 1 + for j, (kind, value) in enumerate(self.raw_stages): + if kind == "H": + need(type(value) is int and 0 <= value < m + and j + 1 < len(self.raw_stages) and self.raw_stages[j + 1][0] == "D", + "native bit and following compensation") + self.reference_stages.append((kind, value)) + else: + need(kind == "D" and len(value) == size and all(z != 0 for z in value), + "complete nonzero reference page") + scale = 2 if j and self.raw_stages[j - 1][0] == "H" else 1 + page = self.e.Vec(scale * z for z in value) + need(all(z.norm2() <= (6 * self.gamma) ** 2 + and (1 / z).norm2() <= (200 * self.gamma) ** 2 for z in page), + "native direct/reciprocal coefficient envelope") + self.setup_height = max(self.setup_height, *(self.e.qheight(p) + for z in page for p in (z.re, z.im))) + self.reference_stages.append((kind, page)) + self.setup_allowance = 32 * (self.input_height + self.gamma_bits + 8) + need(self.setup_height <= self.setup_allowance, "reduced rational setup height") + self.stages = [(kind, value if kind == "H" else self.e.Vec(self.e.rounded(z, self.bits) for z in value)) + for kind, value in self.reference_stages] + need(len(self.stages) == self.stage_count + and sum(kind == "D" for kind, _ in self.stages) == self.diagonal_count + and sum(kind == "H" for kind, _ in self.stages) == self.h_count, "full native stage counts") + self.coefficient_width = 0 + u = Q(1, 1 << self.bits) + for (kind, exact), (_, page) in zip(self.reference_stages, self.stages): + if kind == "D": + for a, z in zip(exact, page): + need((a - z).norm2() <= (2 * u) ** 2, "native diagonal approximation") + for part in (z.re, z.im): + integer = part * (1 << self.bits) + need(integer.denominator == 1, "native coefficient grid") + self.coefficient_width = max(self.coefficient_width, abs(integer.numerator).bit_length() + 1) + self.coefficient_bits = self.bits + self.gamma_bits + 5 + need(self.coefficient_width <= self.coefficient_bits, "signed nonunitary coefficient width") + self.coefficient_error = 4 * e0 * self.diagonal_count * c ** 2 * u * r + self.rounding_error = 4 * e0 * c * self.stage_count * size * u + self.native_error = self.coefficient_error + self.rounding_error + self.total_error = self.upstream_error * r + self.native_error + need(self.diagonal_count * c ** 2 * 2 * u < Q(1, 2) + and self.native_error < Q(1, 1 << (target + 1)) + and self.total_error <= Q(1, 1 << target), "complete conditioned target budget") + self.amplitude = 4 * c * (r + 1) + self.temporary_amplitude = self.gamma * self.amplitude + self.integer_bits = self.e.ceil_log2(self.amplitude) + 2 + self.temporary_bits = 2 * self.bits + self.integer_bits + self.gamma_bits + 6 + self.stage_seal = stage_hash(source, self.stages) + self.reference_stage_seal = stage_hash(source, self.reference_stages) + self.contract_seal = self.contract_hash() + + def contract_hash(self): + return budget_hash(self.source, {name: getattr(self, name) for name in self.CONTRACT_FIELDS}) + + def scalar(self, value): + need(value != 0 and value.norm2() <= self.gamma ** 2 + and (1 / value).norm2() <= self.gamma ** 2, "logical scalar and reciprocal admission") + + def gate_word(self, gate): + E, source = self.e, self.source + if gate[0][0] == gate[1][1] == 1 and gate[1][0] == 0: + kind, value = "upper", gate[0][1] + need(value.norm2() <= self.gamma ** 2, "upper shear coefficient admission") + word = source.ns["gateword"](kind, value, self.gamma) + elif gate[0][0] == gate[1][1] == 1 and gate[0][1] == 0: + kind, value = "lower", gate[1][0] + need(value.norm2() <= self.gamma ** 2, "lower shear coefficient admission") + word = source.ns["gateword"](kind, value, self.gamma) + elif gate[0][1] == gate[1][0] == 0: + kind = "diagonal" + values = gate[0][0], gate[1][1] + for value in values: + self.scalar(value) + word = source.ns["gateword"]("diag", values) + else: + kind, word = source.gate_word(gate) + need(len(word) == 13 and sum(op == "H" for op in word) == 6 + and source.word_endpoint(word) == E.Mat(gate), "actual six-H source gate endpoint") + return kind, word + + def execute(self, vector, capture=False): + need(self.contract_hash() == self.contract_seal, "unchanged native precision and magnitude contract") + if self.context is not None: + self.context.validate() + return super().execute(vector, capture) + + def record(self): + need(self.contract_hash() == self.contract_seal, "unchanged recorded native contract") + if self.context is not None: + self.context.validate() + exact_values = dict(prefix=self.prefix, endpoint=self.endpoint, compiled=self.compiled_bound, + radius=self.radius, upstream=self.upstream_error, + coefficient=self.coefficient_error, rounding=self.rounding_error, + native=self.native_error, total=self.total_error, amplitude=self.amplitude) + tau = Q(1, 1 << self.target) + return dict(size=self.size, logical_layers=self.logical_layers, word_count=self.word_count, + kinds=dict(sorted(self.kinds.items())), layer_mapping_sha256=sha256(canonical(self.layers)).hexdigest(), + logical_plan_sha256=self.reference_hash, reference_stages_sha256=self.reference_stage_seal, + stored_stages_sha256=self.stage_seal, tensor_ast_sha256=self.source.tensor_ast, + counts=dict(primitive_stages=self.stage_count, diagonal_stages=self.diagonal_count, + hadamard_stages=self.h_count), + budget=dict(prefix_upper=reported(self.source, self.prefix), + endpoint_upper=reported(self.source, self.endpoint), + gamma_bits=self.gamma_bits, compiled_prefix_upper=reported(self.source, self.compiled_bound), + radius=str(self.radius), target_bits=self.target, incoming_bits=self.incoming, + data_grid_bits=self.bits, coefficient_bits=self.coefficient_bits, + data_bits=self.bits + self.integer_bits, temporary_bits=self.temporary_bits, + coefficient_error_upper=reported(self.source, self.coefficient_error, tau / 2), + rounding_error_upper=reported(self.source, self.rounding_error, tau / 2), + native_error_upper=reported(self.source, self.native_error, tau / 2), + total_error_upper=reported(self.source, self.total_error, tau), + amplitude_upper=reported(self.source, self.amplitude), + exact_budget_sha256=budget_hash(self.source, exact_values)), + setup=dict(input_component_height=self.input_height, reduced_coefficient_height=self.setup_height, + reduced_height_allowance=self.setup_allowance, + observed_stored_coefficient_bits=self.coefficient_width)) + + +class FullRoot: + def __init__(self, source, matrix, radius, target, indices=None): + self.source, self.e = source, source.e + self.matrix = self.e.Mat(matrix) + self.matrix_seal = source.matrix_hash(self.matrix) + self.lam, self.q, self.leading = source.full.admission(self.matrix) + budget = source.full.allocation(len(matrix), self.lam, self.q, radius, target, 0) + self.builder = source.full.Full(budget, source.raw) + self.plan = self.builder.plan(self.matrix, indices) + self.inverse_plan = self.builder.inverse(self.plan) + self.bound = dict(self.builder.bound) + self.budget_seal = budget_hash(source, self.bound) + self.plan_seal, self.inverse_seal = plan_hash(source, self.plan), plan_hash(source, self.inverse_plan) + self.radius, self.target = Q(radius), target + self.indices = list(range(len(matrix))) if indices is None else list(indices) + self.request_seal = self.request_state() + need(source.matrix_hash(self.matrix) == self.matrix_seal, "borrowed exact root preserved during preparation") + + def request_state(self): + return (self.lam, self.q, self.radius, self.target, tuple(self.indices), tuple(self.leading)) + + def validate(self): + need(self.request_state() == self.request_seal + and self.lam == self.bound["lam"] and self.q == self.bound["inverse_bound"] + and self.radius == self.bound["radius"] and self.target == self.bound["target_bits"] + and self.source.matrix_hash(self.matrix) == self.matrix_seal + and budget_hash(self.source, self.builder.bound) == self.budget_seal + and budget_hash(self.source, self.bound) == self.budget_seal + and plan_hash(self.source, self.plan) == self.plan_seal + and plan_hash(self.source, self.inverse_plan) == self.inverse_seal, + "same admitted root, stored plan and full prefix contract") + + def lower(self, inverse=False, incoming=0): + need(type(inverse) is bool, "declared full direction") + self.validate() + defect = self.bound["coefficient_error"] + endpoint = (self.q if inverse else self.lam) + defect + need(endpoint <= self.bound["execution_bound"], "endpoint bound fits the retained full envelope") + return NativeProgram(self.source, len(self.matrix), self.inverse_plan if inverse else self.plan, + self.bound["execution_bound"], endpoint, self.radius, self.target, incoming, + upstream_error=defect, context=self) + + def record(self): + self.validate() + return dict(n=len(self.matrix), indices=self.indices, root_sha256=self.matrix_seal, + full_source_ast_sha256=self.builder.source_ast, + full_budget_sha256=self.budget_seal, + lambda_upper=reported(self.source, self.lam), inverse_upper=reported(self.source, self.q), + leading=[dict(size=n, inverse_upper=reported(self.source, value)) for n, value in self.leading], + stored_prefix_upper=reported(self.source, self.bound["execution_bound"]), + upstream_operator_error_upper=reported(self.source, self.bound["coefficient_error"]), + target_bits=self.target, radius=str(self.radius), + forward_plan_sha256=self.plan_seal, inverse_plan_sha256=self.inverse_seal) diff --git a/research/rational-fft-fixed-h/lowering.py b/research/rational-fft-fixed-h/lowering.py new file mode 100644 index 0000000000..5c2adcfd96 --- /dev/null +++ b/research/rational-fft-fixed-h/lowering.py @@ -0,0 +1,362 @@ +"""Source-bound rational FFT lowering to H/2 and rounded diagonal stages. + +Jose Manuel (jmbarrancoml), with substantial OpenAI Codex assistance. +Apache-2.0. Fixed words, Beneš routing and FFT schedule: Vlado Vojdanovski. +Rational phases: Chafik Boukhalfa and the retained PR628 lineage. +""" +import ast +import copy +from fractions import Fraction as Q +from hashlib import sha256 +import json +from pathlib import Path +import sys +from types import ModuleType, SimpleNamespace + +HERE = Path(__file__).resolve().parent +ROOT = HERE.parents[1] +PREFIX = "research/rational-fft-fixed-h/" +EXACT = PREFIX + "vendor/rational_fft/exact.py" +ENGINE = PREFIX + "vendor/rational_fft/engine.py" +PROVENANCE = PREFIX + "vendor/PROVENANCE.json" +CIRCUIT = "research/capped-phase-fusion/evidence/inplace_algebra_ultra/check_circuits.py" +TENSOR = "research/capped-phase-fusion/evidence/inplace_algebra_ultra/check_tensor_compiler.py" +ORACLE = "research/cubic-cyclic-inverse/vendor/interval_kernel.py" +DENSE_PREFIX = PREFIX + "vendor/full_source/" +DENSE_CONTEXT = {"research/unconditional-batched-crt/vendor/GAUGEFREE-INTERIORS.md", + "research/unconditional-batched-crt/vendor/SHALLOW-TOEPLITZ.md"} +K = 1 << 14 + + +def need(ok, message): + if not ok: + raise ValueError(message) + + +def canonical(value): + return json.dumps(value, sort_keys=True, separators=(",", ":")).encode() + + +def norm2(vector): + return sum((z.norm2() for z in vector), Q(0)) + + +def authenticate(overrides=None): + rows = json.loads((HERE / "SOURCES.json").read_text())["files"] + raw, seals = {}, [] + for row in rows: + if row["usage"] != "input": + continue + path = row["path"] + value = (overrides or {}).get(path) + if value is None: + value = (ROOT / path).read_bytes() + need(sha256(value).hexdigest() == row["sha256"], "source SHA256: " + path) + raw[path] = value + seals.append(dict(path=path, sha256=row["sha256"])) + provenance = json.loads(raw[PROVENANCE]) + for row in provenance["files"]: + need(sha256(raw[row["retained_path"]]).hexdigest() == row["sha256"], + "unchanged published FFT producer") + return raw, seals + + +class Sources: + def __init__(self, overrides=None, include_dense=False): + need(type(include_dense) is bool, "declared source inventory") + self.raw, all_seals = authenticate(overrides) + self.seals = [row for row in all_seals if include_dense + or not (row["path"].startswith(DENSE_PREFIX) or row["path"] in DENSE_CONTEXT)] + self.e = ModuleType("_fixed_h_retained_exact") + exec(compile(self.raw[EXACT], "", "exec", optimize=0), + self.e.__dict__) + self.engine = ModuleType("_fixed_h_retained_engine") + previous = sys.modules.get("exact") + sys.modules["exact"] = self.e + try: + exec(compile(self.raw[ENGINE], "", "exec", optimize=0), + self.engine.__dict__) + finally: + if previous is None: + del sys.modules["exact"] + else: + sys.modules["exact"] = previous + names = {"benes", "D", "u3", "id2", "id3", "canonical", "gateword", + "matching_permutation", "local_layer"} + tree = ast.parse(self.raw[TENSOR]) + tree.body = [copy.deepcopy(node) for node in tree.body + if isinstance(node, ast.FunctionDef) and node.name in names] + need({node.name for node in tree.body} == names, "complete fixed-word compiler closure") + self.tensor_ast = sha256(ast.dump(tree).encode()).hexdigest() + C = self.e.C + scalar = SimpleNamespace(sympify=C.wrap, simplify=C.wrap, Integer=lambda n: C(n), + Rational=lambda a, b=1: C(Q(a, b)), I=C(0, 1)) + self.ns = dict(s=scalar) + exec(compile(tree, "", "exec", optimize=0), self.ns) + + def matrix_hash(self, matrix): + return sha256((str((len(matrix), len(matrix[0]) if matrix else 0)) + + self.e.vector_hash(x for row in matrix for x in row)).encode()).hexdigest() + + def word_endpoint(self, word): + C, Mat = self.e.C, self.e.Mat + matrix = Mat([[1, 0], [0, 1]]) + for operation in word: + left = Mat([[1, 1], [1, -1]]) if operation == "H" else Mat( + [[operation[0], 0], [0, operation[1]]]) + matrix = Mat([[sum((left[i][k] * matrix[k][j] for k in range(2)), C()) + for j in range(2)] for i in range(2)]) + return matrix + + def butterfly_word(self, a, b, phase, adjoint=False): + C = self.e.C + need(type(a) is Q and type(b) is Q and a >= Q(1, 2) and b >= Q(1, 2) + and a * a + b * b == 1 and phase.norm2() == 1, "bounded rational unit splitter") + r, t = C(b / a), phase * (b / a) + word = self.ns["canonical"](self.ns["u3"](t) + [self.ns["D"](a, -phase / a)] + + list(reversed(self.ns["u3"](r)))) + if adjoint: + word = [op if op == "H" else tuple(z.conjugate() for z in op) + for op in reversed(word)] + return word + + def gate_word(self, gate): + Mat = self.e.Mat + if gate == Mat([[0, 1], [1, 0]]): + kind, word = "swap", self.ns["gateword"]("swap") + elif gate[0][0].im == gate[1][0].im == 0 and min(gate[0][0].re, gate[1][0].re) >= Q(1, 2): + a, b = gate[0][0].re, gate[1][0].re + phase = gate[0][1] / b + kind, word = "butterfly", self.butterfly_word(a, b, phase) + else: + need(gate[0][0].im == gate[0][1].im == 0, "stored FFT adjoint orientation") + a, b = gate[0][0].re, gate[0][1].re + phase = (gate[1][0] / b).conjugate() if b != 0 else self.e.C() + kind, word = "adjoint", self.butterfly_word(a, b, phase, True) + need(len(word) == 13 and sum(op == "H" for op in word) == 6, + "six-H seven-diagonal shape") + need(self.word_endpoint(word) == gate, "exact endpoint equals the actual stored gate") + return kind, word + + +def stage_hash(source, stages): + rows = [(kind, value if kind == "H" else source.e.vector_hash(value)) for kind, value in stages] + return sha256(canonical(rows)).hexdigest() + + +def plan_hash(source, plan): + return sha256(canonical([[(list(ids), source.matrix_hash(gate)) for ids, gate in layer] + for layer in plan])).hexdigest() + + +class FFT: + def __init__(self, source, size, indices, radius, target, incoming=0, *, adjoint=False, + phase_bits=None, data_bits=None): + need(type(size) is int and size > 0 and size & (size - 1) == 0, "dyadic complete field") + need(indices and all(type(i) is int and 0 <= i < size for i in indices) + and len(set(indices)) == len(indices), "distinct oriented integer FFT wires") + n = len(indices) + need(n & (n - 1) == 0, "dyadic active FFT size") + need(type(radius) in (int, Q) and radius >= 1 and type(target) is int and target >= 1 + and type(incoming) is int and incoming >= 0 and type(adjoint) is bool, + "exact input radius, target, grid and direction") + self.source, self.e = source, source.e + self.size, self.indices, self.radius = size, tuple(indices), Q(radius) + self.target, self.incoming, self.adjoint = target, incoming, adjoint + d, field_d = n.bit_length() - 1, size.bit_length() - 1 + self.original_layers = d + 1 if d else 0 + self.stage_count = self.original_layers * (16 * field_d + 5) + self.diagonal_count = self.original_layers * (8 * field_d + 3) + self.h_count = self.original_layers * (8 * field_d + 2) + required_phase = max(3, target + self.e.ceil_log2(8 * d * (radius + 1))) if d else 3 + self.phase_bits = required_phase if phase_bits is None else phase_bits + need(type(self.phase_bits) is int and self.phase_bits >= required_phase, + "funded phase precision and splitter chart") + required_data = max(incoming, target + self.e.ceil_log2(16 * self.diagonal_count * K ** 2 * (radius + 1)), + target + self.e.ceil_log2(16 * K * self.stage_count * size)) if d else max(incoming, 1) + self.bits = required_data if data_bits is None else data_bits + need(type(self.bits) is int and self.bits >= required_data, "funded native coefficient/data grid") + self.phases = source.engine.Phases(self.phase_bits, source.raw[ORACLE]) + self.builder = source.engine.Builder(self.phases, self.bits, source.raw[CIRCUIT]) + forward = self.builder.namespace["fft"](list(indices)) + self.plan = self.builder.namespace["inverse"](forward) if adjoint else forward + need(len(self.plan) == self.original_layers, "retained FFT matching-layer count") + self.reference_hash = plan_hash(source, self.plan) + self.raw_stages, self.layers, self.words = [], [], [] + for layer in self.plan: + pairs, words, kinds, used = [], [], [], set() + for ids, gate in layer: + need(len(ids) == 2 and len(set(ids)) == 2 and not used.intersection(ids) + and all(type(i) is int and 0 <= i < size for i in ids), "matching wire ownership") + need(self.e.unitary(gate), "unitary completed FFT gate") + kind, word = source.gate_word(gate) + used.update(ids) + pairs.append(tuple(ids)) + words.append(word) + kinds.append(kind) + free = [i for i in range(size) if i not in used] + for j in range(0, len(free), 2): + pairs.append(tuple(free[j:j + 2])) + words.append(source.ns["gateword"]("diag", (1, 1))) + kinds.append("idle") + need(len(pairs) == len(words) == size // 2, "complete matching and inactive words") + permutation = source.ns["matching_permutation"](size, pairs) + network = source.ns["benes"](permutation) + self.check_network(permutation, network) + stages = source.ns["local_layer"](size, pairs, words) + need(len(stages) == 16 * field_d + 5, "retained source local-stage count") + begin = len(self.raw_stages) + self.raw_stages.extend(stages) + self.layers.append(dict(begin=begin, end=len(self.raw_stages), pairs=pairs, kinds=kinds, + permutation=permutation, route_layers=len(network))) + self.words.append(words) + self.reference_stages = [] + for i, (kind, value) in enumerate(self.raw_stages): + if kind == "H": + need(type(value) is int and 0 <= value < field_d + and i + 1 < len(self.raw_stages) and self.raw_stages[i + 1][0] == "D", + "Hadamard bit and immediate compensation stage") + self.reference_stages.append(("H", value)) + else: + need(kind == "D" and len(value) == size and all(z != 0 for z in value), + "complete nonzero coefficient table") + multiplier = 2 if i and self.raw_stages[i - 1][0] == "H" else 1 + page = self.e.Vec(multiplier * z for z in value) + need(all(z.norm2() <= 36 for z in page), "bounded native coefficient modulus") + need(all(self.e.qheight(part) <= 8 * self.phase_bits + 64 + for z in page for part in (z.re, z.im)), "rational coefficient setup height") + self.reference_stages.append(("D", page)) + self.stages = [(kind, value if kind == "H" else self.e.Vec(self.e.rounded(z, self.bits) for z in value)) + for kind, value in self.reference_stages] + self.coefficient_width = 0 + for kind, page in self.stages: + if kind == "D": + for z in page: + for part in (z.re, z.im): + integer = part * (1 << self.bits) + need(integer.denominator == 1, "stored coefficient grid") + self.coefficient_width = max(self.coefficient_width, abs(integer.numerator).bit_length() + 1) + need(self.coefficient_width <= self.bits + 5, "signed coefficient record width") + need(len(self.stages) == self.stage_count + and sum(kind == "D" for kind, _ in self.stages) == self.diagonal_count + and sum(kind == "H" for kind, _ in self.stages) == self.h_count, "complete native counts") + u = Q(1, 1 << self.bits) + self.phase_error = 2 * d * Q(1, 1 << self.phase_bits) * self.radius + self.coefficient_error = 4 * self.diagonal_count * K ** 2 * u * self.radius + self.rounding_error = 4 * K * self.stage_count * size * u + self.native_error = self.coefficient_error + self.rounding_error + self.total_error = self.phase_error + self.native_error + need(self.diagonal_count * K ** 2 * 2 * u < Q(1, 2) + and self.total_error < 3 * Q(1, 1 << target) / 4, "complete FFT precision budget") + self.amplitude = 4 * K * (self.radius + 1) + self.integer_bits = self.e.ceil_log2(self.amplitude) + 2 + self.temporary_bits = 2 * self.bits + self.integer_bits + 6 + self.stage_seal = stage_hash(source, self.stages) + self.reference_stage_seal = stage_hash(source, self.reference_stages) + for (kind, exact), (_, stored) in zip(self.reference_stages, self.stages): + if kind == "D": + need(all((a - b).norm2() <= (2 * u) ** 2 for a, b in zip(exact, stored)), + "quantized diagonal operator error") + + def check_network(self, permutation, network): + data = list(range(self.size)) + for layer in network: + used, bits = set(), set() + for i, j, flip in layer: + need(type(flip) is bool and i not in used and j not in used + and 0 <= i < self.size and 0 <= j < self.size, "source switch ownership") + mask = i ^ j + need(mask > 0 and mask & (mask - 1) == 0, "one-bit switch") + bits.add(mask) + used.update((i, j)) + if flip: + data[i], data[j] = data[j], data[i] + need(len(bits) == 1 and len(used) == self.size, "complete common-bit routing layer") + need(all(data[permutation[i]] == i for i in range(self.size)), "source routing orientation") + + def execute(self, vector, capture=False): + E, scale = self.e, 1 << self.bits + state = E.Vec(vector) + need(len(state) == self.size and norm2(state) <= self.radius ** 2, + "declared FFT input dimension and Euclidean radius") + need(all((part * (1 << self.incoming)).denominator == 1 + for z in state for part in (z.re, z.im)), "declared incoming data grid") + need(plan_hash(self.source, self.plan) == self.reference_hash + and stage_hash(self.source, self.reference_stages) == self.reference_stage_seal + and stage_hash(self.source, self.stages) == self.stage_seal, + "unchanged stored FFT and prepared stages") + borrowed = E.vector_hash(state) + stats = dict(real_products=0, complex_roundings=0, h_pair_calls=0, + maximum_data_numerator_bits=0, maximum_temporary_numerator_bits=0, + maximum_temporary_fractional_bits=0) + traces = [] + + def observe(value): + integer = value * (1 << (2 * self.bits)) + need(integer.denominator == 1, "exact native arithmetic temporary grid") + width = abs(integer.numerator).bit_length() + 1 + stats["maximum_temporary_numerator_bits"] = max(stats["maximum_temporary_numerator_bits"], width) + stats["maximum_temporary_fractional_bits"] = max( + stats["maximum_temporary_fractional_bits"], value.denominator.bit_length() - 1) + need(width <= self.temporary_bits + and abs(value) <= 16 * getattr(self, "temporary_amplitude", self.amplitude), + "native arithmetic temporary allowance") + return value + + def store(value): + result = E.rounded(value, self.bits) + stats["complex_roundings"] += 1 + for part in (result.re, result.im): + integer = part * scale + need(integer.denominator == 1, "stored native common grid") + width = abs(integer.numerator).bit_length() + 1 + stats["maximum_data_numerator_bits"] = max(stats["maximum_data_numerator_bits"], width) + need(width <= self.bits + self.integer_bits, "signed data record width") + return result + + need(norm2(state) <= self.amplitude ** 2, "initial magnitude envelope") + for number, (kind, value) in enumerate(self.stages): + if kind == "D": + output = [] + for a, x in zip(value, state): + ar, ai = observe(a.re * x.re), observe(a.re * x.im) + br, bi = observe(a.im * x.re), observe(a.im * x.im) + stats["real_products"] += 4 + output.append(store(E.C(observe(ar - bi), observe(ai + br)))) + state = E.Vec(output) + else: + step = 1 << value + for i in range(self.size): + if not i & step: + a, b = state[i], state[i + step] + sr, si = observe(a.re + b.re), observe(a.im + b.im) + dr, di = observe(a.re - b.re), observe(a.im - b.im) + state[i] = store(E.C(observe(sr / 2), observe(si / 2))) + state[i + step] = store(E.C(observe(dr / 2), observe(di / 2))) + stats["h_pair_calls"] += 1 + need(norm2(state) <= self.amplitude ** 2, "actual native prefix magnitude") + if capture: + traces.append(dict(stage=number, kind=kind, output_sha256=E.vector_hash(state))) + need(borrowed == E.vector_hash(E.Vec(vector)), "borrowed input preserved") + need(stats["complex_roundings"] == self.stage_count * self.size + and stats["h_pair_calls"] == self.h_count * self.size // 2, + "all active and inactive rounding sites charged") + return state, stats, traces + + def record(self): + return dict(field_size=self.size, active_indices=list(self.indices), adjoint=self.adjoint, + radius=str(self.radius), target_bits=self.target, incoming_bits=self.incoming, + phase_bits=self.phase_bits, data_bits=self.bits, + observed_coefficient_bits=self.coefficient_width, + counts=dict(logical_layers=self.original_layers, primitive_stages=self.stage_count, + diagonal_stages=self.diagonal_count, hadamard_stages=self.h_count), + layers=self.layers, logical_plan_sha256=self.reference_hash, + reference_stages_sha256=self.reference_stage_seal, stored_stages_sha256=self.stage_seal, + source_original_ast_sha256=self.builder.original_ast, + source_adapted_ast_sha256=self.builder.adapted_ast, + tensor_ast_sha256=self.source.tensor_ast, phases=self.phases.record(), + budget=dict(prefix_bound=K, phase_error=str(self.phase_error), + coefficient_error=str(self.coefficient_error), rounding_error=str(self.rounding_error), + total_error=str(self.total_error), amplitude_bound=str(self.amplitude), + coefficient_bits=self.bits + 5, data_bits=self.bits + self.integer_bits, + temporary_bits=self.temporary_bits)) diff --git a/research/rational-fft-fixed-h/vendor/PROVENANCE.json b/research/rational-fft-fixed-h/vendor/PROVENANCE.json new file mode 100644 index 0000000000..22cf8bdcb9 --- /dev/null +++ b/research/rational-fft-fixed-h/vendor/PROVENANCE.json @@ -0,0 +1,21 @@ +{ + "schema_version": 1, + "original_pull_request": "https://github.com/CrocSwap/integer-mult-bounds/pull/4205", + "original_commit": "78c72e563173f3822fcb4efed9ceeccf05630258", + "files": [ + { + "original_path": "research/rational-unitary-fft/engine.py", + "retained_path": "research/rational-fft-fixed-h/vendor/rational_fft/engine.py", + "git_blob": "e9f2f2b08dc2bc4c646e47dca92059351ca89384", + "sha256": "9e36aee0835e6b1b62997433581f88b5a434ccab8ae43e5c25d1deaaf1ebf653" + }, + { + "original_path": "research/rational-unitary-fft/exact.py", + "retained_path": "research/rational-fft-fixed-h/vendor/rational_fft/exact.py", + "git_blob": "aae41e6acfb2e5a4fafdd85727891eb746b57b42", + "sha256": "c9b89d7357210af69865690ba12e695366961f8b501a8fab9d0fefe3fb30c479" + } + ], + "attribution": "Jose Manuel (jmbarrancoml), with substantial OpenAI Codex assistance. Phase projection follows Chafik Boukhalfa/PR628; the FFT/shear schedule follows Vlado Vojdanovski/PR1673. Classical Fourier/FFT methods and the retained lineage remain credited. Apache-2.0.", + "scope": "Unchanged published exact arithmetic and rational FFT producer. The new loader authenticates their bytes and declared dependencies before executing them. Retention does not constitute independent acceptance of the predecessor's proofs." +} diff --git a/research/rational-fft-fixed-h/vendor/full_source/FULL-PROOF.md b/research/rational-fft-fixed-h/vendor/full_source/FULL-PROOF.md new file mode 100644 index 0000000000..e650f46863 --- /dev/null +++ b/research/rational-fft-fixed-h/vendor/full_source/FULL-PROOF.md @@ -0,0 +1,374 @@ +# Full Toeplitz recursion with a polynomial conditioning envelope + +Jose Manuel (`jmbarrancoml`), with substantial OpenAI Codex assistance. +Apache-2.0. + +## 1. Domain, result and predecessor scope + +Let \(N=2^d\), and let \(T\in\mathbb Q(i)^{N\times N}\) be Toeplitz. +Write \(T_b=T[0:b,0:b]\). I assume +\[ +\Lambda,Q\ge1,\qquad \|T\|_2\le\Lambda,\qquad +\|T_b^{-1}\|_2\le Q\quad(1\le b\le N). +\tag{1} +\] +The full inverse bound alone would not suffice: the Toeplitz matrix +\(\left[\begin{smallmatrix}0&1\\1&0\end{smallmatrix}\right]\) is invertible +but has a singular first section. I include it as an admission control. +I allow indefinite and non-Hermitian matrices satisfying (1). + +I execute the retained `full_plan` on exactly \(N\) logical complex wires, +using the completed generator circuit and the same stored gates in reverse +for inversion. For input norm at most \(R\ge1\), incoming grid +\(F_0\in\mathbb Z_{\ge0}\) and target \(2^{-t}\), \(t\in\mathbb Z_{\ge1}\), +I prove forward and inverse error at most \(2^{-t}\) +with online widths +\[ +O(F_0+t+\log N+\log\Lambda+\log Q+\log R). +\tag{2} +\] +The circuit has \(O(\log^4N)\) layers and \(O(N\log^4N)\) gates. +I give explicit constants below. + +Vlado/PR1673 supplied the recursive compiler and displacement identities. +The retained `GAUGEFREE-INTERIORS.md`, Theorem 4.1, supplies a +conditioning-dependent guard of +\(O(\log N\log(2+\Lambda/\mu))\) for its accretive domain. +I use the scalar unit-LDU factors of legal completed-call sets to obtain +a joint polynomial envelope in dimension and conditioning. This +sharpens the sufficient guard's conditioning dependence under (1). +I leave the selected multiplication exponent unchanged. + +ToikY's [PR4209](https://github.com/CrocSwap/integer-mult-bounds/pull/4209) +at `ccc33bfde8baadef7979fac81e28c0a6f6f278c4` supplies sharper local +Schur closure under a row/column near-identity majorant. I retain that +distinct credit. I do not infer its hypothesis from (1) or use its +sharper constants in the present allocation. + +### Executable admission + +`admission()` computes exact inverses of the same root's leading sections. +It returns \(\lambda_c,Q_c\), each at least one, using the maximum row +and column sums of \(|\Re z|+|\Im z|\). These bound the corresponding +operator norms and also support the executable component-envelope guards. +The verifier binds them to the unchanged root passed to `Full.plan()`. +That lower-level factory does not repeat the complete admission. + +If (1) is supplied as operator-norm information, then +\[ +\lambda_c\le\sqrt{2N}\Lambda,\qquad +Q_c\le\sqrt{2N}Q. +\tag{3} +\] +I use the executable envelopes in its explicit constants. The inflation +in (3) preserves (2). The written operator lemma below uses any valid +\(\Lambda,Q\); passing uninflated operator bounds directly to the +component-guarded factory need not satisfy its stronger guards. + +Exact admission and Schur preparation are offline rational computations. +I do not equate their denominators with online data words or claim a +complete physical preparation/delivery tariff. + +## 2. Contiguous descendants, shears and rank + +The node on logical interval \([a,b)\) has reference matrix +\[ +M_{a,b}=S(a,b)=T_b/T_a, +\] +with the empty prefix convention at \(a=0\). Its children are \(S(a,c)\) +and \(S(c,b)\), where \(c=(a+b)/2\), by the Schur quotient identity. +Consequently +\[ +\|M_{a,b}\|\le U:=\Lambda+\Lambda^2Q,\qquad +\|M_{a,b}^{-1}\|\le Q. +\tag{4} +\] +The inverse is a principal block of \(T_b^{-1}\); the direct bound follows +from the Schur formula and restrictions of the original \(T\). + +At a split, write \(M=\left[\begin{smallmatrix}A&B\\C&D\end{smallmatrix}\right]\), +\(K=A^{-1}B\), \(L=CA^{-1}\). Let \(R_a\) select the coordinates +\([a,c)\) from \([0,c)\). The block inverse identities give +\[ +K=R_aT_c^{-1}T[0:c,c:b],\qquad +L=T[c:b,0:c]T_c^{-1}R_a^T. +\] +Thus, with \(\chi=\Lambda Q\), +\[ +\|K\|,\|L\|\le\chi. +\tag{5} +\] +The component-envelope versions follow under the executable admission. +The source computes \(L\) as `solve(A.T, C.T).T`. I preserve ordinary +transpose in this complex identity. + +The retained `SHALLOW-TOEPLITZ.md`, Sections 4–6, gives commutator rank +at most \(2+\text{depth}\) for the Schur nodes. Its rectangular identities +bound the shear commutator ranks by \(2\,\text{depth}+5\) and +\(2\,\text{depth}+7\). Conversion to +\(E_X=X-ZXZ^T\) adds at most one rank. I use the uniform allowance +\[ +r_*=2d+8. +\tag{6} +\] +The exact rational producer retains that rank; I do not use the numerical +source's floating SVD threshold as a universal rank argument. +The executable also checks exact generator rank against separate +elimination in its finite fixtures. + +## 3. Uniform scalar-LDU bounds + +Write \(M_J=\ell_J\delta_Ju_J\) for the scalar unit-LDU factors of a +contiguous descendant. All its leading sections are further \(S(a,c)\), +so (4) also applies during scalar elimination. +For \(A_k=M_J[0:k,0:k]\), the last row of +\(A_k^{-1}M_J[0:k,k:]\) is the corresponding unit-upper row tail. +The last column of \(M_J[k:,0:k]A_k^{-1}\) is the unit-lower column tail. +Each therefore has norm at most +\[ +UQ=\chi+\chi^2\le2\chi^2. +\] +The Frobenius estimate therefore permits +\[ +\|\ell_J\|,\|u_J\|\le W:=3N\chi^2. +\] +The pivot diagonal has norm at most \(U\) and inverse norm at most \(Q\). +From +\[ +u_J^{-1}=M_J^{-1}\ell_J\delta_J,\qquad +\ell_J^{-1}=\delta_Ju_JM_J^{-1} +\] +I obtain +\[ +\|\ell_J^{-1}\|,\|u_J^{-1}\|\le V:=6N\chi^4, +\tag{7} +\] +using \(QU\le2\chi^2\). A diagonal containing some completed pivots +and otherwise ones has direct and inverse norms at most \(2\chi^2\). + +## 4. Exact downsets and the upper frontier + +At an internal node \(v\), denote the upper and lower atomic calls by +\[ +R_v=\begin{pmatrix}I&K_v\\0&I\end{pmatrix},\qquad +B_v=\begin{pmatrix}I&0\\L_v&I\end{pmatrix}. +\] +The source executes \(R_v\), its two child programs in parallel, then +\(B_v\). A leaf executes its scalar pivot. +The atomic partial order places \(R_v\) before both children and both +children before \(B_v\). Incomparable calls act on disjoint intervals +and commute. I count identity shears as calls even when they emit zero layers. + +For a downset \(I\), let \(Q(I)\) be its product in any legal linearization. +At a node, the restricted downset is empty; or it contains \(R_v\) and +child downsets \(I_0,I_1\) but not \(B_v\); or it contains the whole node. +For an empty node set all three factors to identity. At a completed leaf +set them to \(I,d_i,I\). Otherwise define +\[ +\begin{aligned} +\mathcal L(I)&=B_v^\epsilon + \operatorname{diag}(\mathcal L(I_0),\mathcal L(I_1)),\\ +\mathcal D(I)&=\operatorname{diag}(\mathcal D(I_0),\mathcal D(I_1)),\\ +\mathcal U(I)&=\operatorname{diag}(\mathcal U(I_0),\mathcal U(I_1))R_v, +\end{aligned} +\tag{8} +\] +where \(\epsilon=1\) exactly when the lower call has completed. +Then both children are complete. Substituting the child identities gives +\[ +\boxed{Q(I)=\mathcal L(I)\mathcal D(I)\mathcal U(I).} +\tag{9} +\] +This induction never commutes overlapping factors. + +Uniqueness of scalar unit-LDU gives +\[ +\ell_v=B_v\operatorname{diag}(\ell_{v0},\ell_{v1}),\quad +\delta_v=\operatorname{diag}(\delta_{v0},\delta_{v1}),\quad +u_v=\operatorname{diag}(u_{v0},u_{v1})R_v. +\] +Hence \(\mathcal L(I)\) is the direct sum of full lower factors on the +maximal completed subtrees, with identity elsewhere. The partial diagonal +contains completed pivots and otherwise ones. + +For the upper frontier, stop at a node whose upper call has not completed; +otherwise descend to its two children. Retain singleton leaves. +The resulting intervals partition the coordinates. With +\(U_{\mathcal F}=\bigoplus_{v\in\mathcal F}u_v\), induction gives +\[ +u_{\rm root}=U_{\mathcal F}\mathcal U(I),\qquad +\mathcal U(I)=U_{\mathcal F}^{-1}u_{\rm root}. +\tag{10} +\] +At an unvisited node the identity is \(u_v=u_vI\); at a visited node, +substitute the two child identities and multiply by \(R_v\). + +Equations (7)–(10) bound the lower factor and its inverse by \(V\), +the partial diagonal and its inverse by \(2\chi^2\), and the upper +factor and its inverse by \(VW\). Thus +\[ +\boxed{\|Q(I)\|,\|Q(I)^{-1}\| +\le216N^3\chi^{12}<\mathcal B:=256N^3\chi^{12}.} +\tag{11} +\] +I use a fixed number of polynomial factors, rather than multiplying +conditioning factors through the recursion depth. + +## 5. Approximate atomic calls and physical cuts + +There are \(J=3N-2\) atomic calls: two shears per internal node and one +scalar per leaf. Approximate each by operator error at most \(\eta\). +For any legal linearization, set +\[ +G_j=Q_j^{-1}(\widehat A_j-A_j)Q_{j-1},\qquad +\theta=J\mathcal B^2\eta. +\] +The downset bound (11) applies to its prefixes. +If \(\theta\le1/4\), the product of \(I+G_j\), or its inverse, has norm +below two and distance from identity at most \(2\theta\). +Using \(\Lambda,Q\le\mathcal B\), I obtain +\[ +\|\widehat T-T\|,\|\widehat T^{-1}-T^{-1}\| +\le e:=2J\mathcal B^3\eta. +\tag{12} +\] +Completed approximate prefixes and inverse prefixes have norm at most +\(2\mathcal B\). Remaining suffixes and inverse suffixes have norm at most +\(2L_0\mathcal B\), where \(L_0=\max(1,\Lambda,Q)\). +I apply the exact LDU identities to the ideal calls, not to the +approximate endpoints. + +At a physical gate cut, completed calls form a downset \(I\). +Started but unfinished calls form a pairwise incomparable set \(A\): +the source finishes an upper call before starting its descendants and +finishes those descendants before its lower call. Active supports are +therefore disjoint, even when sibling ranks and lengths differ. + +Write each active stored endpoint as \(\widehat A_a=T_aF_a\), with executed +prefix \(F_a\) and remaining suffix \(T_a\). A completed call executed +after an active call began must act on a disjoint subtree, so I may commute +only that disjoint work. I obtain the exact identities +\[ +C_{\rm prefix}=(\bigoplus_{a\in A}F_a\oplus I)\widehat Q(I),\qquad +C_{\rm suffix}=\widehat Q_{\rm remaining} + (\bigoplus_{a\in A}T_a\oplus I). +\tag{13} +\] +The set \(I\cup A\) is also a downset, so I may order \(I\), then \(A\), +then the remaining calls in the suffix argument. +A direct sum contributes the largest fragment norm. I do not require +aligned child macro boundaries. A zero-layer identity completes at its +position in the source schedule and contributes no active fragment. + +## 6. Child, scalar and data bounds + +The component envelope used by a shear's exact producer satisfies +\(B_E\le4N\chi\): its displacement norm is at most \(2\chi\), and a +component row/column sum is at most \(\sqrt{2m}\) times that norm. +Substituting \(m\le N\), \(r\le r_*\), \(B_E\le4N\chi\) in +[GENERATOR-PROOF.md](GENERATOR-PROOF.md), put +\[ +\alpha_*=8r_*N^2\chi,\quad K_{G,*}=2N+64N^2\chi, +\] +\[ +H_*=\max(3200N^3,2+2K_{G,*})(10+4\alpha_*+2K_{G,*}). +\] +The generator's radius-one allocation supplies an unrounded operator +defect. I discard its local data grid and execute the flattened full +program on one global grid. + +A scalar pivot has \(1/Q\le|a|\le U\). I round it at \(s\) bits, +with \(\eta=2^{-s}\le1/(2Q)\), and invert the same stored value. +Then \(|\widehat a|\ge1/(2Q)\), its direct norm is at most \(U+1\), +and its inverse norm is at most \(2Q\). +Consequently +\[ +\mathcal H=\max(H_*,U+1,2Q),\qquad +K_{\rm phys}=2L_0\mathcal B\mathcal H +\tag{14} +\] +bounds physical prefixes, suffixes and their inverses. +An arbitrary segment needs the conservative bound \(K_{\rm phys}^2\). + +The unchanged source schedule permits +\[ +D_*=1+2dr_*(8d^2+32d+18),\qquad +G_*=N[1+dr_*(4d^2+24d+22)]. +\tag{15} +\] +I sum two generator calls per level and the \(N\) scalar leaves; parallel +siblings contribute their maximum depth. +Choose +\[ +s=t+\lceil\log_2(4J\mathcal B^3(R+1))\rceil,\quad +F=\max(F_0,t+\lceil\log_2(2D_*K_{\rm phys}N)\rceil). +\tag{16} +\] +These choices ensure the perturbation and scalar margins above. +One layer's rounding error is at most \(N2^{-F}\). Thus either direction +has error at most +\[ +eR+D_*K_{\rm phys}N2^{-F}\le2^{-t}, +\tag{17} +\] +and intermediate amplitude at most +\[ +A=K_{\rm phys}R+D_*K_{\rm phys}^2N2^{-F}. +\tag{18} +\] +The inverse direction uses the exact stored inverse program on its own +input satisfying the same \(R,F_0\) contract. I claim no exact rounded +forward/inverse round trip. + +If \(\widehat a=(a+ib)/2^s\) is a scalar leaf, its reciprocal +\(2^s(a-ib)/(a^2+b^2)\) has component height at most +\(2s+2\lceil\log_2(U+2)\rceil+8\). +I take the maximum of this bound and the generator children's declared +coefficient bounds, all fixed before online execution. +With \(I=\max(1,\lceil\log_2\max(1,A)\rceil)+2\), signed grid words fit +\(F+I\) bits and constant-size gate temporaries fit +\[ +2F+8C_{\rm bits}+2I+32. +\] +The constants in (6), (11), (14) and (15) are polynomial in \(N,\Lambda,Q\). +The child precision formulas and (16) therefore give (2), including the +executable envelope inflation (3). + +## 7. Source binding, finite evidence and limits + +I authenticate and execute the unmodified `full_plan`, `parallel` and +`inverse` ASTs. I supply exact rational matrix slicing, ordinary transpose, +solve, subtraction and multiplication; wrap recursive calls only to +record their exact reference matrices; and bind `lowdisp_shear` to the +completed generator program. I preserve source assertions in its +dependencies under optimized Python. + +The verifier checks direct full matrices and separate inverse references, +exact stored inversion, leading-section bounds, scalar-LDU identities, +legal downsets and upper-frontier identities. It enumerates all downsets +through size four and checks source-timed downsets on larger fixtures. +Selected physical cuts compare both identities in (13), including +simultaneous nonempty active fragments. Timing checks explicitly cover +unequal sibling lengths and zero-layer calls. + +I test indefinite and half-shift Gaussian matrices, complex inputs, +small and large scalar pivots, reversed wire labels, a nonzero 512-bit +incoming fixture and a size-16 identity with zero-layer shear calls. +The finite Frobenius inequalities include the dimension factor \(N\); +they corroborate consequences of the written operator proof rather than +independently establish its sharper bounds. +Exact dense matrices remain checker oracles outside the online +representation claim. Observed heights describe completed rational +components rather than exhaustive `Fraction` internals. + +I preserve the earlier FFT, triangular and generator receipts and the +fixed-head scope of James Reed's finite budget/source-seal audit at +`fdcdc269`. That audit excludes gate replay and universal-proof acceptance. +I credit Vlado Vojdanovski/PR1673 and the retained gauge-free source lineage +for the original recursive, displacement and conditioning arguments; +Chafik Boukhalfa/PR628 for certified phases; Cinna Davis for the later caller; +ToikY for the distinct same-rho Schur theorem; and classical LU/Schur methods. +Internal agents belong to the author process. + +Complete Gaussian assembly, physical coefficient setup/delivery, fixed-tape +costs and independent mathematical acceptance remain open. diff --git a/research/rational-fft-fixed-h/vendor/full_source/GENERATOR-PROOF.md b/research/rational-fft-fixed-h/vendor/full_source/GENERATOR-PROOF.md new file mode 100644 index 0000000000..8bef383eec --- /dev/null +++ b/research/rational-fft-fixed-h/vendor/full_source/GENERATOR-PROOF.md @@ -0,0 +1,277 @@ +# Complete rounded bounded-generator circuits + +Jose Manuel (`jmbarrancoml`), with substantial OpenAI Codex assistance. +Apache-2.0. + +## 1. Input, output and inherited generator identity + +I consider a rational complex matrix \(T\) of dyadic size \(m=2^d\). +Let \(Z\) be the nilpotent lower shift and \(E=T-ZTZ^T\), with rank \(r\). +The target is the dirty shear +\[ +S(T)(x,y)=(x,y+Tx) +\] +on exactly \(2m\) existing complex wires. I require no zero bank. +The input lies on an \(F_0\)-bit dyadic grid and has norm at most \(R\ge1\); +the requested output error is \(2^{-t}\), with integer \(t\ge1\). +I use a rational bound \(B_E\ge\|E\|_2\). The executable takes the maximum +of the row and column component sums, which bounds both induced norms. + +I execute the frozen PR4198 rational producer. For each nonzero orthogonal +residual \(v_j\) of \(E\), it chooses a dyadic \(c_j\) with +\(\|v_j\|\le c_j<2\|v_j\|\), then emits +\[ +g_j=v_j/c_j,\qquad +a_j=\frac{c_j}{v_j^*v_j}v_j^*E. +\] +The orthogonal projectors sum to the identity on the column space of \(E\). +Thus +\[ +E=\sum_{j=1}^r g_ja_j,\quad +\|g_j\|\le1,\quad \|a_j\|<2B_E,\quad +T=\sum_{j=1}^r L(g_j)U(a_j). +\tag{1} +\] +The last identity follows by telescoping +\(T=\sum_{k=0}^{m-1}Z^kE(Z^T)^k\). The \(a_j\) are coefficient rows; +I apply no further conjugation. + +I retain the producer's scale +\[ +\Gamma_j=\min\{2^e:2^e\ge\max(1,2\ell(a_j))\},\quad +\ell(a)=\sum_k(|\Re a_k|+|\Im a_k|), +\] +and put \(P_j=I+U(a_j)/\Gamma_j\), \(G_j=\Gamma_jL(g_j)\). +Then \(P_j\) satisfies the component-l1 admission used in +[TRIANGLE-PROOF.md](TRIANGLE-PROOF.md), and +\(\|P_j\|\le3/2,\ \|P_j^{-1}\|\le2\). +The source realizes the upper triangular map by reversing its source wires. + +I use the rational bounds +\[ +\alpha=2rmB_E,\qquad K_G=2m+16mB_E. +\tag{2} +\] +Indeed \(\|L(g_j)U(a_j)\|\le2mB_E\), and +\[ +\Gamma_j\sqrt m + <2\sqrt m+8\sqrt2\,mB_E\le K_G. +\] +The scaled shear's coefficient component sum is at most +\(\sqrt2K_G<2K_G\). I budget the approximation of \(G_j\) itself; +approximating \(L(g_j)\) first and multiplying by an uncharged +\(\Gamma_j\) would not meet this contract. + +## 2. Paired compensation and additive defects + +I retain the source's four calls: +\[ +P_j;\quad S(G_j);\quad P_j^{-1};\quad S(-G_j). +\tag{3} +\] +I compile \(\widehat P_j\) and reuse its exact stored inverse. +For the positive shear I use the certified FFT construction. +For its compensation I reuse the same unitary tables and negate the +stored spectral multipliers, preserving the source's cyclic-then-negacyclic +order. The two unrounded shear endpoints are exactly +\(S(\widehat G_j)\) and \(S(-\widehat G_j)\). + +Negating a nearest-grid value gives a nearest-grid value of the negative +exact quantity, so this paired rule retains the same multiplier-error +bound. At half-grid ties it need not agree with independent rounding +toward positive infinity. The verifier includes that discrepancy as a +control. I do not replace the negative child with a reversed whole +positive-shear program. + +For \(\|\widehat P_j-P_j\|,\|\widehat G_j-G_j\|\le\eta\le1/4\), +the exact stored-coefficient program has endpoint +\[ +(x,y)\longmapsto(x,y+\widehat G_j(\widehat P_j-I)x). +\] +Write \(T_j=G_j(P_j-I)\) and +\(\widehat T_j=\widehat G_j(\widehat P_j-I)\). Then +\[ +\|\widehat T_j-T_j\| +\le K_G\eta+(1/2+\eta)\eta\le(K_G+1)\eta. +\] +The completed terms restore the same source bank, so their endpoints +compose to \(S(\sum_j\widehat T_j)\). Consequently +\[ +\boxed{\|\sum_j\widehat T_j-T\|\le +\Delta:=r(K_G+1)\eta.} +\tag{4} +\] +I need no generic product-of-perturbed-matrices bound across completed +generator terms. This additive estimate uses exact stored-coefficient +restoration; data-rounding errors receive a separate bound below. + +## 3. Checkpoint prefixes and suffixes + +Assume \(\Delta\le1\). Any preceding or remaining sum of approximate terms +has norm at most \(a_*=\alpha+1\). +For the current term write \(p=\widehat P_j\), \(h=\widehat G_j\), +\(A=\sum_{i2/3\), while the endpoint error is below \(1/6\). +For the true rotated phase \(v_0\), put +\[ +\tau=\Im(v_0)/(1+\Re(v_0)). +\] +The denominator is at least \(3/2\). Compute quotient extrema at the box +corners, round their midpoint to the nearest \((P+4)\)-grid value \(t_0\), +and form +\[ +\sigma=C\frac{1-t_0^2+2it_0}{1+t_0^2}. +\] +Then \(|\sigma|=1\) exactly, \(|t_0|<1\), and +\[ +|\sigma-e^{2\pi ik/n}|\le4E_B+2^{-P-4}. +\] +This follows from the quotient bound \(2E_B\) and the stereographic +derivative bound \(2/(1+t^2)\le2\). + +The executable checks the **maximum endpoint** component-\(\ell^1\) +distance from \(\sigma\) to the entire oracle box. That distance is at +most \(10E_B+2^{-P-3}<\nu\). A minimum distance to the box would not +certify the approximation. + +The parameter denominator is at most \(2^{P+4}\). Unit-phase components +therefore have numerator/denominator height at most \(2P+9\), with +slightly larger allowances in the code. The oracle's offline exact +Machin sums, factorials and phase-description arithmetic use polynomial +time and space in \(P\) and the input description length. A long external +phase numerator is an input cost, not a free oracle argument. + +## 3. Exactly unitary butterflies and FFT approximation + +Use the certified phase \(a+ib\) approximating \(e^{i\pi/4}\), and a +certified twiddle \(w\). Define +\[ +B(a,b,w)= +\begin{pmatrix}a&bw\\b&-aw\end{pmatrix} += +\begin{pmatrix}a&b\\b&-a\end{pmatrix}\operatorname{diag}(1,w). +\] +Since \(a^2+b^2=|w|^2=1\), this matrix is exactly unitary. +Its adjoint is +\(\operatorname{diag}(1,\bar w)\left[\begin{smallmatrix}a&b\\b&-a\end{smallmatrix}\right]\). +It is generally not the same expression with \(w\) merely conjugated. + +The splitter difference from the ideal \(1/\sqrt2\) splitter has norm +\(|(a+ib)-e^{i\pi/4}|\). Hence each butterfly differs from its ideal gate +by at most \(2\nu\). Disjoint gates give the same layer bound. I retain +the source's exact bit reversal and its \(d\) butterfly layers, so +\[ +\boxed{\|\widetilde F-F\|_2\le2d\nu,\qquad +F=\mathrm{DFT}_{-}/\sqrt n.} +\] +We telescope products of unitary layers, so no amplification factor +appears. Every exact FFT prefix and suffix is unitary. Reverse the layers +and adjoint the stored gates to obtain \(\widetilde F^*\); do not generate +a new approximate inverse phase table. + +## 4. Spectral preparation and Toeplitz approximation + +For a kernel \(k\), the ideal circulant eigenvalues are the +**unnormalized** DFT +\[ +\lambda_\ell=\sum_j k_j e^{-2\pi i\ell j/n}. +\] +I approximate each root independently by a certified rational unit, +sum exactly offline, and round each multiplier component once to the +nearest \(S\)-grid value. If \(\|k\|_1\le M_k\), then +\[ +\|\widehat\lambda-\lambda\|_\infty\le M_k\nu+v,\qquad +\|\widehat\lambda\|_\infty\le M_k+v. +\] +The second bound uses exact unit modulus before rounding. It does not +add a spurious phase-error magnitude factor. + +I store the rounded multipliers before online use. The exact offline +sum can have a denominator whose size grows with \(n\); carrying that +sum into data gates would violate the stated coefficient-width bound. + +The actual circulant child applies +\[ +\widetilde C_k=\widetilde F^*\operatorname{diag}(\widehat\lambda)\widetilde F, +\] +so +\[ +\|\widetilde C_k-C_k\|_2 +\le M_k\nu+v+2M_k\|\widetilde F-F\|_2. +\] +Exact unitarity alone does not make \(\widetilde F\) an exact diagonalizer. + +For \(D=\operatorname{diag}(e^{\pi ij/n})\), the source decomposition is +\[ +T=C_g+D^*C_{Dh}D. +\] +I generate each diagonal entry of \(\widetilde D\) independently. Thus +\(\|\widetilde D-D\|_2\le\nu\), without a powers' accumulated error or +an \(nP\)-bit phase representation. The actual schedule is the first +circulant shear, the common twist \(\widetilde D\) on both banks, the +circulant shear prepared from \(\widetilde D h\), and the inverse twist. + +Twisting preserves \(\|h\|_1\). Charge \(M_h\nu\) for changing the +second kernel and \(2M_h\nu\) for its two exterior twists. Combining +these terms with the two circulant estimates gives +\[ +\|\widetilde T-T\|_2 +\le(4d+1)M_g\nu+(4d+4)M_h\nu+2v +\le\boxed{\delta}. +\] +In exact rational arithmetic the complete map is +\((x,y)\mapsto(x,y+\widetilde T x)\). The input matrix may be zero +or singular because the pointwise shear gates remain invertible. + +## 5. Prefix bounds and rounded execution + +Every non-spectral gate stays within one bank and is unitary. +Each of the two spectral layers has orientation +\[ +S(A)=\begin{pmatrix}I&0\\A&I\end{pmatrix}. +\] +For a bank unitary \(U=\operatorname{diag}(U_x,U_y)\), +\[ +US(A)U^*=S(U_yAU_x^*). +\] +Move the bank unitaries to the ends of any contiguous gate segment. +The remaining shear matrix is a sum of at most two unitarily transformed +spectral matrices. Their norms are at most \(M_g+v\) and \(M_h+v\). +Consequently every segment, inverse segment, prefix and suffix has norm +at most +\[ +\boxed{H=1+M_g+M_h+2v.} +\] +Masked partial spectral layers satisfy the same estimate. +I do not multiply a shear bound once per stage. + +Nearest component rounding has error at most \(u/2\); a complex output +has error at most \(u/\sqrt2\). A disjoint layer updates at most \(2n\) +complex wires, so its error vector has norm at most \(\sqrt n\,u\). +The source schedule has at most \(L=4d+8\) layers, including empty +bit-reversal layers. Transporting each local residual through its +remaining exact suffix gives +\[ +\|\widehat z-\widetilde S z_0\|_2 +\le LH\sqrt n\,u\le LHnu. +\] +Every pre-round and rounded state has amplitude at most +\[ +A=H(R+Lnu). +\] +The source-bank error receives no coefficient-approximation term, +because the unrounded rational circuit restores that bank exactly. +Fixed-grid execution can still have nonzero source-restoration error; +the verifier includes a small-grid witness. + +The precision choices in Section 1 close the requested budget: +\[ +4(d+1)M\nu R<2^{-t}/4,\quad 2vR<2^{-t}/4,\quad +LHnu<2^{-t}/2. +\] +Here \(v\le1/32\) and \(H 0 and n & (n - 1) == 0 and matrix.shape == (n, n), "dyadic full Toeplitz shape") + need(all(matrix[i][j] == (matrix[i - j][0] if i >= j else matrix[0][j - i]) + for i in range(n) for j in range(n)), "exact full Toeplitz admission") + leading = [] + for b in range(1, n + 1): + inverse = solve(matrix[:b, :b], Matrix(eye(b))) + leading.append((b, norm_upper(inverse))) + return max(Q(1), norm_upper(matrix)), max(Q(1), *(v for _, v in leading)), leading + + +def allocation(n, lam, inverse_bound, radius, target, incoming): + need(type(n) is int and n > 0 and n & (n - 1) == 0, "full dyadic dimension") + need(all(type(x) in (int, Q) and x >= 1 for x in (lam, inverse_bound, radius)), + "exact full conditioning and radius bounds") + need(type(target) is int and target >= 1 and type(incoming) is int and incoming >= 0, + "full target and incoming precision") + lam, inverse_bound, radius = (Q(x) for x in (lam, inverse_bound, radius)) + d, chi = n.bit_length() - 1, lam * inverse_bound + descendant = lam + lam ** 2 * inverse_bound + b = 256 * n ** 3 * chi ** 12 + calls, rank = 3 * n - 2, 2 * d + 8 + alpha, kg = 8 * rank * n ** 2 * chi, 2 * n + 64 * n ** 2 * chi + generator = max(Q(3200 * n ** 3), 2 + 2 * kg) * (10 + 4 * alpha + 2 * kg) + fragment = max(generator, descendant + 1, 2 * inverse_bound) + execution = 2 * max(Q(1), lam, inverse_bound) * b * fragment + layers = 1 + 2 * d * rank * (8 * d * d + 32 * d + 18) + gates = n * (1 + d * rank * (4 * d * d + 24 * d + 22)) + s = target + ceil_log2(4 * calls * b ** 3 * (radius + 1)) + eta = Q(1, 1 << s) + theta = calls * b ** 2 * eta + data = max(incoming, target + ceil_log2(2 * layers * execution * n)) + coefficient = 2 * calls * b ** 3 * eta + rounding = layers * execution * n * Q(1, 1 << data) + need(theta <= Q(1, 4) and eta <= 1 / (2 * inverse_bound) + and coefficient * radius + rounding <= Q(1, 1 << target), "complete full precision budget") + amplitude = execution * radius + layers * execution ** 2 * n * Q(1, 1 << data) + integer = max(1, ceil_log2(max(Q(1), amplitude))) + 2 + scalar_height = 2 * s + 2 * ceil_log2(descendant + 2) + 8 + return dict(n=n, lam=lam, inverse_bound=inverse_bound, radius=radius, + target_bits=target, incoming_bits=incoming, conditioning=chi, + descendant_bound=descendant, completed_prefix_bound=b, + atomic_calls=calls, rank_bound=rank, displacement_bound=4 * n * chi, + generator_fragment_bound=generator, fragment_bound=fragment, + execution_bound=execution, segment_bound=execution ** 2, + layer_bound=layers, gate_bound=gates, child_target_bits=s, eta=eta, + theta=theta, data_bits=data, coefficient_error=coefficient, + rounding_error=rounding, total_error=coefficient * radius + rounding, + amplitude_bound=amplitude, integer_bits=integer, + scalar_height_bound=scalar_height, coefficient_height_bound=scalar_height) + + +class Plan(list): + def __init__(self, layers, atoms): + super().__init__(layers) + self.atoms = atoms + + def __add__(self, other): + need(isinstance(other, Plan), "typed full-plan composition") + return Plan(list(self) + list(other), [dict(a) for a in self.atoms] + + [dict(a, start=a["start"] + len(self), stop=a["stop"] + len(self)) + for a in other.atoms]) + + +class Full: + def __init__(self, bound, raw): + self.bound, self.raw = dict(bound), raw + self.producer = load_producer(raw) + self.phase_pool, self.nodes = {}, {} + selected = {"full_plan", "parallel", "inverse"} + tree = ast.parse(raw[SOURCE + "check_circuits.py"]) + tree.body = [copy.deepcopy(node) for node in tree.body + if isinstance(node, ast.FunctionDef) and node.name in selected] + need({node.name for node in tree.body} == selected, "complete full source function set") + self.source_ast = sha256(ast.dump(tree).encode()).hexdigest() + self.stats = {"full_nodes": 0} + self.namespace = dict(STATS=self.stats, + np=SimpleNamespace(linalg=SimpleNamespace(solve=solve, inv=inverse_gate)), + lowdisp_shear=self.shear, diag=self.leaf) + exec(compile(tree, "", "exec", optimize=0), + self.namespace) + self.source_full = self.namespace["full_plan"] + self.source_parallel = self.namespace["parallel"] + self.source_inverse = self.namespace["inverse"] + self.namespace["parallel"] = self.parallel + self.namespace["full_plan"] = self.node + + def interval(self, indices): + positions = [self.positions[i] for i in indices] + a, b = min(positions), max(positions) + 1 + need(positions == list(range(a, b)), "ordered contiguous logical interval") + return a, b + + def node(self, indices, matrix): + interval = self.interval(indices) + need(interval not in self.nodes, "one exact Schur reference per node") + self.nodes[interval] = dict(indices=list(indices), matrix=Matrix(matrix)) + need(norm_upper(matrix) <= self.bound["descendant_bound"], "component descendant bound") + result = self.source_full(indices, matrix) + self.nodes[interval]["layers"] = len(result) + return result + + def leaf(self, indices, values): + need(len(indices) == len(values) == 1, "scalar full leaf") + intended = values[0] + need(intended.norm2() >= 1 / self.bound["inverse_bound"] ** 2 + and intended.norm2() <= self.bound["descendant_bound"] ** 2, "full pivot margins") + stored = rounded(intended, self.bound["child_target_bits"]) + need((stored - intended).norm2() <= self.bound["eta"] ** 2 + and stored.norm2() >= 1 / (4 * self.bound["inverse_bound"] ** 2), "stored full pivot margin") + gate = Mat([[stored]]) + need(stored * inverse_gate(gate)[0][0] == 1, "same stored full scalar inverse") + a, b = self.interval(indices) + atom = dict(id=f"D:{a}", kind="leaf", interval=(a, b), indices=list(indices), + start=0, stop=1, ideal=Mat([[intended]]), + local_plan=[[((0,), gate)]], error=self.bound["eta"], + coefficient_height=self.bound["scalar_height_bound"]) + return Plan([[(tuple(indices), gate)]], [atom]) + + def shear(self, target, source, matrix): + m = len(source) + positions = sorted(self.positions[i] for i in source + target) + a, b = min(positions), max(positions) + 1 + need(positions == list(range(a, b)) and b - a == 2 * m, "owned full shear interval") + upper = self.positions[target[0]] < self.positions[source[0]] + need([self.positions[i] for i in target] == list(range(a if upper else a + m, + a + m if upper else b)) + and [self.positions[i] for i in source] == list(range(a + m if upper else a, + b if upper else a + m)), + "source full shear orientation") + need(norm_upper(matrix) <= self.bound["conditioning"], "full K/L conditioning") + items, displacement, be, heights = prepare(matrix, self.producer) + need(len(items) <= self.bound["rank_bound"] and be <= self.bound["displacement_bound"], + "retained displacement rank and envelope") + child, triangle = generator_allocation(m, len(items), be, 1, + self.bound["child_target_bits"], 0) + need(child["coefficient_error"] <= self.bound["eta"] + and child["execution_bound"] <= self.bound["generator_fragment_bound"], + "full atomic generator contract") + p = child["phase_bits"] + if p not in self.phase_pool: + self.phase_pool[p] = Phases(p, self.raw[ORACLE]) + program = Program(child, triangle, self.raw, self.phase_pool[p]) + local = program.plan(items) + mapping = source + target + plan = [[(tuple(mapping[i] for i in ids), gate) for ids, gate in layer] for layer in local] + kind = "upper" if upper else "lower" + atom = dict(id=f"{'R' if upper else 'L'}:{a}:{b}", kind=kind, + interval=(a, b), indices=list(mapping), start=0, stop=len(plan), + ideal=Mat(matrix), local_plan=local, error=child["coefficient_error"], + coefficient_height=child["coefficient_height_bound"], + program=program, generator_items=items, displacement=displacement, + producer_heights=heights) + return Plan(plan, [atom]) + + def parallel(self, left, right): + left_ids = {i for a in left.atoms for i in a["indices"]} + right_ids = {i for a in right.atoms for i in a["indices"]} + need(not left_ids.intersection(right_ids), "disjoint full sibling supports") + return Plan(self.source_parallel(left, right), + [dict(a) for a in left.atoms] + [dict(a) for a in right.atoms]) + + def plan(self, matrix, indices=None): + n = len(matrix) + need(n == self.bound["n"], "allocated full dimension") + indices = list(range(n)) if indices is None else list(indices) + need(all(type(i) is int for i in indices) and sorted(indices) == list(range(n)), + "full root wire permutation") + self.positions = {wire: pos for pos, wire in enumerate(indices)} + result = self.node(indices, Matrix(matrix)) + need(self.stats["full_nodes"] == 2 * n - 1 and len(result.atoms) == 3 * n - 2 + and len({a["id"] for a in result.atoms}) == len(result.atoms), "complete full recursion") + need(len(result) <= self.bound["layer_bound"] + and sum(map(len, result)) <= self.bound["gate_bound"], "full source schedule bounds") + self.bound["coefficient_height_bound"] = max(a["coefficient_height"] for a in result.atoms) + self.bound["temporary_bits"] = (2 * self.bound["data_bits"] + + 8 * self.bound["coefficient_height_bound"] + 2 * self.bound["integer_bits"] + 32) + return result + + def inverse(self, plan): + return self.source_inverse(plan) diff --git a/research/rational-fft-fixed-h/vendor/full_source/generator.py b/research/rational-fft-fixed-h/vendor/full_source/generator.py new file mode 100644 index 0000000000..d761d3c88b --- /dev/null +++ b/research/rational-fft-fixed-h/vendor/full_source/generator.py @@ -0,0 +1,233 @@ +"""Exact generators connected to the retained compensated gate program. + +Jose Manuel (jmbarrancoml), with substantial OpenAI Codex assistance. +Apache-2.0. Macro/compiler: Vlado/PR1673; certified phases: Chafik/PR628. +""" +import ast +import copy +from fractions import Fraction as Q +from hashlib import sha256 +import json + +from exact import C, Mat, Vec, ceil_log2, component_l1, inverse_gate, need, norm2 +from engine import Builder, Phases, split_kernel +from triangular import Triangle, allocation as triangle_allocation + +SOURCE = "research/capped-phase-fusion/evidence/inplace_algebra_ultra/" +PRODUCER = "research/rational-unitary-fft/vendor/stein_rational.py" +PROVENANCE = "research/rational-unitary-fft/vendor/STEIN-PROVENANCE.json" +ORACLE = "research/cubic-cyclic-inverse/vendor/interval_kernel.py" + + +def load_producer(raw): + record = json.loads(raw[PROVENANCE]) + need(record["source_commit"] == "e457d5e912a27307d725827aeca7198169320b84" + and record["source_git_blob"] == "71b240a460a185ee2fe926cd189b56aeac6f997f" + and sha256(raw[PRODUCER]).hexdigest() == record["sha256"], + "original producer identity and retained bytes") + namespace = {"__name__": "_frozen_rational_stein"} + exec(compile(raw[PRODUCER], "", "exec", optimize=0), namespace) + return namespace + + +def lower(vector): + return Mat([[vector[i - j] if i >= j else C() for j in range(len(vector))] + for i in range(len(vector))]) + + +def prepare(matrix, producer): + matrix = Mat(matrix) + m = len(matrix) + need(m > 0 and m & (m - 1) == 0 and matrix.shape == (m, m), "dyadic generator matrix") + displacement = Mat([[matrix[i][j] - (matrix[i - 1][j - 1] if i and j else C()) + for j in range(m)] for i in range(m)]) + old_c = producer["C"] + old_e = producer["Matrix"]([[old_c(x.re, x.im) for x in row] for row in displacement]) + old_items, heights = producer["generators"](old_e) + items = [dict(column=item["column"], gamma=item["gamma"], + scale=item["scale"], scale_exponent=item["scale_exponent"], + g=Vec(C(x.re, x.im) for x in item["g"]), + row=Vec(C(x.re, x.im) for x in item["row"])) for item in old_items] + bound_e = max(max(component_l1(row) for row in displacement), + max(component_l1(column) for column in zip(*displacement))) + need(all(norm2(item["g"]) <= 1 and norm2(item["row"]) < 4 * bound_e ** 2 + and 2 * component_l1(item["row"]) <= item["gamma"] for item in items), + "actual bounded rational generator output") + return items, displacement, bound_e, heights + + +def allocation(m, rank, bound_e, radius, target, incoming): + need(type(m) is int and m > 0 and m & (m - 1) == 0, "dyadic macro size") + need(type(rank) is int and 0 <= rank <= m and type(bound_e) in (int, Q) + and bound_e >= 0, "generator count and exact displacement bound") + need(type(radius) in (int, Q) and radius >= 1 and type(target) is int + and target >= 1 and type(incoming) is int and incoming >= 0, "global generator input") + d = m.bit_length() - 1 + kg, alpha = 2 * m + 16 * m * bound_e, 2 * rank * m * bound_e + checkpoint = 10 + 4 * alpha + 2 * kg + child = max(Q(3200 * m ** 3), 2 + 2 * kg) + execution = child * checkpoint + depth = rank * (8 * d * d + 32 * d + 18) + gates = rank * m * (4 * d * d + 24 * d + 22) + s = target + ceil_log2(4 * (rank + 1) * (kg + 1) * (radius + 1)) + eta = Q(1, 1 << s) + triangle = triangle_allocation(m, 1, s) + envelope = max(Q(6), 2 * kg) + p = max(triangle["phase_bits"], s + ceil_log2(32 * (d + 1) * (envelope + 1))) + mult = max(triangle["multiplier_bits"], s + 4) + triangle = dict(triangle, phase_bits=p, multiplier_bits=mult, + coefficient_height_bound=max(4 * p + 24, mult + 5, + 2 * triangle["leaf_bits"] + 16)) + data = max(incoming, target + ceil_log2(max(Q(1), 2 * depth * execution * m))) + u = Q(1, 1 << data) + defect = rank * (kg + 1) * eta + rounding = depth * execution * m * u + need(defect <= 1 and defect * radius + rounding <= Q(1, 1 << target), + "complete compensated-generator target budget") + amplitude = execution * radius + depth * execution ** 2 * m * u + integer = max(1, ceil_log2(max(Q(1), amplitude))) + 2 + height = max(4 * p + 24, mult + ceil_log2(envelope + 2) + 2, + 2 * triangle["leaf_bits"] + 16) + return dict(m=m, rank=rank, displacement_bound=bound_e, radius=radius, + target_bits=target, incoming_bits=incoming, kg=kg, alpha=alpha, + checkpoint_bound=checkpoint, child_bound=child, execution_bound=execution, + segment_bound=execution ** 2, layer_bound=depth, gate_bound=gates, + child_target_bits=s, eta=eta, phase_bits=p, multiplier_bits=mult, + leaf_bits=triangle["leaf_bits"], data_bits=data, + multiplier_envelope=envelope, coefficient_error=defect, + rounding_error=rounding, total_error=defect * radius + rounding, + amplitude_bound=amplitude, integer_bits=integer, + coefficient_height_bound=height, + temporary_bits=2 * data + 8 * height + 2 * integer + 32), triangle + + +def macro_template(source): + tree = ast.parse(source) + function = next(n for n in tree.body if isinstance(n, ast.FunctionDef) + and n.name == "lowdisp_shear") + loop = next(n for n in function.body if isinstance(n, ast.For)) + start = next(i for i, n in enumerate(loop.body) if isinstance(n, ast.Assign) + and isinstance(n.targets[0], ast.Name) and n.targets[0].id == "p") + tail = copy.deepcopy(loop.body[start:]) + need([type(n).__name__ for n in tail] == + ["Assign", "AugAssign", "Assign", "AugAssign", "AugAssign", "AugAssign", "AugAssign"], + "retained seven-statement macro") + identity = sha256(ast.dump(ast.Module(body=tail, type_ignores=[])).encode()).hexdigest() + need(identity == "2a2e0f16a779597fa555857307451831bace11e66861fdc3caf26d2373a21210", + "unchanged source macro identity") + function = ast.FunctionDef( + name="emit", args=ast.arguments(posonlyargs=[], args=[ast.arg(arg=n) for n in + ("target", "source", "g", "h", "gamma")], vararg=None, kwonlyargs=[], + kw_defaults=[], kwarg=None, defaults=[]), + body=[ast.Assign(targets=[ast.Name(id="result", ctx=ast.Store())], + value=ast.List(elts=[], ctx=ast.Load()))] + tail + + [ast.Return(value=ast.Name(id="result", ctx=ast.Load()))], + decorator_list=[], type_params=[]) + module = ast.fix_missing_locations(ast.Module(body=[function], type_ignores=[])) + return compile(module, "", "exec", optimize=0), identity + + +class PairedNegative(Builder): + def __init__(self, positive, source): + super().__init__(positive.phases, positive.multiplier_bits, source) + self.positive = positive + + def spectrum(self, kernel): + row = self.positive.spectra[len(self.spectra)] + need(Vec(kernel) == Vec(-x for x in row["kernel"]), "paired negative spectral kernel") + result = Vec(-x for x in row["multipliers"]) + bound = self.bounds[len(self.spectra)] + need(bound == row["bound"], "paired spectral envelope") + v = Q(1, 1 << self.multiplier_bits) + need(all(x.norm2() <= (bound + v) ** 2 for x in result), "negative stored multiplier bound") + self.spectra.append(dict(kernel=Vec(kernel), multipliers=result, bound=bound)) + return result + + +class Program: + def __init__(self, bound, triangle_bound, raw, phase_oracle=None): + self.bound, self.triangle_bound, self.raw = bound, triangle_bound, raw + self.phases = phase_oracle if phase_oracle is not None else Phases(bound["phase_bits"], raw[ORACLE]) + need(self.phases.bits == bound["phase_bits"], "common generator phase precision") + self.source = list(range(bound["m"])) + self.target = list(range(bound["m"], 2 * bound["m"])) + self.stages, self.position, self.calls = [], 0, [] + self.pending = None + self.last_triangle = None + self.triangles, self.shears = [], [] + code, self.macro_ast = macro_template(raw[SOURCE + "check_circuits.py"]) + namespace = dict(triangle_plan=self.triangle, toeplitz_shear=self.shear, + inverse=self.inverse, lower=lower) + exec(code, namespace) + self.emit = namespace["emit"] + + def record(self, kind, plan, **details): + self.calls.append(kind) + self.stages.append(dict(kind=kind, start=self.position, stop=self.position + len(plan), + **details)) + self.position += len(plan) + return plan + + def triangle(self, indices, coefficients): + builder = Triangle(self.triangle_bound, self.raw[SOURCE + "check_circuits.py"], + self.raw[ORACLE], self.phases) + local_indices = [self.source.index(i) for i in indices] + local = builder.plan(coefficients, local_indices) + plan = [[(tuple(self.source[i] for i in ids), gate) for ids, gate in layer] for layer in local] + need(builder.bound["coefficient_error"] <= self.bound["eta"], "triangle operator allowance") + self.last_triangle = plan + self.triangles.append(builder) + return self.record("triangle", plan, symbol=Vec(coefficients), indices=list(indices), + operator_error=builder.bound["coefficient_error"]) + + def inverse(self, plan): + need(plan is self.last_triangle, "inverse of the same stored triangle") + result = [[(ids, inverse_gate(gate)) for ids, gate in layer] for layer in reversed(plan)] + return self.record("inverse", result) + + def shear(self, target, source, matrix): + need(target == self.target and source == self.source, "literal complete dirty banks") + if self.pending is None: + builder = Builder(self.phases, self.bound["multiplier_bits"], + self.raw[SOURCE + "check_circuits.py"]) + local = builder.plan(matrix) + self.pending = (builder, matrix, local) + kind = "positive" + else: + positive, original, positive_plan = self.pending + need(matrix == -1 * original, "literal negative compensation matrix") + builder = PairedNegative(positive, self.raw[SOURCE + "check_circuits.py"]) + local = builder.plan(matrix) + m = len(source) + need(len(local) == len(positive_plan), "same compensation layer schedule") + for a, b in zip(positive_plan, local): + need(len(a) == len(b), "same compensation gate schedule") + for (ids, g), (other_ids, h) in zip(a, b): + crossing = len(ids) == 2 and ids[0] < m <= ids[1] + expected = Mat([[1, 0], [-g[1][0], 1]]) if crossing else g + need(ids == other_ids and h == expected, "paired coefficients preserve source order") + self.pending = None + kind = "compensation" + g, h = split_kernel(matrix) + envelope = component_l1(g) + component_l1(h) + defect = (4 * len(source).bit_length() * envelope * self.phases.nu + + Q(2, 1 << self.bound["multiplier_bits"])) + need(envelope <= self.bound["multiplier_envelope"] and defect <= self.bound["eta"], + "scaled shear absolute operator allowance") + self.shears.append(builder) + mapping = source + target + plan = [[(tuple(mapping[i] for i in ids), gate) for ids, gate in layer] for layer in local] + return self.record(kind, plan, matrix=matrix, operator_error=defect, envelope=envelope) + + def plan(self, items): + result = [] + for item in items: + before = Vec(item["row"]) + result += self.emit(self.target, self.source, item["g"], item["row"], item["gamma"]) + need(before == item["row"], "fresh division preserves generator row") + need(self.pending is None and self.calls == + ["triangle", "positive", "inverse", "compensation"] * len(items), + "complete retained compensated call sequence") + need(self.position == len(result) <= self.bound["layer_bound"] + and sum(map(len, result)) <= self.bound["gate_bound"], "full program gate/layer budget") + return result diff --git a/research/rational-fft-fixed-h/vendor/full_source/triangular.py b/research/rational-fft-fixed-h/vendor/full_source/triangular.py new file mode 100644 index 0000000000..2bd8de01a7 --- /dev/null +++ b/research/rational-fft-fixed-h/vendor/full_source/triangular.py @@ -0,0 +1,200 @@ +"""Source-bound rational parity recursion with a common precision allocation. + +Jose Manuel (jmbarrancoml), with substantial OpenAI Codex assistance. +Apache-2.0. Parity recursion and its stability premise: Vlado/PR1673. +""" +import ast +import copy +from fractions import Fraction as Q +from hashlib import sha256 +from types import SimpleNamespace + +from exact import (C, Mat, Vec, ceil_log2, component_l1, inverse_gate, need, + rounded) +from engine import Builder, Phases, split_kernel + + +def allocation(n, radius, target, incoming=0): + need(type(n) is int and n > 0 and n & (n - 1) == 0, "dyadic triangle size") + need(type(radius) in (int, Q) and radius >= 1, "exact triangle radius") + need(type(target) is int and target >= 1, "positive triangle target") + need(type(incoming) is int and incoming >= 0, "incoming triangle grid") + d = n.bit_length() - 1 + q, depth, bound = 2 * d + 1, 4 * d * d + 12 * d + 1, 100 * n ** 3 + s = target + ceil_log2(4 * q * bound ** 3 * (radius + 1)) + p = s + ceil_log2(48 * (d + 1)) + multiplier = s + 2 + data = max(incoming, target + ceil_log2(64 * depth * bound * n)) + epsilon, u = Q(1, 1 << s), Q(1, 1 << data) + theta = q * bound ** 2 * epsilon + coefficient = 2 * q * bound ** 3 * epsilon + rounding = 32 * depth * bound * n * u + need(theta <= Q(1, 4) and coefficient * radius + rounding <= Q(1, 1 << target), + "complete triangular error allocation") + amplitude = 32 * bound * radius + 1024 * depth * bound ** 2 * n * u + integer = max(1, ceil_log2(max(Q(1), amplitude))) + 2 + height = max(4 * p + 24, multiplier + 5, 2 * s + 16) + return dict(n=n, radius=radius, target_bits=target, incoming_bits=incoming, + ideal_bound=bound, macro_layers=q, layer_bound=depth, + gate_bound=n * (2 * d * d + 8 * d + 1), + leaf_bits=s, phase_bits=p, multiplier_bits=multiplier, + data_bits=data, epsilon=epsilon, theta=theta, + coefficient_error=coefficient, rounding_error=rounding, + total_error=coefficient * radius + rounding, + physical_prefix_suffix_bound=32 * bound, + segment_bound=1024 * bound ** 2, amplitude_bound=amplitude, + integer_bits=integer, coefficient_height_bound=height, + temporary_bits=2 * data + 8 * height + 2 * integer + 32) + + +class Plan(list): + """Actual layers plus consecutive, disjoint-support macro partitions.""" + def __init__(self, layers, blocks): + super().__init__(layers) + self.blocks = blocks + + def __add__(self, other): + need(isinstance(other, Plan), "typed triangular composition") + offset = len(self) + blocks = [dict(block) for block in self.blocks] + blocks += [dict(start=b["start"] + offset, stop=b["stop"] + offset, + pieces=list(b["pieces"])) for b in other.blocks] + return Plan(list(self) + list(other), blocks) + + +class Concatenate: + def __getitem__(self, parts): + values = [] + for part in parts: + values.extend(part if isinstance(part, (list, tuple)) else [part]) + return Vec(values) + + +def convolution(a, b): + return Vec(sum((a[j] * b[k - j] for j in range(len(a)) + if 0 <= k - j < len(b)), C()) + for k in range(len(a) + len(b) - 1)) + + +class AdaptTriangle(ast.NodeTransformer): + def __init__(self): + self.changes = [] + + def visit_Assign(self, node): + if (len(node.targets) == 1 and isinstance(node.targets[0], ast.Subscript) + and isinstance(node.targets[0].slice, ast.Constant) + and node.targets[0].slice.value == "max_triangle_coefficient"): + self.changes.append("exact-component-diagnostic") + node.value = ast.parse( + "max(STATS['max_triangle_coefficient'], component_l1(f))", + mode="eval").body + return node + return self.generic_visit(node) + + def visit_Return(self, node): + if (isinstance(node.value, ast.Call) and isinstance(node.value.func, ast.Name) + and node.value.func.id == "diag"): + self.changes.append("stored-rounded-leaf") + node.value.func.id = "leaf_diag" + return self.generic_visit(node) + + +class Triangle: + def __init__(self, bound, circuit_source, oracle_source, phase_oracle=None): + self.bound, self.circuit_source = bound, circuit_source + self.phases = phase_oracle if phase_oracle is not None else Phases(bound["phase_bits"], oracle_source) + need(self.phases.bits == bound["phase_bits"], "common triangular phase precision") + base = Builder(self.phases, bound["multiplier_bits"], circuit_source) + self.source_parallel = base.namespace["parallel"] + self.source_inverse = base.namespace["inverse"] + self.children, self.leaves = [], [] + selected = {"toep", "lower", "invseries", "conv", "triangle_plan"} + tree = ast.parse(circuit_source) + tree.body = [copy.deepcopy(node) for node in tree.body + if isinstance(node, ast.FunctionDef) and node.name in selected] + need({node.name for node in tree.body} == selected, "complete parity source") + self.original_ast = sha256(ast.dump(tree).encode()).hexdigest() + adapter = AdaptTriangle() + tree = ast.fix_missing_locations(adapter.visit(tree)) + need(adapter.changes == ["exact-component-diagnostic", "stored-rounded-leaf"], + "exact declared parity adaptations") + self.adapted_ast = sha256(ast.dump(tree).encode()).hexdigest() + facade = SimpleNamespace(asarray=lambda a, dtype=None: Vec(a), + array=lambda a, dtype=None: Mat(a), + zeros=lambda n, dtype=None: Vec([0] * n), + convolve=convolution, r_=Concatenate()) + self.stats = {"triangle_nodes": 0, "max_triangle_coefficient": Q(0)} + self.namespace = dict(np=facade, STATS=self.stats, component_l1=component_l1, + leaf_diag=self.leaf, parallel=self.parallel, + toeplitz_shear=self.shear) + exec(compile(tree, "", "exec", optimize=0), + self.namespace) + + def leaf(self, indices, values): + need(len(indices) == len(values) == 1, "scalar triangular leaf") + wanted = values[0] + stored = rounded(wanted, self.bound["leaf_bits"]) + epsilon = self.bound["epsilon"] + need(Q(1, 4) <= wanted.norm2() <= 4, "ideal scalar leaf margin") + need((stored - wanted).norm2() <= epsilon ** 2 + and stored.norm2() >= Q(1, 16), "stored scalar leaf margin") + gate = Mat([[stored]]) + inverse = inverse_gate(gate) + need(stored * inverse[0][0] == 1, "same stored scalar reciprocal") + piece = dict(kind="leaf", ids=list(indices), ideal=Mat([[wanted]]), + error=epsilon, stored=stored, reciprocal=inverse[0][0]) + self.leaves.append(piece) + return Plan([[(tuple(indices), gate)]], + [dict(start=0, stop=1, pieces=[piece])]) + + def shear(self, target, source, matrix): + n = len(source) + need(len(set(source + target)) == 2 * n, "owned disjoint parity banks") + child = Builder(self.phases, self.bound["multiplier_bits"], self.circuit_source) + local = child.plan(matrix) + g, h = split_kernel(matrix) + envelope = component_l1(g) + component_l1(h) + v = Q(1, 1 << self.bound["multiplier_bits"]) + defect = 4 * (n.bit_length()) * envelope * self.phases.nu + 2 * v + need(envelope + 2 * v <= 6 and defect <= self.bound["epsilon"], + "bounded parity child and allocated operator defect") + mapping = source + target + layers = [[(tuple(mapping[i] for i in ids), gate) for ids, gate in layer] + for layer in local] + piece = dict(kind="shear", ids=list(mapping), ideal=matrix, error=defect, + envelope=envelope, source_ast=child.original_ast, + adapted_ast=child.adapted_ast, local_layers=len(local)) + self.children.append(piece) + return Plan(layers, [dict(start=0, stop=len(layers), pieces=[piece])]) + + def parallel(self, left, right): + need([(b["start"], b["stop"]) for b in left.blocks] + == [(b["start"], b["stop"]) for b in right.blocks], + "equal sibling macro boundary partitions") + blocks = [] + for a, b in zip(left.blocks, right.blocks): + pieces = a["pieces"] + b["pieces"] + ids = [i for piece in pieces for i in piece["ids"]] + need(len(ids) == len(set(ids)), "disjoint sibling macro supports") + blocks.append(dict(start=a["start"], stop=a["stop"], pieces=pieces)) + return Plan(self.source_parallel(left, right), blocks) + + def plan(self, symbol, indices=None): + symbol = Vec(symbol) + n = len(symbol) + need(n == self.bound["n"], "allocated triangle dimension") + need(component_l1(symbol - Vec([1] + [0] * (n - 1))) <= Q(1, 2), + "component-l1 near-identity admission") + indices = list(range(n)) if indices is None else list(indices) + need(len(indices) == n and sorted(indices) == list(range(n)), "triangle wire permutation") + result = self.namespace["triangle_plan"](indices, symbol) + need(self.stats["triangle_nodes"] == 2 * n - 1, "complete parity recursion tree") + need(len(result.blocks) == self.bound["macro_layers"] + and len(result) <= self.bound["layer_bound"] + and sum(map(len, result)) <= self.bound["gate_bound"], "complete schedule bounds") + return result + + def inverse(self, plan): + # Forward pieces certify forward defects. Inverse endpoint bounds use + # the inverse-product proof, not relabelled forward macro certificates. + return self.source_inverse(plan) diff --git a/research/rational-fft-fixed-h/vendor/full_source/vendor/STEIN-PROVENANCE.json b/research/rational-fft-fixed-h/vendor/full_source/vendor/STEIN-PROVENANCE.json new file mode 100644 index 0000000000..c602fb5bb3 --- /dev/null +++ b/research/rational-fft-fixed-h/vendor/full_source/vendor/STEIN-PROVENANCE.json @@ -0,0 +1,14 @@ +{ + "schema_version": 1, + "source_pr": "https://github.com/CrocSwap/integer-mult-bounds/pull/4198", + "source_commit": "e457d5e912a27307d725827aeca7198169320b84", + "source_path": "research/rational-stein-generators/rational.py", + "source_git_blob": "71b240a460a185ee2fe926cd189b56aeac6f997f", + "sha256": "3065fc5d6d22ec44de533b08fb0c639603449456383668e561dc46629e641502", + "retained_path": "research/rational-unitary-fft/vendor/stein_rational.py", + "byte_comparison": "I compared the retained bytes with git show at the original source commit before the local archive commit.", + "usage": "Execute the exact bounded rational generator producer as a coefficient-preparation dependency of the complete gate program.", + "proof": "https://github.com/CrocSwap/integer-mult-bounds/blob/e457d5e912a27307d725827aeca7198169320b84/research/rational-stein-generators/PROOF.md", + "credit": "Jose Manuel (jmbarrancoml), with substantial OpenAI Codex assistance. The producer refines Vlado Vojdanovski/PR1673's rational Gram-Schmidt route; Cinna Davis supplied the later caller. Apache-2.0.", + "history_boundary": "The original sibling PR commit is not an ancestor of PR4205. SOURCES.json will pin these local retained bytes at their archive commit, separately from this original-source identity. The later integration does not claim original sibling-branch ancestry." +} diff --git a/research/rational-fft-fixed-h/vendor/full_source/vendor/stein_rational.py b/research/rational-fft-fixed-h/vendor/full_source/vendor/stein_rational.py new file mode 100644 index 0000000000..71b240a460 --- /dev/null +++ b/research/rational-fft-fixed-h/vendor/full_source/vendor/stein_rational.py @@ -0,0 +1,196 @@ +"""Exact dyadically scaled Stein generators over Q(i). + +Jose Manuel (jmbarrancoml), with substantial OpenAI Codex assistance. +Apache-2.0. The construction refines the retained rational Gram-Schmidt route. +""" +from fractions import Fraction as Q + +MAX_HEIGHT = 0 + + +def need(ok, message): + if not ok: + raise ValueError(message) + + +def observe(x): + global MAX_HEIGHT + x = Q(x) + MAX_HEIGHT = max(MAX_HEIGHT, abs(x.numerator).bit_length(), x.denominator.bit_length()) + return x + + +class C: + __slots__ = ("re", "im") + + def __init__(self, re=0, im=0): + need(type(re) in (int, Q) and type(im) in (int, Q), "exact complex rational components") + self.re, self.im = observe(re), observe(im) + + @staticmethod + def wrap(x): + return x if isinstance(x, C) else C(x) + + def __add__(self, other): + other = C.wrap(other) + return C(self.re + other.re, self.im + other.im) + + __radd__ = __add__ + + def __neg__(self): + return C(-self.re, -self.im) + + def __sub__(self, other): + return self + -C.wrap(other) + + def __rsub__(self, other): + return C.wrap(other) + -self + + def __mul__(self, other): + other = C.wrap(other) + return C(self.re * other.re - self.im * other.im, + self.re * other.im + self.im * other.re) + + __rmul__ = __mul__ + + def conjugate(self): + return C(self.re, -self.im) + + def norm2(self): + return observe(self.re ** 2 + self.im ** 2) + + def __truediv__(self, other): + other = C.wrap(other) + denominator = other.norm2() + need(denominator != 0, "nonzero complex divisor") + product = self * other.conjugate() + return C(product.re / denominator, product.im / denominator) + + def __rtruediv__(self, other): + return C.wrap(other) / self + + def __bool__(self): + return bool(self.re or self.im) + + def __eq__(self, other): + if not isinstance(other, C) and type(other) not in (int, Q): + return False + other = C.wrap(other) + return self.re == other.re and self.im == other.im + + def record(self): + return [str(self.re), str(self.im)] + + +class Vector(list): + def __truediv__(self, scalar): + return Vector(C.wrap(x) / scalar for x in self) + + +class Matrix(list): + def __init__(self, rows): + super().__init__([[C.wrap(x) for x in row] for row in rows]) + + def __rmul__(self, scalar): + return Matrix([[C.wrap(x) * scalar for x in row] for row in self]) + + +def inner(left, right): + return sum((C.wrap(x).conjugate() * C.wrap(y) for x, y in zip(left, right)), C()) + + +def norm2(vector): + return observe(sum((C.wrap(x).norm2() for x in vector), Q(0))) + + +def component_l1(vector): + return observe(sum((abs(C.wrap(x).re) + abs(C.wrap(x).im) for x in vector), Q(0))) + + +def pow2(exponent): + return Q(1 << exponent) if exponent >= 0 else Q(1, 1 << -exponent) + + +def norm_scale(squared): + need(squared > 0, "positive residual squared norm") + exponent = (squared.numerator.bit_length() - squared.denominator.bit_length()) // 2 + scale = pow2(exponent) + while scale * scale < squared: + exponent += 1 + scale *= 2 + while scale * scale / 4 >= squared: + exponent -= 1 + scale /= 2 + need(squared <= scale * scale < 4 * squared, "dyadic norm bracket") + observe(scale) + return scale, exponent + + +def upper_dyadic(value): + need(value > 0, "positive dyadic upper target") + exponent = value.numerator.bit_length() - value.denominator.bit_length() + result = pow2(exponent) + while result < value: + result *= 2 + while result / 2 >= value: + result /= 2 + return observe(result) + + +def height(matrix): + return max([1] + [max(abs(x.numerator).bit_length(), x.denominator.bit_length()) + for row in matrix for entry in row + for x in (C.wrap(entry).re, C.wrap(entry).im)]) + + +def generators(displacement): + """Return column g, explicit coefficient ROW r, and a dyadic admission scale.""" + global MAX_HEIGHT + displacement = Matrix(displacement) + size = len(displacement) + need(size > 0 and all(len(row) == size for row in displacement), "nonempty square displacement") + input_height = height(displacement) + MAX_HEIGHT = input_height + columns = [[displacement[i][j] for i in range(size)] for j in range(size)] + residuals, result = [], [] + for column_index, column in enumerate(columns): + residual = list(column) + for old, squared in residuals: + factor = inner(old, residual) / squared + residual = [x - factor * y for x, y in zip(residual, old)] + squared = norm2(residual) + if squared == 0: + continue + scale, exponent = norm_scale(squared) + g = Vector(x / scale for x in residual) + row = Vector(inner(residual, c) * (scale / squared) for c in columns) + gamma = upper_dyadic(max(Q(1), 2 * component_l1(row))) + need(norm2(g) <= 1 and 2 * component_l1(row) <= gamma, "bounded generator admission") + residuals.append((residual, squared)) + result.append({"column": column_index, "residual": residual, "squared_norm": squared, + "scale": scale, "scale_exponent": exponent, "g": g, "row": row, + "gamma": gamma}) + return result, {"input_height": input_height, "max_observed_height": MAX_HEIGHT} + + +def lower(vector): + return Matrix([[vector[i - j] if i >= j else C() for j in range(len(vector))] + for i in range(len(vector))]) + + +def transpose(matrix): + return Matrix(zip(*matrix)) + + +def frobenius2(matrix): + return observe(sum((C.wrap(x).norm2() for row in matrix for x in row), Q(0))) + + +def norm_bound(matrix): + """Rational upper bound for the induced 2-norm.""" + return max(max(component_l1(row) for row in matrix), + max(component_l1(column) for column in zip(*matrix))) + + +def exact_record(matrix): + return [[C.wrap(x).record() for x in row] for row in matrix] diff --git a/research/rational-fft-fixed-h/vendor/rational_fft/engine.py b/research/rational-fft-fixed-h/vendor/rational_fft/engine.py new file mode 100644 index 0000000000..e9f2f2b08d --- /dev/null +++ b/research/rational-fft-fixed-h/vendor/rational_fft/engine.py @@ -0,0 +1,258 @@ +"""Certified phases and an authenticated rational adaptation of the FFT child. + +Jose Manuel (jmbarrancoml), with substantial OpenAI Codex assistance. +Apache-2.0. Phase projection follows Chafik/PR628; the schedule is Vlado/PR1673. +""" +import ast +import copy +from fractions import Fraction as Q +from hashlib import sha256 +from math import gcd, isqrt +from types import SimpleNamespace + +from exact import (C, Mat, Vec, adjoint, ceil_log2, component_l1, inverse_gate, + need, qheight, round_q, rounded) + + +class Phases: + def __init__(self, bits, source): + need(type(bits) is int and bits >= 1, "positive phase precision") + self.bits, self.nu = bits, Q(1, 1 << bits) + self.working_bits = 2 * (bits + 32) + self.oracle = {"__name__": "_authenticated_phase_oracle"} + exec(compile(source, "", "exec", optimize=0), self.oracle) + need(self.oracle["pi_interval"].cache_info().currsize == 0 + and self.oracle["kernel"].cache_info().currsize == 0, "fresh precision namespace") + self.oracle["BITS"] = self.working_bits + self.oracle["SCALE"] = 1 << self.working_bits + self.values, self.boxes, self.certificates = {}, {}, {} + + @staticmethod + def key(denominator, numerator): + need(type(denominator) is int and denominator > 0 and type(numerator) is int, + "integer phase description") + numerator %= denominator + divisor = gcd(denominator, numerator) + return denominator // divisor, numerator // divisor + + def box(self, denominator, numerator): + key = self.key(denominator, numerator) + if key not in self.boxes: + exact = {(1, 0): C(1), (2, 1): C(-1), (4, 1): C(0, 1), (4, 3): C(0, -1)} + if key in exact: + z = exact[key] + self.boxes[key] = ((z.re, z.re), (z.im, z.im)) + else: + self.boxes[key] = self.oracle["root_interval"](*key) + return self.boxes[key] + + def phase(self, denominator, numerator): + key = self.key(denominator, numerator) + if key in self.values: + return self.values[key] + box = self.box(*key) + if box[0][0] == box[0][1] and box[1][0] == box[1][1]: + sigma, error, parameter = C(box[0][0], box[1][0]), Q(0), None + else: + choices = [] + for center in (C(1), C(-1), C(0, 1), C(0, -1)): + conjugate = center.conjugate() + point = self.oracle["point"] + rotated = self.oracle["cmul"](box, (point(conjugate.re), point(conjugate.im))) + choices.append((rotated[0][0], center, rotated)) + _, center, rotated = max(choices, key=lambda item: item[0]) + need(rotated[0][0] >= Q(1, 2), "certified unit-phase chart") + parameters = [imag / (1 + real) for imag in rotated[1] for real in rotated[0]] + parameter = round_q((min(parameters) + max(parameters)) / 2, self.bits + 4) + need(abs(parameter) < 1, "bounded stereographic parameter") + sigma = center * C((1 - parameter ** 2) / (1 + parameter ** 2), + 2 * parameter / (1 + parameter ** 2)) + # Maximum endpoint distance, not minimum distance to a rectangle. + error = sum(max(abs(endpoint - value) for endpoint in interval) + for interval, value in zip(box, (sigma.re, sigma.im))) + need(sigma.norm2() == 1 and error <= self.nu, "certified exact rational unit phase") + need(max(qheight(sigma.re), qheight(sigma.im)) <= 2 * self.bits + 12, + "unit-phase representation bound") + self.values[key] = sigma + self.certificates[key] = {"denominator": key[0], "numerator": key[1], + "unit": sigma.record(), "error_bound": str(error), + "parameter": None if parameter is None else str(parameter)} + return sigma + + def record(self): + return {"phase_bits": self.bits, "oracle_bits": self.working_bits, + "certificates": [self.certificates[key] for key in sorted(self.certificates)]} + + +def split_kernel(matrix): + n = len(matrix) + need(n > 0 and n & (n - 1) == 0 and matrix.shape == (n, n), "square dyadic tile") + column, row = matrix[:, 0], matrix[0, :] + need(all(matrix[i][j] == (column[i - j] if i >= j else row[j - i]) + for i in range(n) for j in range(n)), "Toeplitz input") + g, h = Vec([0] * n), Vec([0] * n) + g[0] = h[0] = column[0] / 2 + for j in range(1, n): + g[j] = (column[j] + row[n - j]) / 2 + h[j] = (column[j] - row[n - j]) / 2 + return g, h + + +def budget(matrix, radius, target, input_fractional_bits=0): + g, h = split_kernel(matrix) + n, depth = len(matrix), len(matrix).bit_length() - 1 + mg, mh = component_l1(g), component_l1(h) + total = mg + mh + need(type(target) is int and target >= 1 and type(radius) in (int, Q) and radius >= 1, + "target and exact input radius") + need(type(input_fractional_bits) is int and input_fractional_bits >= 0, + "declared incoming fractional precision") + layers = 4 * depth + 8 + phase_bits = target + ceil_log2(16 * (depth + 1) * (total + 1) * (radius + 1)) + multiplier_bits = target + ceil_log2(8 * (radius + 1)) + data_bits = max(input_fractional_bits, target + ceil_log2(2 * layers * n * (total + 2))) + return certificate(n, mg, mh, radius, target, phase_bits, multiplier_bits, data_bits, + input_fractional_bits) + + +def certificate(n, mg, mh, radius, target, phase_bits, multiplier_bits, data_bits, + input_fractional_bits=0): + need(type(n) is int and n > 0 and n & (n - 1) == 0, "dyadic certificate dimension") + need(all(type(x) is int for x in (target, phase_bits, multiplier_bits, data_bits)) + and all(type(x) in (int, Q) for x in (mg, mh, radius)), "exact budget fields") + need(type(input_fractional_bits) is int and 0 <= input_fractional_bits <= data_bits, + "preserve incoming data grid") + need(min(phase_bits, multiplier_bits, data_bits, target) >= 1 and min(mg, mh) >= 0 + and radius >= 1, "positive precisions and valid envelopes") + depth, layers = n.bit_length() - 1, 4 * (n.bit_length() - 1) + 8 + nu, v, u = (Q(1, 1 << b) for b in (phase_bits, multiplier_bits, data_bits)) + total = mg + mh + coefficient_error = 4 * (depth + 1) * total * nu + 2 * v + prefix_bound = 1 + total + 2 * v + rounding_error = layers * prefix_bound * n * u + error = coefficient_error * radius + rounding_error + need(error <= Q(1, 1 << target), "complete target error budget") + need(layers * n * u <= 1, "bounded rounding accumulation") + amplitude = prefix_bound * (radius + layers * n * u) + integer_bits = max(1, ceil_log2(max(Q(1), amplitude))) + 2 + coefficient_height = max(4 * phase_bits + 24, + multiplier_bits + ceil_log2(total + 2) + 2) + temporary_bits = 2 * data_bits + 8 * coefficient_height + 2 * integer_bits + 32 + return dict(n=n, depth=depth, layer_bound=layers, mg=mg, mh=mh, radius=radius, + target_bits=target, phase_bits=phase_bits, multiplier_bits=multiplier_bits, + data_bits=data_bits, input_fractional_bits=input_fractional_bits, + coefficient_error=coefficient_error, prefix_bound=prefix_bound, + rounding_error=rounding_error, total_error=error, amplitude_bound=amplitude, + integer_bits=integer_bits, coefficient_height_bound=coefficient_height, + temporary_bits=temporary_bits) + + +class Adapt(ast.NodeTransformer): + def __init__(self): + self.function = None + self.changes = [] + + def visit_FunctionDef(self, node): + old = self.function + self.function = node.name + result = self.generic_visit(node) + self.function = old + return result + + def visit_Assign(self, node): + if len(node.targets) == 1 and isinstance(node.targets[0], ast.Name): + name = node.targets[0].id + if self.function == "fft" and name == "w": + self.changes.append("fft-twiddle") + node.value = ast.Call(func=ast.Name(id="phase", ctx=ast.Load()), + args=[ast.Name(id="size", ctx=ast.Load()), + ast.UnaryOp(op=ast.USub(), operand=ast.Name(id="j", ctx=ast.Load()))], + keywords=[]) + return node + if self.function == "toeplitz_shear" and name == "gamma": + self.changes.append("independent-twist-coefficients") + return ast.Pass() + if self.function == "toeplitz_shear" and name == "d": + self.changes.append("twist-vector") + node.value = ast.Call(func=ast.Name(id="certified_twist_vector", ctx=ast.Load()), + args=[ast.Name(id="n", ctx=ast.Load())], keywords=[]) + return node + return self.generic_visit(node) + + def visit_BinOp(self, node): + if (self.function == "fft" and isinstance(node.op, ast.Div) + and isinstance(node.right, ast.Call) and isinstance(node.right.func, ast.Attribute) + and node.right.func.attr == "sqrt"): + self.changes.append("unitary-splitter") + return ast.Call(func=ast.Name(id="butterfly", ctx=ast.Load()), + args=[ast.Name(id="w", ctx=ast.Load())], keywords=[]) + return self.generic_visit(node) + + +class Builder: + def __init__(self, phase_oracle, multiplier_bits, circuit_source): + self.phases, self.multiplier_bits = phase_oracle, multiplier_bits + self.spectra, self.twist_vectors, self.bounds = [], [], [] + selected = {"diag", "inverse", "parallel", "fft", "circ_shear", "toeplitz_shear"} + tree = ast.parse(circuit_source) + tree.body = [copy.deepcopy(n) for n in tree.body + if isinstance(n, ast.FunctionDef) and n.name in selected] + need({node.name for node in tree.body} == selected, "complete source function set") + self.original_ast = sha256(ast.dump(tree).encode()).hexdigest() + adapter = Adapt() + tree = ast.fix_missing_locations(adapter.visit(tree)) + need(adapter.changes == ["fft-twiddle", "unitary-splitter", + "independent-twist-coefficients", "twist-vector"], + "exact declared coefficient adaptation") + self.adapted_ast = sha256(ast.dump(tree).encode()).hexdigest() + self.np = SimpleNamespace(array=lambda rows, dtype=None: Mat(rows), + zeros=lambda n, dtype=None: Vec([0] * n), + fft=SimpleNamespace(fft=self.spectrum), + linalg=SimpleNamespace(inv=inverse_gate)) + self.namespace = {"np": self.np, "phase": self.phases.phase, + "butterfly": self.butterfly, "certified_twist_vector": self.twists} + exec(compile(tree, "", "exec", optimize=0), self.namespace) + + def butterfly(self, w): + splitter = self.phases.phase(8, 1) + a, b = splitter.re, splitter.im + return Mat([[a, b * w], [b, -a * w]]) + + def twists(self, n): + result = Vec(self.phases.phase(2 * n, j) for j in range(n)) + self.twist_vectors.append(result) + return result + + def spectrum(self, kernel): + n = len(kernel) + bound = self.bounds[len(self.spectra)] + exact = Vec(sum((kernel[j] * self.phases.phase(n, -i * j) for j in range(n)), C()) + for i in range(n)) + result = Vec(rounded(value, self.multiplier_bits) for value in exact) + v = Q(1, 1 << self.multiplier_bits) + need(all(value.norm2() <= bound ** 2 for value in exact), "unit-root multiplier bound") + need(all(value.norm2() <= (bound + v) ** 2 for value in result), "rounded multiplier bound") + need(all((x - y).norm2() <= v ** 2 for x, y in zip(result, exact)), "multiplier rounding") + self.spectra.append({"kernel": Vec(kernel), "multipliers": result, "bound": bound}) + return result + + def plan(self, matrix): + g, h = split_kernel(matrix) + n = len(matrix) + self.bounds = [component_l1(g), component_l1(h)] + self.spectra, self.twist_vectors = [], [] + result = self.namespace["toeplitz_shear"](list(range(n, 2 * n)), list(range(n)), matrix) + need(len(self.spectra) == 2 and len(self.twist_vectors) == 1, "two circulants and one twist table") + need(self.spectra[0]["kernel"] == g, "first source kernel") + need(self.spectra[1]["kernel"] == self.twist_vectors[0] * h, "twisted second source kernel") + need(all(a.norm2() == b.norm2() for a, b in zip(self.spectra[1]["kernel"], h)), + "twisting preserves kernel magnitudes") + return result + + +def inverse_sqrt_box(n, bits): + scale = 1 << bits + lower = isqrt((1 << (2 * bits)) // n) + if lower * lower * n == 1 << (2 * bits): + return Q(lower, scale), Q(lower, scale) + return Q(lower, scale), Q(lower + 1, scale) diff --git a/research/rational-fft-fixed-h/vendor/rational_fft/exact.py b/research/rational-fft-fixed-h/vendor/rational_fft/exact.py new file mode 100644 index 0000000000..aae41e6acf --- /dev/null +++ b/research/rational-fft-fixed-h/vendor/rational_fft/exact.py @@ -0,0 +1,242 @@ +"""Exact complex rational arithmetic for coefficient and rounded-gate certificates. + +Jose Manuel (jmbarrancoml), with substantial OpenAI Codex assistance. +Apache-2.0. +""" +from fractions import Fraction as Q +from hashlib import sha256 + +TRACK = None + + +def need(ok, message): + if not ok: + raise ValueError(message) + + +def qheight(value): + return max(abs(value.numerator).bit_length(), value.denominator.bit_length()) + + +class C: + __slots__ = ("re", "im") + + def __init__(self, re=0, im=0): + need(type(re) in (int, Q) and type(im) in (int, Q), "exact rational components") + self.re, self.im = Q(re), Q(im) + if TRACK is not None: + TRACK["max_observed_height"] = max(TRACK["max_observed_height"], + qheight(self.re), qheight(self.im)) + + @staticmethod + def wrap(value): + return value if isinstance(value, C) else C(value) + + def __add__(self, other): + other = C.wrap(other) + return C(self.re + other.re, self.im + other.im) + + __radd__ = __add__ + + def __neg__(self): + return C(-self.re, -self.im) + + def __sub__(self, other): + return self + -C.wrap(other) + + def __rsub__(self, other): + return C.wrap(other) + -self + + def __mul__(self, other): + other = C.wrap(other) + return C(self.re * other.re - self.im * other.im, + self.re * other.im + self.im * other.re) + + __rmul__ = __mul__ + + def conjugate(self): + return C(self.re, -self.im) + + def norm2(self): + return self.re * self.re + self.im * self.im + + def __truediv__(self, other): + other = C.wrap(other) + denominator = other.norm2() + need(denominator > 0, "nonzero divisor") + product = self * other.conjugate() + return C(product.re / denominator, product.im / denominator) + + def __rtruediv__(self, other): + return C.wrap(other) / self + + def __eq__(self, other): + if not isinstance(other, C) and type(other) not in (int, Q): + return False + other = C.wrap(other) + return self.re == other.re and self.im == other.im + + def record(self): + return [str(self.re), str(self.im)] + + +class Vec(list): + def __init__(self, values): + super().__init__(C.wrap(x) for x in values) + + def __mul__(self, other): + if isinstance(other, C) or type(other) in (int, Q): + return Vec(x * other for x in self) + need(len(self) == len(other), "elementwise vector dimensions") + return Vec(x * y for x, y in zip(self, other)) + + __rmul__ = __mul__ + + def __sub__(self, other): + need(len(self) == len(other), "vector subtraction dimensions") + return Vec(x - y for x, y in zip(self, other)) + + def copy(self): + return Vec(self) + + def __truediv__(self, scalar): + return Vec(x / scalar for x in self) + + +class Mat(list): + def __init__(self, rows): + super().__init__([C.wrap(x) for x in row] for row in rows) + + def __rmul__(self, scalar): + return Mat([[scalar * x for x in row] for row in self]) + + @property + def shape(self): + return (len(self), len(list.__getitem__(self, 0)) if self else 0) + + def __getitem__(self, key): + if not isinstance(key, tuple): + return list.__getitem__(self, key) + row, column = key + if isinstance(row, slice): + rows = list.__getitem__(self, row) + return Mat([r[column] for r in rows]) if isinstance(column, slice) else Vec(r[column] for r in rows) + selected = list.__getitem__(self, row) + return Vec(selected[column]) if isinstance(column, slice) else selected[column] + + +def eye(n): + return Mat([[int(i == j) for j in range(n)] for i in range(n)]) + + +def adjoint(matrix): + return Mat([[matrix[j][i].conjugate() for j in range(len(matrix))] + for i in range(len(matrix[0]))]) + + +def multiply(left, right): + return Mat([[sum((x * right[k][j] for k, x in enumerate(row)), C()) + for j in range(len(right[0]))] for row in left]) + + +def unitary(matrix): + return multiply(adjoint(matrix), matrix) == eye(len(matrix)) + + +def inverse_gate(matrix): + if unitary(matrix): + return adjoint(matrix) + if matrix.shape == (1, 1): + return Mat([[1 / matrix[0][0]]]) + if matrix.shape == (2, 2) and matrix[0][0] == matrix[1][1] == 1 and matrix[0][1] == 0: + return Mat([[1, 0], [-matrix[1][0], 1]]) + raise ValueError("unsupported nonunitary inverse gate") + + +def round_q(value, bits): + scale = 1 << bits + shifted = value * scale + Q(1, 2) + return Q(shifted.numerator // shifted.denominator, scale) + + +def rounded(value, bits): + return C(round_q(value.re, bits), round_q(value.im, bits)) + + +def norm2(vector): + return sum((x.norm2() for x in vector), Q(0)) + + +def component_l1(vector): + return sum((abs(x.re) + abs(x.im) for x in vector), Q(0)) + + +def ceil_log2(value): + value = Q(value) + need(value > 0, "positive logarithm bound") + result = value.numerator.bit_length() - value.denominator.bit_length() + power = Q(1 << result) if result >= 0 else Q(1, 1 << -result) + while power < value: + result += 1 + power *= 2 + while power / 2 >= value: + result -= 1 + power /= 2 + return result + + +def upper(value, bits=160): + scale = 1 << bits + scaled = value * scale + return Q(-((-scaled.numerator) // scaled.denominator), scale) + + +def vector_hash(vector): + digest = sha256() + for entry in vector: + for value in (entry.re, entry.im): + for integer in (value.numerator, value.denominator): + raw = abs(integer).to_bytes(max(1, (abs(integer).bit_length() + 7) // 8), "big") + digest.update(bytes([integer < 0])) + digest.update(len(raw).to_bytes(8, "big")) + digest.update(raw) + return digest.hexdigest() + + +def apply(plan, vector, bits=None): + global TRACK + state = Vec(vector) + old_track = TRACK + stats = {"max_observed_height": 0, "max_grid_numerator_bits": 0, + "max_reduced_fractional_bits": 0, "rounded_outputs": 0} + if bits is not None: + need(bits >= 1, "positive data precision") + for entry in state: + for value in (entry.re, entry.im): + need((value * (1 << bits)).denominator == 1, "input already on data grid") + TRACK = stats + try: + for layer in plan: + used = set() + for ids, gate in layer: + need(len(ids) == len(gate) in (1, 2), "one- or two-wire gate") + need(not used.intersection(ids) and len(set(ids)) == len(ids), "disjoint layer") + need(all(0 <= i < len(state) for i in ids), "wire bounds") + used.update(ids) + before = [state[i] for i in ids] + outputs = [sum((a * x for a, x in zip(row, before)), C()) for row in gate] + for i, value in zip(ids, outputs): + if bits is not None: + value = rounded(value, bits) + stats["rounded_outputs"] += 1 + for part in (value.re, value.im): + grid_integer = part * (1 << bits) + need(grid_integer.denominator == 1, "rounded output grid") + stats["max_grid_numerator_bits"] = max(stats["max_grid_numerator_bits"], + abs(grid_integer.numerator).bit_length()) + stats["max_reduced_fractional_bits"] = max(stats["max_reduced_fractional_bits"], + part.denominator.bit_length() - 1) + state[i] = value + finally: + TRACK = old_track + return state, stats diff --git a/research/rational-fft-fixed-h/verify.py b/research/rational-fft-fixed-h/verify.py new file mode 100644 index 0000000000..7637cc85f6 --- /dev/null +++ b/research/rational-fft-fixed-h/verify.py @@ -0,0 +1,308 @@ +"""Check retained fixed-H words, routing and complete rounded FFT executions. + +Jose Manuel (jmbarrancoml), with substantial OpenAI Codex assistance. +Apache-2.0. This is finite evidence for the separately written proof. +""" +import argparse +from fractions import Fraction as Q +from hashlib import sha256 +import json +from pathlib import Path +import sys + +HERE = Path(__file__).resolve().parent +sys.path.insert(0, str(HERE)) +from lowering import FFT, K, Sources, TENSOR, canonical, need, norm2, stage_hash + + +def rejected(operation, label): + try: + operation() + except (ValueError, AssertionError): + return label + raise ValueError("control unexpectedly passed: " + label) + + +def exact_native(source, stages, vector, raw=False): + E = source.e + state = E.Vec(vector) + for kind, value in stages: + if kind == "D": + need(len(value) == len(state), "complete coefficient page") + state = E.Vec(a * x for a, x in zip(value, state)) + else: + step = 1 << value + for i in range(len(state)): + if not i & step: + a, b = state[i], state[i + step] + denominator = 1 if raw else 2 + state[i], state[i + step] = (a + b) / denominator, (a - b) / denominator + return state + + +def difference(a, b): + need(len(a) == len(b), "complete output comparison") + return [x - y for x, y in zip(a, b)] + + +def upper_report(source, value, allowance): + need(0 <= value <= allowance, "reported exact error allowance") + if allowance == 0: + return "0" + bits = max(160, source.e.ceil_log2(1 / allowance) + 16) + return str(min(allowance, source.e.upper(value, bits))) + + +def scale_interval(interval, scalar): + ends = [x * scalar for x in interval] + return min(ends), max(ends) + + +def add_interval(a, b): + return a[0] + b[0], a[1] + b[1] + + +def ideal_squared_error(program, vector, output): + source, indices = program.source, program.indices + n = len(indices) + normalization = source.engine.inverse_sqrt_box(n, program.phases.working_bits) + boxes = {i: ((z.re, z.re), (z.im, z.im)) for i, z in enumerate(vector)} + for row, physical_row in enumerate(indices): + real, imaginary = (Q(0), Q(0)), (Q(0), Q(0)) + for col, physical_col in enumerate(indices): + phase = program.phases.box(n, (1 if program.adjoint else -1) * row * col) + z = vector[physical_col] + real = add_interval(real, add_interval(scale_interval(phase[0], z.re), + scale_interval(phase[1], -z.im))) + imaginary = add_interval(imaginary, add_interval(scale_interval(phase[0], z.im), + scale_interval(phase[1], z.re))) + intervals = [] + for part in (real, imaginary): + corners = [a * b for a in part for b in normalization] + intervals.append((min(corners), max(corners))) + boxes[physical_row] = tuple(intervals) + error = Q(0) + for i, z in enumerate(output): + for value, interval in zip((z.re, z.im), boxes[i]): + error += max(abs(value - interval[0]), abs(value - interval[1])) ** 2 + return error + + +def basis_check(program): + E, source = program.e, program.source + columns = [] + for j in range(program.size): + vector = E.Vec(int(i == j) for i in range(program.size)) + expected = E.apply(program.plan, vector)[0] + actual = exact_native(source, program.reference_stages, vector) + need(actual == expected, "complete exact native matrix equals retained FFT") + columns.append(actual) + return E.vector_hash(z for column in columns for z in column) + + +def one_case(source, label, size, indices, target=18, incoming=0, radius=1, adjoint=False): + incoming = max(incoming, size.bit_length() + 2) + program = FFT(source, size, indices, radius, target, incoming, adjoint=adjoint) + E = source.e + vector = E.Vec(E.C(Q((3 * i) % 7 - 3, 8 * size), Q((5 * i) % 7 - 3, 8 * size)) + for i in range(size)) + if radius >= 2: + vector = E.Vec(z * (Q(radius) / 2) for z in vector) + if incoming >= 512: + fine = next((i for i in range(size) if i not in indices), 0) + vector[fine] = E.C(Q(1, 1 << incoming)) + expected = E.apply(program.plan, vector)[0] + exact = exact_native(source, program.reference_stages, vector) + need(exact == expected, "exact complete forward or stored-adjoint lowering") + output, statistics, trace = program.execute(vector, capture=True) + rational_error = norm2(difference(output, expected)) + need(rational_error <= program.native_error ** 2, "complete native error against rational FFT") + ideal_error = ideal_squared_error(program, vector, output) + need(ideal_error <= program.total_error ** 2, "complete independent-formula Fourier enclosure") + matrix_hash = basis_check(program) if size <= 8 else None + spectators = [i for i in range(size) if i not in indices] + spectator_error = norm2([output[i] - vector[i] for i in spectators]) + record = program.record() + record.update(label=label, input_sha256=E.vector_hash(vector), output_sha256=E.vector_hash(output), + rational_target_sha256=E.vector_hash(expected), exact_matrix_sha256=matrix_hash, + rounded_rational_error_squared_upper=upper_report(source, rational_error, program.native_error ** 2), + ideal_error_squared_upper=upper_report(source, ideal_error, program.total_error ** 2), + nonzero_rounded_error=rational_error > 0, + nonzero_spectator_error=spectator_error > 0, + spectator_error_squared_upper=upper_report(source, spectator_error, program.native_error ** 2), + statistics=statistics, trace=trace) + print("CHECK_FIXED_H_FFT", label, flush=True) + return record, program, vector, output + + +def segment_envelopes(source, word): + E, matrices = source.e, [] + for i, operation in enumerate(word): + if operation == "H": + matrices.append(E.Mat([[Q(1, 2), Q(1, 2)], [Q(1, 2), Q(-1, 2)]])) + else: + scale = 2 if i and word[i - 1] == "H" else 1 + matrices.append(E.Mat([[scale * operation[0], 0], [0, scale * operation[1]]])) + forward, inverse, count = Q(0), Q(0), 0 + for first in range(len(matrices)): + product = E.Mat([[1, 0], [0, 1]]) + for matrix in matrices[first:]: + product = E.Mat([[sum((matrix[i][k] * product[k][j] for k in range(2)), E.C()) + for j in range(2)] for i in range(2)]) + a, b = product[0] + c, d = product[1] + determinant = a * d - b * c + need(determinant != 0, "invertible exact native segment") + inv = E.Mat([[d / determinant, -b / determinant], [-c / determinant, a / determinant]]) + forward = max(forward, sum(z.norm2() for row in product for z in row)) + inverse = max(inverse, sum(z.norm2() for row in inv for z in row)) + count += 1 + need(max(forward, inverse) <= 4096 ** 2, "independent Frobenius segment envelope") + return dict(segments=count, forward_frobenius_squared_upper=str(E.upper(forward)), + inverse_frobenius_squared_upper=str(E.upper(inverse))) + + +def word_cases(source): + E, records = source.e, [] + for a, b in ((Q(3, 5), Q(4, 5)), (Q(33, 65), Q(56, 65))): + for w in (E.C(1), E.C(-1), E.C(0, 1), E.C(3, 4) / 5): + for inverse in (False, True): + word = source.butterfly_word(a, b, w, inverse) + wanted = E.Mat([[a, b * w], [b, -a * w]]) + if inverse: + wanted = E.adjoint(wanted) + need(source.word_endpoint(word) == wanted, "standalone six-H butterfly identity") + coefficients = [z for op in word if op != "H" for z in op] + need(all(z.norm2() <= 9 and (1 / z).norm2() <= 400 ** 2 for z in coefficients), + "raw direct and reciprocal coefficient bounds") + records.append(dict(a=str(a), b=str(b), phase=w.record(), adjoint=inverse, + endpoint_sha256=source.matrix_hash(wanted), + envelopes=segment_envelopes(source, word))) + return records + + +def controls(source, cases): + E, results = source.e, [] + program = next(row[1] for row in cases if row[1].size == 4 and not row[1].adjoint + and len(row[1].indices) == 4) + row = next(row for row in cases if row[1] is program) + vector = row[2] + expected = E.apply(program.plan, vector)[0] + wrong = exact_native(source, program.raw_stages, vector) + results.append(rejected(lambda: need(wrong == expected, "missing H0 compensation"), + "omitted-native-scale-restoration")) + a, b, w = Q(3, 5), Q(4, 5), E.C(3, 4) / 5 + actual = E.Mat([[a, b * w], [b, -a * w]]) + old_word = source.ns["gateword"]("fft", w) + results.append(rejected(lambda: need(source.word_endpoint(old_word) == actual, "old FFT shape"), + "unnormalized-source-fft-word-is-not-rational-butterfly")) + wrong_adjoint = source.butterfly_word(a, b, w.conjugate()) + results.append(rejected(lambda: need(source.word_endpoint(wrong_adjoint) == E.adjoint(actual), "adjoint orientation"), + "twiddle-conjugation-alone-is-not-the-adjoint")) + pairs = [(3, 1), (0, 2)] + words = [source.butterfly_word(a, b, w), source.ns["gateword"]("diag", (1, 1))] + routed = source.ns["local_layer"](4, pairs, words) + expected_partial = E.apply([[(pairs[0], actual)]], vector)[0] + need(exact_native(source, routed, vector, raw=True) == expected_partial, "oriented noninvolutive route") + sorted_route = source.ns["local_layer"](4, [tuple(sorted(p)) for p in pairs], words) + results.append(rejected(lambda: need(exact_native(source, sorted_route, vector, raw=True) + == expected_partial, "sorted pair changes gate"), "discarded-pair-orientation")) + all_one = [(E.C(1), E.C(1)) if op != "H" else op for op in words[1]] + bad_idle = source.ns["local_layer"](4, pairs, [words[0], all_one]) + results.append(rejected(lambda: need(exact_native(source, bad_idle, vector, raw=True) + == expected_partial, "idle word needs compensation"), "all-one-idle-padding")) + missing_half = list(routed) + half_index = next(i for i, (kind, value) in enumerate(missing_half) + if kind == "D" and all(z == Q(1, 2) for z in value)) + missing_half[half_index] = ("D", E.Vec([1] * 4)) + results.append(rejected(lambda: need(exact_native(source, missing_half, vector, raw=True) + == expected_partial, "switch half"), "omitted-switch-half")) + results.append(rejected(lambda: source.ns["u3"](E.C(0)), "source-zero-shear-assertion-active")) + results.append(rejected(lambda: source.ns["gateword"]("diag", (0, 1)), "source-zero-diagonal-assertion-active")) + results.append(rejected(lambda: FFT(source, 4, [0, 0], 1, 18), "duplicate-active-wire")) + results.append(rejected(lambda: FFT(source, 4, [False, 1], 1, 18), "boolean-wire-id")) + results.append(rejected(lambda: FFT(source, 4, [-1, 1], 1, 18), "negative-wire-id")) + results.append(rejected(lambda: FFT(source, 4, [0, 4], 1, 18), "out-of-range-wire-id")) + results.append(rejected(lambda: FFT(source, 3, [0, 1], 1, 18), "nondyadic-field")) + results.append(rejected(lambda: FFT(source, 4, [0, 1, 2], 1, 18), "nondyadic-active-FFT")) + results.append(rejected(lambda: FFT(source, 4, [0, 1], 1.0, 18), "floating-radius")) + results.append(rejected(lambda: FFT(source, 4, [0, 1], 1, 18, phase_bits=3), "underfunded-phase-precision")) + results.append(rejected(lambda: FFT(source, 4, [0, 1], 1, 18, data_bits=8), "underfunded-native-grid")) + bad_input = E.Vec(vector) + bad_input[0] = E.C(Q(1, 1 << (program.incoming + 1))) + results.append(rejected(lambda: program.execute(bad_input), "underdeclared-incoming-grid")) + bad_input[0] = E.C(2) + results.append(rejected(lambda: program.execute(bad_input), "underdeclared-input-radius")) + first_diagonal = next(i for i, (kind, _) in enumerate(program.stages) if kind == "D") + old_page = program.stages[first_diagonal] + replacement = E.Vec(old_page[1]) + replacement[0] += 1 + program.stages[first_diagonal] = ("D", replacement) + results.append(rejected(lambda: program.execute(vector), "changed-prepared-coefficient-page")) + program.stages[first_diagonal] = old_page + results.append(rejected(lambda: Sources({TENSOR: source.raw[TENSOR] + b"\n# changed\n"}), + "source-byte-corruption")) + fine = next(row for row in cases if row[0]["label"] == "spectator-fine") + need(fine[0]["nonzero_spectator_error"], "inactive rounding has a nonzero error witness") + results.append("inactive-word-error-is-paid") + need(canonical({"count": 1}) != canonical({"count": True}), "typed scientific receipt") + return results + ["typed-receipt"] + + +def run(quick=False): + source = Sources() + specifications = [ + ("singleton", 1, [0], 18, 0, 1, False), + ("two-forward", 2, [0, 1], 18, 0, 1, False), + ("two-adjoint", 2, [0, 1], 18, 0, 1, True), + ("spectator-fine", 4, [3, 1], 18, 512, 1, False), + ("four-forward", 4, [0, 1, 2, 3], 18, 0, 1, False), + ("four-reversed-adjoint", 4, [3, 2, 1, 0], 24, 0, 1, True), + ("eight-permuted", 8, [5, 1, 7, 3, 0, 6, 2, 4], 22, 0, 1, False), + ("eight-partial-adjoint", 8, [6, 1, 4, 3], 20, 0, 1, True), + ("sixteen-forward", 16, list(range(16)), 24, 0, 1, False), + ("thirtytwo-adjoint", 32, list(reversed(range(32))), 20, 0, 1, True), + ("sixtyfour-forward", 64, list(range(64)), 18, 0, 1, False), + ("large-radius", 8, list(range(8)), 18, 0, 1 << 20, False), + ("idle-field", 4, [2], 18, 0, 1, True), + ] + if quick: + specifications = specifications[:4] + cases = [one_case(source, *spec) for spec in specifications] + if quick: + return dict(cases=[row[0] for row in cases]) + need(any(row[0]["nonzero_rounded_error"] for row in cases), "nonzero native rounding evidence") + return dict(schema_version=1, scope="retained rational FFT and stored-adjoint fixed-H stage lowering", + sources=source.seals, tensor_ast_sha256=source.tensor_ast, + word_cases=word_cases(source), cases=[row[0] for row in cases], + controls=controls(source, cases), + source_FFT_execution=True, source_Benes_execution=True, native_stage_arithmetic=True, + physical_root_execution=False, physical_tape_execution=False, + complete_tensor_execution=False, dyadic_programs_exact_inverse_pair=False, + independent_transcendental_backend=False, universal_proof_kernel=False) + + +def main(): + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--write-result", action="store_true") + parser.add_argument("--quick", action="store_true") + args = parser.parse_args() + need(not (args.quick and args.write_result), "diagnostic replay cannot replace the full receipt") + result = run(args.quick) + if args.quick: + print(json.dumps(dict(status="PASS_FIXED_H_DIAGNOSTIC", cases=len(result["cases"])))) + return + path = HERE / "RESULT.json" + if args.write_result: + path.write_text(json.dumps(result, indent=2) + "\n") + else: + need(canonical(result) == canonical(json.loads(path.read_text())), "complete typed native FFT receipt") + print(json.dumps(dict(status="PASS_RATIONAL_FFT_FIXED_H", cases=len(result["cases"]), + word_cases=len(result["word_cases"]), controls=len(result["controls"]), + native_stages=sum(row["counts"]["primitive_stages"] for row in result["cases"]), + rounded_outputs=sum(row["statistics"]["complex_roundings"] for row in result["cases"]), + result_sha256=sha256(canonical(result)).hexdigest()), sort_keys=True)) + + +if __name__ == "__main__": + main() diff --git a/research/rational-fft-fixed-h/verify_dense.py b/research/rational-fft-fixed-h/verify_dense.py new file mode 100644 index 0000000000..f291736487 --- /dev/null +++ b/research/rational-fft-fixed-h/verify_dense.py @@ -0,0 +1,323 @@ +"""Check complete conditioned Toeplitz programs on the fixed-H grammar. + +Jose Manuel (jmbarrancoml), with substantial OpenAI Codex assistance. +Apache-2.0. The retained compiler and source identities belong to their +credited authors. Finite checks support the separate written argument. +""" +import argparse +from fractions import Fraction as Q +from hashlib import sha256 +import json +from pathlib import Path +import sys + +HERE = Path(__file__).resolve().parent +sys.path.insert(0, str(HERE)) +from dense import DenseSources, FullRoot, NativeProgram +from lowering import DENSE_PREFIX, canonical, need, norm2 +from verify import basis_check, difference, exact_native, rejected, upper_report + + +def identity(E, n): + return E.Mat([[int(i == j) for j in range(n)] for i in range(n)]) + + +def multiply(E, a, b): + return E.Mat([[sum((a[i][k] * b[k][j] for k in range(len(b))), E.C()) + for j in range(len(b[0]))] for i in range(len(a))]) + + +def inverse(E, matrix): + """Separate exact elimination, with both inverse identities checked.""" + n = len(matrix) + rows = [list(row) + list(identity(E, n)[i]) for i, row in enumerate(matrix)] + for column in range(n): + pivot = next(i for i in range(column, n) if rows[i][column] != 0) + rows[column], rows[pivot] = rows[pivot], rows[column] + value = rows[column][column] + rows[column] = [z / value for z in rows[column]] + for i in range(n): + if i != column: + value = rows[i][column] + rows[i] = [a - value * b for a, b in zip(rows[i], rows[column])] + result = E.Mat([row[n:] for row in rows]) + need(multiply(E, matrix, result) == multiply(E, result, matrix) == identity(E, n), + "two exact inverse reference identities") + return result + + +def apply_matrix(E, matrix, vector, indices): + result = E.Vec(vector) + for i, physical in enumerate(indices): + result[physical] = sum((matrix[i][j] * vector[address] + for j, address in enumerate(indices)), E.C()) + return result + + +def operator(E, plan, n): + columns = [E.apply(plan, E.Vec(int(i == j) for i in range(n)))[0] for j in range(n)] + return E.Mat(zip(*columns)) + + +def check_execution(program, vector, reference): + E = program.e + stored = E.apply(program.plan, vector)[0] + output, statistics, trace = program.execute(vector) + need(not trace, "compact execution omits per-stage hashes") + native_error = norm2(difference(output, stored)) + upstream_error = norm2(difference(stored, reference)) + total_error = norm2(difference(output, reference)) + need(native_error <= program.native_error ** 2, "new rounded execution against stored logical program") + need(upstream_error <= (program.upstream_error * program.radius) ** 2, + "source coefficient defect against same exact root") + need(total_error <= program.total_error ** 2, "complete native-to-root error") + return dict(input_sha256=E.vector_hash(vector), stored_output_sha256=E.vector_hash(stored), + output_sha256=E.vector_hash(output), reference_sha256=E.vector_hash(reference), + native_squared_error_upper=upper_report(program.source, native_error, program.native_error ** 2), + upstream_squared_error_upper=upper_report(program.source, upstream_error, + (program.upstream_error * program.radius) ** 2), + total_squared_error_upper=upper_report(program.source, total_error, program.total_error ** 2), + nonzero_native_error=native_error > 0, statistics=statistics) + + +def full_case(source, label, matrix, reverse=False, fine=False): + E, n = source.e, len(matrix) + indices = list(reversed(range(n))) if reverse else list(range(n)) + root = FullRoot(source, matrix, 4, 18, indices) + inverse_matrix = inverse(E, matrix) + incoming = 512 if fine else n.bit_length() + 6 + vector = E.Vec(E.C(Q((3 * i) % 7 - 3, 8 * n), Q((5 * i) % 7 - 2, 8 * n)) + for i in range(n)) + if fine: + vector[indices[0]] += E.C(Q(1, 1 << incoming)) + record = dict(label=label, **root.record(), directions=[]) + for backward, reference_matrix in ((False, matrix), (True, inverse_matrix)): + program = root.lower(backward, incoming) + exact = apply_matrix(E, reference_matrix, vector, indices) + row = dict(inverse=backward, **program.record(), execution=check_execution(program, vector, exact)) + if n <= 2: + row["complete_exact_basis_sha256"] = basis_check(program) + record["directions"].append(row) + print("CHECK_NATIVE_FULL", label, flush=True) + return record + + +def logical_envelopes(source, plan, n): + E, maximum = source.e, Q(1) + for stop in range(len(plan) + 1): + matrix = operator(E, plan[:stop], n) + maximum = max(maximum, source.full.norm_upper(matrix), source.full.norm_upper(inverse(E, matrix))) + return maximum, max(Q(1), source.full.norm_upper(operator(E, plan, n))) + + +def generic_cases(source): + E = source.e + specifications = [ + ("mixed-oriented", 4, [[((3, 0), E.Mat([[1, 0], [E.C(1, 2), 1]])), + ((2,), E.Mat([[E.C(2, -1)]]))], []]), + ("zero-and-tiny", 2, [[((0, 1), E.Mat([[1, 0], [0, 1]]))], + [((1, 0), E.Mat([[1, 0], [E.C(Q(1, 1 << 80), Q(-1, 1 << 81)), 1]]))]]), + ("one-wire-idle", 1, [[]]), + ("empty-program", 4, []), + ] + rows = [] + for label, n, plan in specifications: + prefix, endpoint = logical_envelopes(source, plan, n) + program = NativeProgram(source, n, plan, prefix, endpoint, 1, 18, 512) + vector = E.Vec(E.C(Q(i + 1, 16 * n), Q(-1, 32 * n)) for i in range(n)) + vector[-1] += E.C(Q(1, 1 << 512)) + reference = E.apply(plan, vector)[0] + need(exact_native(source, program.reference_stages, vector) == reference, + "complete generic native endpoint") + rows.append(dict(label=label, **program.record(), complete_exact_basis_sha256=basis_check(program), + execution=check_execution(program, vector, reference))) + return rows + + +def segment_cases(source): + E, records = source.e, [] + for gamma in (1, 16, 1 << 32): + values = (E.C(0), E.C(gamma), E.C(0, -gamma), E.C(Q(1, 1 << 80), Q(-1, 1 << 81))) + for kind in ("upper", "lower"): + for value in values: + word = source.ns["gateword"](kind, value, gamma) + wanted = E.Mat([[1, value], [0, 1]]) if kind == "upper" else E.Mat([[1, 0], [value, 1]]) + need(source.word_endpoint(word) == wanted, "zero-safe exact shear endpoint") + matrices = [] + for i, op in enumerate(word): + if op == "H": + matrices.append(E.Mat([[Q(1, 2), Q(1, 2)], [Q(1, 2), Q(-1, 2)]])) + else: + scale = 2 if i and word[i - 1] == "H" else 1 + matrices.append(E.Mat([[scale * op[0], 0], [0, scale * op[1]]])) + direct, reciprocal, count = Q(0), Q(0), 0 + for start in range(len(matrices)): + product = identity(E, 2) + for matrix in matrices[start:]: + product = multiply(E, matrix, product) + direct = max(direct, source.full.norm_upper(product)) + reciprocal = max(reciprocal, source.full.norm_upper(inverse(E, product))) + count += 1 + need(max(direct, reciprocal) <= 4096 * gamma, "native shear segments and inverse segments") + records.append(dict(kind=kind, gamma=str(gamma), value=value.record(), segments=count, + direct_component_envelope=str(direct), inverse_component_envelope=str(reciprocal))) + return records + + +def endpoint_control(source): + E, q = source.e, 1 << 16 + plan = [[((0, 1), E.Mat([[Q(1, q), 0], [0, q]]))], + [((0, 1), identity(E, 2))], + [((0, 1), E.Mat([[q * q, 0], [0, Q(1, q * q)]]))]] + prefix, endpoint = logical_envelopes(source, plan, 2) + need(prefix == endpoint == q, "exact scaled logical prefix and endpoint bounds") + program = NativeProgram(source, 2, plan, prefix, endpoint, 1, 18, 8) + stop = program.layers[1]["begin"] + 5 + columns = [exact_native(source, program.reference_stages[:stop], E.Vec(int(i == j) for i in range(2))) + for j in range(2)] + partial = E.Mat(zip(*columns)) + c = 2 * q * q + 1 + need(partial == E.Mat([[Q(-5, 4 * c * q), 0], [0, q]]), "actual zero-shear intermediate prefix") + suffix = multiply(E, operator(E, plan, 2), inverse(E, partial)) + suffix_bound = source.full.norm_upper(suffix) + need(suffix_bound == Q(4 * c * q * q, 5) > program.compiled_bound, + "unit-endpoint suffix allowance fails on admitted nonunitary program") + need(suffix_bound <= endpoint * program.compiled_bound, "corrected endpoint suffix allowance") + vector = E.Vec([Q(1, 16), Q(-1, 32)]) + return dict(q=q, logical_prefix=str(prefix), endpoint=str(endpoint), native_prefix=str(program.compiled_bound), + intermediate_sha256=source.matrix_hash(partial), suffix_norm=str(suffix_bound), + ratio_to_invalid_allowance=str(suffix_bound / program.compiled_bound), + execution=check_execution(program, vector, E.apply(plan, vector)[0])) + + +def controls(source): + E, rows = source.e, [] + rows.append(rejected(lambda: FullRoot(source, E.Mat([[0, 1], [1, 0]]), 1, 18), "singular-leading-section")) + rows.append(rejected(lambda: FullRoot(source, E.Mat([[1, 2], [3, 4]]), 1, 18), "non-Toeplitz-root")) + root = FullRoot(source, E.Mat([[E.C(Q(1, 3), Q(1, 7))]]), 1, 18) + forward, backward = root.lower(False, 8), root.lower(True, 8) + vector = E.Vec([Q(1, 16)]) + a = root.plan[0][0][1][0][0] + need(root.inverse_plan[0][0][1][0][0] == 1 / a and a.conjugate() != 1 / a, + "stored reciprocal differs from global adjoint") + rows.append("stored-inverse-requires-reciprocals") + rows.append(rejected(lambda: forward.execute(E.Vec([Q(1, 1 << 9)])), "underdeclared-incoming-grid")) + rows.append(rejected(lambda: forward.execute(E.Vec([2])), "underdeclared-radius")) + rows.append(rejected(lambda: root.lower(False, True), "boolean-incoming-grid")) + rows.append(rejected(lambda: NativeProgram(source, 2, [[((0, 1), E.Mat([[1, 2], [3, 4]]))]], + 16, 16, 1, 18), "unsupported-two-wire-gate")) + rows.append(rejected(lambda: NativeProgram(source, 2, [[((0, 0), identity(E, 2))]], + 1, 1, 1, 18), "duplicate-wire-ownership")) + rows.append(rejected(lambda: NativeProgram(source, 2, [[]], 1, 1, 1, 18, data_bits=8), "underfunded-native-grid")) + rows.append(rejected(lambda: NativeProgram(source, 1, [[((0,), E.Mat([[0]]))]], 1, 1, 1, 18), "zero-scalar")) + rows.append(rejected(lambda: source.ns["u3"](E.C(0)), "unshifted-zero-shear-division")) + z = E.C(1, 2) + lower = source.ns["gateword"]("lower", z, 4) + wrong = [op if op == "H" else tuple(a.conjugate() for a in op) for op in lower] + need(source.word_endpoint(wrong) != source.word_endpoint(lower), "conjugate transpose changes complex lower shear") + rows.append("ordinary-transpose-is-required") + old = root.matrix[0][0] + root.matrix[0][0] += 1 + rows.append(rejected(lambda: forward.execute(vector), "changed-admitted-root")) + root.matrix[0][0] = old + old = root.builder.bound["execution_bound"] + root.builder.bound["execution_bound"] += 1 + rows.append(rejected(lambda: backward.execute(vector), "changed-full-prefix-contract")) + root.builder.bound["execution_bound"] = old + for name in ("lam", "q", "radius", "target"): + old = getattr(root, name) + setattr(root, name, old + 1) + rows.append(rejected(lambda: root.lower(), "changed-root-" + name)) + setattr(root, name, old) + old = root.indices[:] + root.indices[0] = 1 + rows.append(rejected(lambda: root.lower(), "changed-root-wire-order")) + root.indices = old + old = backward.endpoint + backward.endpoint = 1 + rows.append(rejected(lambda: backward.execute(vector), "changed-native-endpoint-envelope")) + backward.endpoint = old + old = root.plan[0][0] + root.plan[0][0] = ((0,), E.Mat([[1]])) + rows.append(rejected(lambda: forward.execute(vector), "changed-stored-logical-plan")) + root.plan[0][0] = old + old = forward.stages[0] + forward.stages[0] = ("D", E.Vec([1])) + rows.append(rejected(lambda: forward.execute(vector), "changed-prepared-diagonal-page")) + forward.stages[0] = old + rows.append(rejected(lambda: DenseSources({DENSE_PREFIX + "full.py": source.raw[DENSE_PREFIX + "full.py"] + b"\n"}), + "changed-frozen-full-source")) + try: + source.raw[DENSE_PREFIX + "full.py"] = b"changed" + except TypeError: + rows.append("authenticated-byte-map-is-read-only") + else: + raise ValueError("authenticated source map accepted a replacement") + need(canonical({"n": 1}) != canonical({"n": True}), "typed dense receipt") + rows.append("typed-receipt") + return rows + + +def run(quick=False): + source = DenseSources() + E = source.e + fixtures = [("complex-scalar", E.Mat([[E.C(Q(1, 3), Q(1, 7))]]), False, False), + ("small-pivot", E.Mat([[Q(1, 1 << 80)]]), False, False), + ("large-scalar", E.Mat([[E.C(1 << 30, -(1 << 28))]]), False, False)] + for n in (2, 4, 8): + matrix = E.Mat([[E.C(1, Q(1, 16)) if i == j else + E.C(Q(-1 if (i - j) % 2 else 1, 32 * n), Q(i - j, 64 * n)) + for j in range(n)] for i in range(n)]) + fixtures.append(("complex-near-identity-" + str(n), matrix, n >= 4, n == 4)) + for n in (2, 4): + fixtures.append(("indefinite-" + str(n), E.Mat([[1 if i == j else 2 if abs(i - j) == 1 else 0 + for j in range(n)] for i in range(n)]), True, False)) + fixtures.append(("half-shift-Gaussian", E.Mat([[Q(1, 2) ** ((i - j) ** 2 - (i - j)) + for j in range(4)] for i in range(4)]), False, False)) + fixtures.append(("unequal-subtrees", E.Mat([[1 if i == j else Q(1, 16) if abs(i - j) == 4 + else Q(1, 32) if abs(i - j) == 5 else 0 + for j in range(8)] for i in range(8)]), False, False)) + fixtures.append(("zero-call-identity-16", identity(E, 16), True, False)) + if quick: + fixtures = [fixtures[0], fixtures[3]] + cases = [full_case(source, *fixture) for fixture in fixtures] + result = dict(schema_version=1, scope="actual conditioned full Toeplitz and stored inverse fixed-H execution", + sources=source.seals, cases=cases, generic_cases=generic_cases(source), + endpoint_factor_control=endpoint_control(source)) + if not quick: + result.update(word_cases=segment_cases(source), controls=controls(source), + actual_full_source_execution=True, native_stage_arithmetic=True, + physical_root_execution=False, physical_tape_execution=False, + complete_tensor_execution=False, dyadic_programs_exact_inverse_pair=False, + independent_transcendental_backend=False, universal_proof_kernel=False) + need(any(d["execution"]["nonzero_native_error"] for row in cases for d in row["directions"]), + "nonzero complete rounded error witness") + return result + + +def main(): + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--quick", action="store_true") + parser.add_argument("--write-result", action="store_true") + args = parser.parse_args() + need(not (args.quick and args.write_result), "diagnostic cannot replace complete receipt") + result = run(args.quick) + if args.quick: + print("PASS_NATIVE_FULL_DIAGNOSTIC", len(result["cases"])) + return + path = HERE / "DENSE-RESULT.json" + if args.write_result: + path.write_text(json.dumps(result, indent=2) + "\n") + else: + need(canonical(result) == canonical(json.loads(path.read_text())), "complete typed conditioned native receipt") + directions = [d for row in result["cases"] for d in row["directions"]] + print(json.dumps(dict(status="PASS_CONDITIONED_TOEPLITZ_FIXED_H", cases=len(result["cases"]), + directions=len(directions), generic_cases=len(result["generic_cases"]), + word_cases=len(result["word_cases"]), controls=len(result["controls"]), + native_stages=sum(d["counts"]["primitive_stages"] for d in directions), + rounded_outputs=sum(d["execution"]["statistics"]["complex_roundings"] for d in directions), + result_sha256=sha256(canonical(result)).hexdigest()), sort_keys=True)) + + +if __name__ == "__main__": + main()