#!/usr/bin/env python3 """Exhaustive, exact search in a bounded-degree Dirichlet product family. Only sorted tuples of positive factor lengths are visited. This removes permutations but imposes no additional condition on the product family. After one common denominator is chosen for each Fourier vector, every candidate comparison is an integer comparison. The search is not a search over all integer polynomials; see exact_toeplitz.py for that task. """ from __future__ import annotations import argparse import csv import hashlib import json import math import platform import time from fractions import Fraction from pathlib import Path from typing import Iterator from exact_toeplitz import dirichlet_product, normalized_fourier, normalized_objective F = Fraction def length_tuples(factors: int, budget: int, smallest: int = 1) -> Iterator[tuple[int, ...]]: if factors == 0: yield () return for first in range(smallest, budget // factors + 1): for rest in length_tuples(factors - 1, budget - first, first): yield (first,) + rest def sparse_binomial_product(lengths: tuple[int, ...]) -> dict[int, int]: result = {0: 1} for length in lengths: shifted = result.copy() for exponent, coefficient in result.items(): shifted[exponent + length] = shifted.get(exponent + length, 0) - coefficient result = {exponent: coefficient for exponent, coefficient in shifted.items() if coefficient} return result def integer_fourier(alpha: F, degree: int) -> tuple[int, list[int]]: fractions = normalized_fourier(alpha, degree) denominator = math.lcm(*(f.denominator for f in fractions)) return denominator, [f.numerator * (denominator // f.denominator) for f in fractions] def autocorrelation(polynomial: dict[int, int]) -> list[tuple[int, int]]: items = sorted(polynomial.items()) result = {0: sum(coefficient * coefficient for _, coefficient in items)} for i, (first, a) in enumerate(items): for second, b in items[i + 1:]: lag = second - first result[lag] = result.get(lag, 0) + 2 * a * b return [(lag, value) for lag, value in result.items() if value] def search_group(parameters: list[F], degree: int) -> list[dict]: if not parameters: return [] factors = int(parameters[0]) if any(lam <= 0 or lam.denominator == 1 or int(lam) != factors for lam in parameters): raise ValueError("Each group must have the same integer part and noninteger positive parameters") start = time.perf_counter() fourier = [integer_fourier(lam - factors, degree + factors) for lam in parameters] scores: list[int | None] = [None] * len(parameters) winners: list[list[tuple[int, ...]]] = [[] for _ in parameters] count = 0 for lengths in length_tuples(factors, degree + factors): correlation = autocorrelation(sparse_binomial_product(lengths)) for i, (_, coefficients) in enumerate(fourier): score = sum(value * coefficients[lag] for lag, value in correlation) if scores[i] is None or score < scores[i]: scores[i] = score winners[i] = [lengths] elif score == scores[i]: winners[i].append(lengths) count += 1 elapsed = time.perf_counter() - start results = [] for i, lam in enumerate(parameters): alpha = lam - factors ratio_alpha = F(scores[i], fourier[i][0]) lengths = winners[i][0] polynomial = dirichlet_product(list(lengths)) ratio_lambda = normalized_objective(polynomial, normalized_fourier(lam, degree)) # Gamma factors cancel under an integral shift, giving an independent # exact comparison between the Q and P coefficient calculations. c0_ratio = F(1) for j in range(factors): c0_ratio *= 2 * (2 * (alpha + j) + 1) / (alpha + j + 1) if ratio_alpha != c0_ratio * ratio_lambda: raise ArithmeticError("sparse-product and original-polynomial objectives disagree") c0_alpha_display = math.gamma(2 * float(alpha) + 1) / math.gamma(float(alpha) + 1) ** 2 results.append({ "lambda": str(lam), "degree_bound": degree, "factor_count": factors, "status": "complete_product_family", "tuples_visited": count, "arithmetic": "integer arithmetic after a common rational Fourier denominator", "lengths": list(lengths), "all_winning_length_tuples": [list(x) for x in winners[i]], "polynomial": polynomial, "polynomial_degree": len(polynomial) - 1, "ratio_to_c0_alpha": str(ratio_alpha), "ratio_to_c0_lambda": str(ratio_lambda), "integral_over_2pi_display": float(ratio_alpha) * c0_alpha_display, "elapsed_seconds_for_group": elapsed, "python_version": platform.python_version(), }) return results def main() -> None: parser = argparse.ArgumentParser(description=__doc__) parser.add_argument("--degree", type=int, required=True) parser.add_argument("--parameters", default="3/2,5/2,7/2,17/4,9/2,11/2", help="comma-separated positive noninteger rational parameters") parser.add_argument("--output", type=Path) parser.add_argument("--csv", type=Path) args = parser.parse_args() parameters = [F(value) for value in args.parameters.split(",")] groups: dict[int, list[F]] = {} for lam in parameters: groups.setdefault(int(lam), []).append(lam) rows = [] for group in groups.values(): rows.extend(search_group(group, args.degree)) rows.sort(key=lambda row: F(row["lambda"])) report = { "scope": "all Dirichlet products with k=floor(lambda) and degree <= N", "script_sha256": hashlib.sha256(Path(__file__).read_bytes()).hexdigest(), "results": rows, } serialized = json.dumps(report, indent=2, sort_keys=True) + "\n" if args.output: args.output.parent.mkdir(parents=True, exist_ok=True) args.output.write_text(serialized) if args.csv: args.csv.parent.mkdir(parents=True, exist_ok=True) columns = ["lambda", "degree_bound", "factor_count", "status", "tuples_visited", "lengths", "polynomial_degree", "ratio_to_c0_lambda", "integral_over_2pi_display", "elapsed_seconds_for_group"] with args.csv.open("w", newline="") as stream: writer = csv.DictWriter(stream, fieldnames=columns) writer.writeheader() for row in rows: writer.writerow({column: row[column] for column in columns}) print(serialized, end="") if __name__ == "__main__": main()