From b1895eb132fe9a67798c1878acec0df46b53ea11 Mon Sep 17 00:00:00 2001 From: Nick Petalas Date: Fri, 30 Oct 2026 13:53:04 +0000 Subject: [PATCH] feat(research): classify common CRT Fourier permutations --- .../crt-fourier-common-permutation/PROOF.md | 145 +++++++++++ .../crt-fourier-common-permutation/README.md | 122 +++++++++ .../SOURCES.json | 30 +++ .../crt-fourier-common-permutation/TASK.json | 90 +++++++ .../VERIFY.json | 37 +++ .../crt-fourier-common-permutation/verify.py | 238 ++++++++++++++++++ 6 files changed, 662 insertions(+) create mode 100644 research/crt-fourier-common-permutation/PROOF.md create mode 100644 research/crt-fourier-common-permutation/README.md create mode 100644 research/crt-fourier-common-permutation/SOURCES.json create mode 100644 research/crt-fourier-common-permutation/TASK.json create mode 100644 research/crt-fourier-common-permutation/VERIFY.json create mode 100644 research/crt-fourier-common-permutation/verify.py diff --git a/research/crt-fourier-common-permutation/PROOF.md b/research/crt-fourier-common-permutation/PROOF.md new file mode 100644 index 0000000000..f87d39e953 --- /dev/null +++ b/research/crt-fourier-common-permutation/PROOF.md @@ -0,0 +1,145 @@ +# One permutation on both Fourier indices + +This contribution was produced by Nick Petalas. Apache-2.0. + +## Contract and theorem + +Let `D >= 1`, and let `s_1,...,s_D >= 2` be pairwise-coprime integers. +Put `N = product_i s_i`, `M_i = N/s_i`, and +`G = product_i Z/(s_i)`. Representatives are canonical integers. +Use the standard positive-exponent, unnormalized Fourier kernels + +```text +K_N(a,b) = exp(2 pi i a b / N), +K_G(x,y) = exp(2 pi i sum_i x_i y_i / s_i). +``` + +The same statements hold with both exponent signs reversed, or with the +usual unitary normalization on both sides. Normalizing just one side does +not give the displayed identity. + +A bijection `f: G -> Z/(N)` satisfies + +```text +K_N(f(x), f(y)) = K_G(x,y) for every x,y in G (1) +``` + +if and only if it has the form + +```text +f(x) = sum_i M_i v_i x_i mod N, +M_i v_i^2 = 1 mod s_i for every i. (2) +``` + +Consequently, such an `f` exists exactly when every `M_i` is a square +unit modulo `s_i`. The number of such bijections is the product of the +numbers of solutions to the local congruences in (2). These are all the +solutions, even when the initial search allows arbitrary nonlinear +permutations. + +This is an exact interface classification. Local roots are fixed to the +standard roots in the displayed kernels. There are no extra diagonal +phases, changes of local primitive root, different permutations on the two +indices, or arbitrary linear input/output transforms. + +## Necessity, including arbitrary bijections + +Equality of the two roots of unity in (1) is equivalent to + +```text +f(x) f(y) = sum_i M_i x_i y_i mod N. (3) +``` + +Since `f` is onto, there is a `y_*` with `f(y_*) = 1`. Taking `x=0` in +(3) gives `f(0)=0`. For any `x,x',y`, subtract the three instances of (3) +at `(x+x',y)`, `(x,y)` and `(x',y)`. The right sides cancel modulo `N`, +including any coordinate wrap, because `M_i s_i=N`. Thus + +```text +(f(x+x') - f(x) - f(x')) f(y) = 0 mod N. +``` + +Take `y=y_*`. It follows that `f` is a group homomorphism. Nonlinear +bijections cannot evade this conclusion. + +Let `e_i` be the unit vector in factor `i`. Since `s_i e_i=0`, additivity +gives `s_i f(e_i)=0 mod N`. Hence `M_i` divides the canonical integer +`f(e_i)`, and `f(e_i)=M_i v_i` for a unique `0<=v_i None: + if not condition: + raise ValueError(reason) + + +def bounded_read(path: Path) -> bytes: + with path.open("rb") as stream: + data = stream.read(MAX_BYTES + 1) + require(len(data) <= MAX_BYTES, "file exceeds checker capacity") + return data + + +def source_hashes() -> dict[str, str]: + manifest = json.loads(bounded_read(PACKAGE / "SOURCES.json")) + rows = manifest["files"] + require(isinstance(rows, list) and len(rows) <= 16, "source inventory") + hashes = {} + for row in rows: + if row["usage"] != "input": + continue + name = row["path"] + require(name in { + "research/native-coprime-grid/README.md", + "research/native-coprime-grid/PROOF.md", + }, "unexpected source path") + data = bounded_read(PACKAGE.parents[1] / name) + digest = hashlib.sha256(data).hexdigest() + require(digest == row["sha256"], "source SHA256: " + name) + header = b"blob " + str(len(data)).encode("ascii") + b"\0" + require(hashlib.sha1(header + data).hexdigest() == row["git_blob"], + "source Git blob: " + name) + hashes[name] = digest + require(len(hashes) == 2, "missing source inputs") + return hashes + + +def shape_check(axes: tuple[int, ...]) -> int: + require(type(axes) is tuple and 1 <= len(axes) <= 4, "axis tuple") + require(all(type(s) is int and 2 <= s <= 257 for s in axes), + "axis bounds") + total = math.prod(axes) + require(total <= 65_536, "product bound") + for i in range(len(axes)): + for j in range(i): + require(math.gcd(axes[i], axes[j]) == 1, "axes not coprime") + return total + + +def root_coefficients(axes: tuple[int, ...]) -> tuple[tuple[int, ...], ...]: + """All f(e_i) tuples, in canonical residue order, from local square roots.""" + total = shape_check(axes) + choices = [] + for s in axes: + complement = total // s + roots = tuple(v for v in range(s) if complement * v * v % s == 1) + choices.append(tuple(complement * v for v in roots)) + return tuple(itertools.product(*choices)) + + +def coordinates(axes: tuple[int, ...]) -> tuple[tuple[int, ...], ...]: + shape_check(axes) + return tuple(itertools.product(*(range(s) for s in axes))) + + +def make_map(axes: tuple[int, ...], coefficients: tuple[int, ...]) -> tuple[int, ...]: + total = shape_check(axes) + require(type(coefficients) is tuple and len(coefficients) == len(axes), + "coefficient tuple") + require(all(type(c) is int and 0 <= c < total for c in coefficients), + "coefficient bounds") + return tuple(sum(c * x for c, x in zip(coefficients, point, strict=True)) + % total for point in coordinates(axes)) + + +def kernel_holds(axes: tuple[int, ...], mapping: tuple[int, ...]) -> bool: + """Literal equality of every phase, with no floating complex arithmetic.""" + total = shape_check(axes) + if (type(mapping) is not tuple or len(mapping) != total + or any(type(x) is not int for x in mapping) + or sorted(mapping) != list(range(total))): + return False + points = coordinates(axes) + weights = tuple(total // s for s in axes) + for ix in range(total): + for iy in range(total): + tensor = sum(w * x * y for w, x, y in + zip(weights, points[ix], points[iy], strict=True)) + if (mapping[ix] * mapping[iy] - tensor) % total: + return False + return True + + +def cyclic_oracle(axes: tuple[int, ...]) -> tuple[set[tuple[int, ...]], int]: + """Enumerate every homomorphism by the image of (1,...,1), not local roots.""" + total = shape_check(axes) + points = coordinates(axes) + # Build the cyclic-coordinate table by direct simultaneous residues. + lookup = {tuple(k % s for s in axes): k for k in range(total)} + require(len(lookup) == total, "cyclic coordinate bijection") + accepted = set() + trials = 0 + for multiplier in range(total): + mapping = tuple(multiplier * lookup[x] % total for x in points) + trials += 1 + if kernel_holds(axes, mapping): + accepted.add(mapping) + return accepted, trials + + +def finite_shapes() -> tuple[tuple[int, ...], ...]: + shapes = [(s,) for s in range(2, 33)] + for count in (2, 3): + for axes in itertools.combinations(range(2, 14), count): + if math.prod(axes) <= 255 and all( + math.gcd(a, b) == 1 for a, b in itertools.combinations(axes, 2)): + shapes.append(axes) + return tuple(shapes) + + +def compare_maps() -> dict[str, int]: + counts = {"shapes": 0, "cyclic_candidates": 0, "accepted_maps": 0, + "literal_phase_entries": 0, "no_map_shapes": 0} + for axes in finite_shapes(): + produced = {make_map(axes, c) for c in root_coefficients(axes)} + expected, trials = cyclic_oracle(axes) + require(produced == expected, "root constructor versus cyclic oracle") + total = math.prod(axes) + for mapping in produced: + require(kernel_holds(axes, mapping), "constructed Fourier identity") + counts["literal_phase_entries"] += total * total + counts["shapes"] += 1 + counts["cyclic_candidates"] += trials + counts["accepted_maps"] += len(produced) + counts["no_map_shapes"] += not produced + return counts + + +def arbitrary_permutations() -> dict[str, int]: + counts = {"full_permutations": 0, "accepted_permutations": 0, + "nonlinear_accepted": 0} + for axes in ((2,), (3,), (4,), (5,), (6,), (2, 3)): + total = math.prod(axes) + expected = {make_map(axes, c) for c in root_coefficients(axes)} + accepted = set() + # Include every translation and every nonlinear bijection, not just f(0)=0. + for mapping in itertools.permutations(range(total)): + counts["full_permutations"] += 1 + if kernel_holds(axes, mapping): + accepted.add(mapping) + counts["accepted_permutations"] += 1 + counts["nonlinear_accepted"] += mapping not in expected + require(accepted == expected, "arbitrary permutation classification") + require(counts["nonlinear_accepted"] == 0, "nonlinear escape") + return counts + + +def examples() -> dict[str, object]: + positive = (13, 17) + coefficients = root_coefficients(positive) + require(len(coefficients) == 4, "13 by 17 count") + for c in coefficients: + require(kernel_holds(positive, make_map(positive, c)), "positive fixture") + require(not root_coefficients((3, 5)), "3 by 5 obstruction") + require(not root_coefficients((3, 7)), "3 by 7 obstruction") + # Standard two-map CRT still works at these excluded common-map shapes. + two_map_entries = 0 + for axes in ((3, 5), (3, 7), (4, 5)): + total = math.prod(axes) + points = coordinates(axes) + left = tuple((total // s) * pow(total // s, -1, s) for s in axes) + right = tuple(total // s for s in axes) + for x in points: + for y in points: + a = sum(c * z for c, z in zip(left, x, strict=True)) % total + b = sum(c * z for c, z in zip(right, y, strict=True)) % total + phase = sum(w * u * v for w, u, v in zip(right, x, y, strict=True)) + require((a * b - phase) % total == 0, "two-map scope control") + two_map_entries += 1 + return {"axes": positive, "coefficients": coefficients, + "two_map_scope_entries": two_map_entries} + + +def rejection_controls() -> int: + invalid = [(), (1,), (True,), (2.0,), (258,), (2, 2), (6, 9), + (2, 3, 5, 7, 11), (251, 257, 3)] + for axes in invalid: + try: + root_coefficients(axes) + except ValueError: + pass + else: + raise ValueError("invalid shape accepted") + axes = (13, 17) + good = make_map(axes, root_coefficients(axes)[0]) + shifted = tuple((x + 1) % 221 for x in good) + nonlinear = list(good) + nonlinear[1], nonlinear[2] = nonlinear[2], nonlinear[1] + require(not kernel_holds(axes, shifted), "translated common map") + require(not kernel_holds(axes, tuple(nonlinear)), "nonlinear swap") + require(not kernel_holds(axes, good[:-1]), "truncated map") + require(not kernel_holds(axes, (0,) * 221), "nonbijective map") + naive = make_map((3, 5), (10, 6)) + require(not kernel_holds((3, 5), naive), "same canonical CRT map") + return len(invalid) + 5 + + +def main() -> None: + inputs = source_hashes() + report = {"status": "PASS_CRT_FOURIER_COMMON_PERMUTATION", + "classification": compare_maps(), + "arbitrary_bijections": arbitrary_permutations(), + "examples": examples(), "rejection_controls": rejection_controls(), + "source_sha256": inputs} + names = ("README.md", "PROOF.md", "TASK.json", "SOURCES.json", "VERIFY.json", "verify.py") + report["package_sha256"] = { + name: hashlib.sha256(bounded_read(PACKAGE / name)).hexdigest() for name in names + } + print(json.dumps(report, indent=2, sort_keys=True)) + + +if __name__ == "__main__": + main()