diff --git a/.github/workflows/dirty-stride-fast-phases.yml b/.github/workflows/dirty-stride-fast-phases.yml new file mode 100644 index 0000000000..473df7e3e1 --- /dev/null +++ b/.github/workflows/dirty-stride-fast-phases.yml @@ -0,0 +1,27 @@ +name: Dirty stride fast phase stream +on: + pull_request: + paths: + - 'research/dirty-stride-fast-phases/**' + - '.github/workflows/dirty-stride-fast-phases.yml' + push: + paths: + - 'research/dirty-stride-fast-phases/**' + - '.github/workflows/dirty-stride-fast-phases.yml' + workflow_dispatch: +permissions: + contents: read +jobs: + verify: + runs-on: ubuntu-latest + timeout-minutes: 20 + strategy: + matrix: + python: ['3.11', '3.14'] + steps: + - uses: actions/checkout@v4 + - uses: actions/setup-python@v5 + with: + python-version: ${{ matrix.python }} + - run: python -I -O -B research/dirty-stride-fast-phases/verify.py + - run: python -I -O -B research/dirty-stride-fast-phases/test_package.py diff --git a/research/dirty-stride-fast-phases/BENCHMARK.json b/research/dirty-stride-fast-phases/BENCHMARK.json new file mode 100644 index 0000000000..dca383548b --- /dev/null +++ b/research/dirty-stride-fast-phases/BENCHMARK.json @@ -0,0 +1,265 @@ +{ + "platform": "macOS-26.6.1-arm64-arm-64bit-Mach-O", + "python": "3.14.7", + "rows": [ + { + "bits": 16, + "length": 32, + "median_seconds": { + "pr1038": 0.042217166163027287, + "stream": 0.1914247078821063 + }, + "trials": { + "pr1038": [ + { + "erase_seconds": 3.3294782042503357e-07, + "phase_divisions": 7056, + "phase_products": 14118, + "prepare_seconds": 0.0006293328478932381, + "seed_divisions": 0, + "seed_peak_bits": null, + "seed_products": 0, + "stream_seconds": 0.04158750036731362, + "tape_moves": 221, + "total_seconds": 0.042217166163027287 + }, + { + "erase_seconds": 2.0815059542655945e-07, + "phase_divisions": 7056, + "phase_products": 14118, + "prepare_seconds": 0.0006140000186860561, + "seed_divisions": 0, + "seed_peak_bits": null, + "seed_products": 0, + "stream_seconds": 0.04078587470576167, + "tape_moves": 221, + "total_seconds": 0.041400082875043154 + }, + { + "erase_seconds": 3.7532299757003784e-07, + "phase_divisions": 7056, + "phase_products": 14118, + "prepare_seconds": 0.0006021247245371342, + "seed_divisions": 0, + "seed_peak_bits": null, + "seed_products": 0, + "stream_seconds": 0.04313837504014373, + "tape_moves": 221, + "total_seconds": 0.04374087508767843 + } + ], + "stream": [ + { + "erase_seconds": 0.08656962495297194, + "phase_divisions": 384, + "phase_products": 966, + "prepare_seconds": 0.08834325009956956, + "seed_divisions": 5178, + "seed_peak_bits": 11000, + "seed_products": 5184, + "stream_seconds": 0.01533125014975667, + "tape_moves": 221, + "total_seconds": 0.19024412520229816 + }, + { + "erase_seconds": 0.09268087521195412, + "phase_divisions": 384, + "phase_products": 966, + "prepare_seconds": 0.09867466613650322, + "seed_divisions": 5178, + "seed_peak_bits": 11000, + "seed_products": 5184, + "stream_seconds": 0.01634766673669219, + "tape_moves": 221, + "total_seconds": 0.20770320808514953 + }, + { + "erase_seconds": 0.08787808287888765, + "phase_divisions": 384, + "phase_products": 966, + "prepare_seconds": 0.08811087487265468, + "seed_divisions": 5178, + "seed_peak_bits": 11000, + "seed_products": 5184, + "stream_seconds": 0.015435750130563974, + "tape_moves": 221, + "total_seconds": 0.1914247078821063 + } + ] + } + }, + { + "bits": 48, + "length": 257, + "median_seconds": { + "pr1038": 0.6899422500282526, + "stream": 1.2866912498138845 + }, + "trials": { + "pr1038": [ + { + "erase_seconds": 2.9103830456733704e-07, + "phase_divisions": 92340, + "phase_products": 184686, + "prepare_seconds": 0.001130416989326477, + "seed_divisions": 0, + "seed_peak_bits": null, + "seed_products": 0, + "stream_seconds": 0.5955920000560582, + "tape_moves": 1796, + "total_seconds": 0.5967227080836892 + }, + { + "erase_seconds": 3.334134817123413e-07, + "phase_divisions": 92340, + "phase_products": 184686, + "prepare_seconds": 0.001091458834707737, + "seed_divisions": 0, + "seed_peak_bits": null, + "seed_products": 0, + "stream_seconds": 1.1057697078213096, + "tape_moves": 1796, + "total_seconds": 1.106861500069499 + }, + { + "erase_seconds": 2.0815059542655945e-07, + "phase_divisions": 92340, + "phase_products": 184686, + "prepare_seconds": 0.001532666850835085, + "seed_divisions": 0, + "seed_peak_bits": null, + "seed_products": 0, + "stream_seconds": 0.6884093750268221, + "tape_moves": 1796, + "total_seconds": 0.6899422500282526 + } + ], + "stream": [ + { + "erase_seconds": 0.59731016587466, + "phase_divisions": 3084, + "phase_products": 7716, + "prepare_seconds": 0.6183066670782864, + "seed_divisions": 7914, + "seed_peak_bits": 25350, + "seed_products": 7920, + "stream_seconds": 0.15272741671651602, + "tape_moves": 1796, + "total_seconds": 1.3683442496694624 + }, + { + "erase_seconds": 0.5393960415385664, + "phase_divisions": 3084, + "phase_products": 7716, + "prepare_seconds": 0.5949774170294404, + "seed_divisions": 7914, + "seed_peak_bits": 25350, + "seed_products": 7920, + "stream_seconds": 0.15231779124587774, + "tape_moves": 1796, + "total_seconds": 1.2866912498138845 + }, + { + "erase_seconds": 0.46511820796877146, + "phase_divisions": 3084, + "phase_products": 7716, + "prepare_seconds": 0.5326998331584036, + "seed_divisions": 7914, + "seed_peak_bits": 25350, + "seed_products": 7920, + "stream_seconds": 0.18261054204776883, + "tape_moves": 1796, + "total_seconds": 1.180428583174944 + } + ] + } + }, + { + "bits": 96, + "length": 2048, + "median_seconds": { + "pr1038": 10.187936124857515, + "stream": 9.092779374681413 + }, + "trials": { + "pr1038": [ + { + "erase_seconds": 4.7497451305389404e-06, + "phase_divisions": 1146600, + "phase_products": 2293206, + "prepare_seconds": 0.0018824171274900436, + "seed_divisions": 0, + "seed_peak_bits": null, + "seed_products": 0, + "stream_seconds": 10.785656167194247, + "tape_moves": 14333, + "total_seconds": 10.787543334066868 + }, + { + "erase_seconds": 1.4170072972774506e-06, + "phase_divisions": 1146600, + "phase_products": 2293206, + "prepare_seconds": 0.0027474593371152878, + "seed_divisions": 0, + "seed_peak_bits": null, + "seed_products": 0, + "stream_seconds": 9.294908415991813, + "tape_moves": 14333, + "total_seconds": 9.297657292336226 + }, + { + "erase_seconds": 2.250075340270996e-06, + "phase_divisions": 1146600, + "phase_products": 2293206, + "prepare_seconds": 0.0022939168848097324, + "seed_divisions": 0, + "seed_peak_bits": null, + "seed_products": 0, + "stream_seconds": 10.185639957897365, + "tape_moves": 14333, + "total_seconds": 10.187936124857515 + } + ], + "stream": [ + { + "erase_seconds": 2.6205718340352178, + "phase_divisions": 24576, + "phase_products": 61446, + "prepare_seconds": 3.9281516657210886, + "seed_divisions": 11802, + "seed_peak_bits": 55767, + "seed_products": 11808, + "stream_seconds": 2.5440558749251068, + "tape_moves": 14333, + "total_seconds": 9.092779374681413 + }, + { + "erase_seconds": 3.426748082973063, + "phase_divisions": 24576, + "phase_products": 61446, + "prepare_seconds": 3.2136184577830136, + "seed_divisions": 11802, + "seed_peak_bits": 55767, + "seed_products": 11808, + "stream_seconds": 2.2742746672593057, + "tape_moves": 14333, + "total_seconds": 8.914641208015382 + }, + { + "erase_seconds": 3.318305708002299, + "phase_divisions": 24576, + "phase_products": 61446, + "prepare_seconds": 4.192794916685671, + "seed_divisions": 11802, + "seed_peak_bits": 55767, + "seed_products": 11808, + "stream_seconds": 1.683004958089441, + "tape_moves": 14333, + "total_seconds": 9.19410558277741 + } + ] + } + } + ], + "scope": "Two opposite phase passes; new seed erasure included; old API has no explicit seed erasure. PR1038 puts all product and division counts in the phase counters. The new seed counters are separate." +} diff --git a/research/dirty-stride-fast-phases/CHECKS.md b/research/dirty-stride-fast-phases/CHECKS.md new file mode 100644 index 0000000000..e9a6a86aa7 --- /dev/null +++ b/research/dirty-stride-fast-phases/CHECKS.md @@ -0,0 +1,30 @@ +# Verification contracts + +The test authoring gate was applied before the checks were added. + +| Contract | Credible failure | Check and distinct coverage | +| --- | --- | --- | +| Uniform state and output error | Wrong second difference, chart, or guard width | Direct quadratic residues and full-angle outward rational intervals | +| Exact stream and seed cleanup | Wrong inverse order or omitted arithmetic uncomputation | Both phase-state registers and seed outputs reach exact zero | +| Dirty operator and inverse | Wrong outer phase direction or target-dependent roundoff omitted | Actual leaf and factor calls with odd lowest bits, large dirt, and geometric-sum oracle | +| Complete access bill | An input read, data move, or inverse pass is not charged | Data banks count literal accesses on two physical tapes and reject copies or slices | +| Linear stream work | Term evaluation remains inside each phase site | Executed product, division, and seed-term counts at three lengths and precisions | +| Reproducible package | Missing local dependency, altered source, or unchecked receipt | Source pins, semantic mutations, and relocated full receipt replay | + +These checks extend the phase-interface boundary. They do not duplicate +the full selected-frontier replay or claim the inherited dirty solver has +been re-proved. No production export or injection hook was added for tests. +The new public phase lifecycle is also used by each production data pass. + +The small phase domain is all residues for moduli 1 through 9 at `W=2` +and `W=24`, with `L=2m+1`. Seven extra cases cover wider precision, +negative or wide numerators, large moduli, and chart boundaries. One long +case uses `a=12345`, `m=65537`, `W=24`, and `L=16384`; all steps execute, +and indices divisible by 257 plus the last index get interval checks. +There is no random seed or open-ended search. Every stated case completes. + +The dirty-word domain has lengths 1, 5, 8, and 17, strides 0, 2, -3, +and 7, both signs, and one or three factors. Each call is followed by +its exact inverse. One-factor cases use no borrowed array; three-factor +cases retain two arbitrary borrowed arrays. This is finite coverage of +the written all-length theorem, not an exhaustive input audit. diff --git a/research/dirty-stride-fast-phases/CLAIM.md b/research/dirty-stride-fast-phases/CLAIM.md new file mode 100644 index 0000000000..a6d60a21e4 --- /dev/null +++ b/research/dirty-stride-fast-phases/CLAIM.md @@ -0,0 +1,30 @@ +# Local route 8 claim + +Base: `abdc6a7002abdbcd9e21bfc5e1a795784deacfdf`. +Branch: `codex/k9-fast-phases`. Start: 2026-10-13 18:04 UTC. + +No pinned sprint issue was supplied. The coordinator records the claim locally. +The current PR metadata lists PR1038 at `99c71d385171480ecc2c43287aa4a04e1027fcd9` +and PR1054 at `6bc675be1554b3a3be771147f21f9dc2c6dae50a`. + +Candidate 1: use two reversible rounded complex phase states, a fixed second +difference, and exact backward reconstruction. Pay logarithmic guard precision. +Candidate 2: use exact modular angle state and blocks of polynomial values. +This needs a paid polynomial table or repeated evaluation at block boundaries. +Candidate 3: prove a work/error tradeoff for noninvertible rounded products. +This does not give the requested reversible implementation. + +Chosen object: the first candidate, for all integer numerators, positive moduli, +finite positive lengths, and working precision at least two bits. Prove a uniform +quadratic drift bound. Use bounded three-shear charts. Prepare two seeds with +exact rational series, then clear their arithmetic state. Reverse the phase +state after every data pass. Keep the inherited borrowed-array count fixed. + +Expected result: an executable scalar phase interface with O(L) stream work, +polynomial seed work, full dirty-target error, and a separate bit/tape bill. +The complete kappa calculation is still complex-supplier bound. This result +can remove a phase-interface expense; it does not change a supplier invoice. +No new kappa is claimed. Expected finish: this research session, after checks. + +No random search is used. PR479, PR948, PR985, PR1038 and PR1054 are +predecessor or scope references. PR1054 is not combined by assumption. diff --git a/research/dirty-stride-fast-phases/DERIVATION.diff b/research/dirty-stride-fast-phases/DERIVATION.diff new file mode 100644 index 0000000000..19eeadc3fc --- /dev/null +++ b/research/dirty-stride-fast-phases/DERIVATION.diff @@ -0,0 +1,182 @@ +--- PR1038/chirp.py ++++ chirp.py +@@ -1,2 +1 @@ +-"""Use reversible quadratic phases around the pinned dirty-array solver.""" +-from dataclasses import dataclass ++"""Use the phase stream around the pinned dirty-array word.""" +@@ -9,159 +8,6 @@ +- 'dirty_phase_predecessor', Path(__file__).resolve().parent/'dirty_stride.py') +-d = importlib.util.module_from_spec(_SPEC) +-sys.modules[_SPEC.name] = d +-_SPEC.loader.exec_module(d) +- +- +-def real_round(value, bits): +- scale = 1 << bits +- value = value*scale + F(1, 2) +- return F(value.numerator//value.denominator, scale) +- +- +-@dataclass +-class Meter(d.Meter): +- phase_evaluations: int = 0 +- phase_passes: int = 0 +- phase_terms: int = 0 +- phase_products: int = 0 +- phase_divisions: int = 0 +- phase_record_moves: int = 0 +- phase_recurrence_steps: int = 0 +- input_reads: int = 0 +- descriptor_bits: int = 0 +- coefficient_bits: int = 0 +- internal_bits: int = 0 +- error_bound: F = F(0) +- initial_bound: F = F(0) +- phase_peak_l1: F = F(0) +- +- +-def pi_approx(bits, meter): +- """Use the Machin identity and round each arctangent term.""" +- value = F(0) +- for denominator, scale in ((5, 16), (239, -4)): +- power = denominator +- for k in range(bits): +- term = real_round(F(1, (2*k+1)*power), bits) +- value += scale*term if k % 2 == 0 else -scale*term +- power *= denominator*denominator +- meter.phase_terms += 1 +- meter.phase_divisions += 1 +- meter.phase_products += 2 +- return value +- +- +-@dataclass(frozen=True) +-class Phase: +- numerator: int +- modulus: int +- bits: int +- internal: int +- pi: F +- meter: Meter +- +- @classmethod +- def prepare(cls, numerator, modulus, bits, meter): +- d.need(type(numerator) is int and type(modulus) is int and modulus > 0, +- 'integer phase descriptor') +- d.need(type(bits) is int and bits >= 2, 'phase precision') +- coefficient_bits = bits+8 +- internal = coefficient_bits + 2*(coefficient_bits+32).bit_length()+20 +- d.need((8*internal+17)*F(1, 1 << internal) +- <= F(1, 1 << (coefficient_bits+2)), 'phase coefficient budget') +- meter.descriptor_bits += abs(numerator).bit_length()+modulus.bit_length()+2 +- meter.coefficient_bits = coefficient_bits +- meter.internal_bits = internal +- return cls(numerator % modulus, modulus, bits, internal, +- pi_approx(internal, meter), meter) +- +- def coefficients(self, residue): +- """Return a quarter-turn and two dyadic shear coefficients.""" +- d.need(type(residue) is int and 0 <= residue < self.modulus, 'phase residue') +- self.meter.phase_evaluations += 1 +- quadrant = (8*residue+self.modulus)//(2*self.modulus) +- remainder = 4*residue-quadrant*self.modulus +- if remainder == 0: +- return quadrant % 4, F(0), F(0) +- angle = real_round(self.pi*remainder/(2*self.modulus), self.internal) +- self.meter.phase_products += 1 +- self.meter.phase_divisions += 1 +- d.need(abs(angle) <= 1, 'reduced angle') +- sine = sine_term = angle +- cosine = cosine_term = F(1) +- for k in range(1, self.internal): +- sine_term = real_round(-sine_term*angle*angle/((2*k)*(2*k+1)), self.internal) +- cosine_term = real_round(-cosine_term*angle*angle/((2*k-1)*(2*k)), self.internal) +- sine += sine_term +- cosine += cosine_term +- self.meter.phase_terms += 2 +- self.meter.phase_products += 4 +- self.meter.phase_divisions += 2 +- tangent = real_round(sine/(1+cosine), self.internal) +- self.meter.phase_divisions += 1 +- tangent = real_round(tangent, self.bits+8) +- sine = real_round(sine, self.bits+8) +- d.need(abs(tangent) <= F(1, 2) and abs(sine) <= 1, 'bounded phase shears') +- return quadrant % 4, tangent, sine +- +- +-def quarter_turn(value, quadrant): +- real, imag = value +- if quadrant == 0: +- return value +- if quadrant == 1: +- return -imag, real +- if quadrant == 2: +- return -real, -imag +- return imag, -real +- +- +-def observe_phase(value, meter): +- if meter is not None: +- meter.phase_peak_l1 = max(meter.phase_peak_l1, abs(value[0])+abs(value[1])) +- +- +-def rotate(value, coefficients, bits, inverse=False, meter=None): +- """Apply three grid shears. The inverse uses the reverse instruction order.""" +- quadrant, tangent, sine = coefficients +- observe_phase(value, meter) +- if inverse: +- real, imag = quarter_turn(value, (-quadrant) % 4) +- real += real_round(tangent*imag, bits) +- observe_phase((real, imag), meter) +- imag -= real_round(sine*real, bits) +- observe_phase((real, imag), meter) +- real += real_round(tangent*imag, bits) +- observe_phase((real, imag), meter) +- result = real, imag +- else: +- real, imag = value +- real -= real_round(tangent*imag, bits) +- observe_phase((real, imag), meter) +- imag += real_round(sine*real, bits) +- observe_phase((real, imag), meter) +- real -= real_round(tangent*imag, bits) +- observe_phase((real, imag), meter) +- result = quarter_turn((real, imag), quadrant) +- if meter is not None: +- meter.phase_products += 3 +- return result +- +- +-def phase_pass(data, phase, inverse): +- """Scan one array. Use quadratic finite differences and no phase table.""" +- meter = phase.meter +- before = meter.tape_moves +- residue = 0 +- difference = phase.numerator +- for i in range(len(data)): +- coefficients = phase.coefficients(residue) +- meter.access(data, i) +- value = rotate(data[i], coefficients, phase.bits, inverse, meter) +- meter.access(data, i, True) +- data[i] = value +- residue = (residue+difference) % phase.modulus +- difference = (difference+2*phase.numerator) % phase.modulus +- meter.phase_recurrence_steps += 1 +- meter.phase_passes += 1 +- meter.phase_record_moves += meter.tape_moves-before ++ 'dirty_fast_phase_stream', Path(__file__).resolve().parent/'phase_stream.py') ++stream = importlib.util.module_from_spec(_SPEC) ++sys.modules[_SPEC.name] = stream ++_SPEC.loader.exec_module(stream) ++d = stream.d ++Meter, Phase, phase_pass = stream.Meter, stream.Phase, stream.phase_pass +@@ -223 +69 @@ +- phase = Phase.prepare(numerator, modulus, bits, meter) ++ phase = Phase.prepare(numerator, modulus, bits, meter, len(x)) +@@ -230,0 +77 @@ ++ phase.clear() +@@ -262 +109 @@ +- phase = Phase.prepare(numerator, modulus, bits, meter) ++ phase = Phase.prepare(numerator, modulus, bits, meter, len(x)) +@@ -269,0 +117 @@ ++ phase.clear() diff --git a/research/dirty-stride-fast-phases/EXPECTED.json b/research/dirty-stride-fast-phases/EXPECTED.json new file mode 100644 index 0000000000..8a86f29c1d --- /dev/null +++ b/research/dirty-stride-fast-phases/EXPECTED.json @@ -0,0 +1,1480 @@ +{ + "controls": [ + "zero_modulus", + "boolean_numerator", + "zero_length", + "boolean_precision", + "second_stream", + "early_seed_clear", + "cleared_seed" + ], + "costs": [ + { + "bits": 16, + "length": 32, + "new_seed_divisions": 5178, + "new_seed_peak_bits": 11000, + "new_seed_products": 5184, + "new_seed_terms": 3432, + "new_state_bits": 48, + "new_stream_divisions": 384, + "new_stream_products": 966, + "old_divisions": 7056, + "old_products": 14118, + "old_terms": 6932, + "tape_moves": 221 + }, + { + "bits": 48, + "length": 257, + "new_seed_divisions": 7914, + "new_seed_peak_bits": 25350, + "new_seed_products": 7920, + "new_seed_terms": 5256, + "new_state_bits": 86, + "new_stream_divisions": 3084, + "new_stream_products": 7716, + "old_divisions": 92340, + "old_products": 184686, + "old_terms": 91316, + "tape_moves": 1796 + }, + { + "bits": 96, + "length": 2048, + "new_seed_divisions": 11802, + "new_seed_peak_bits": 55767, + "new_seed_products": 11808, + "new_seed_terms": 7848, + "new_state_bits": 140, + "new_stream_divisions": 24576, + "new_stream_products": 61446, + "old_divisions": 1146600, + "old_products": 2293206, + "old_terms": 1138412, + "tape_moves": 14333 + } + ], + "phases": [ + { + "bits": 2, + "length": 3, + "max_squared_drift_ratio_log2": null, + "modulus": "1", + "numerator": "0", + "output_digest": "9035ca0f9a502deb17afde3277ce35aee6fc29c57fecc63e905a4abcc5c09fea", + "sample_stride": 1, + "samples": 3, + "seed_peak_bits": 0, + "state_bits": 28 + }, + { + "bits": 2, + "length": 5, + "max_squared_drift_ratio_log2": null, + "modulus": "2", + "numerator": "0", + "output_digest": "3265f76cc76da093be157697d85d1f7c2b56ca2fdc6ee05d051a63f873702854", + "sample_stride": 1, + "samples": 5, + "seed_peak_bits": 0, + "state_bits": 28 + }, + { + "bits": 2, + "length": 5, + "max_squared_drift_ratio_log2": -201, + "modulus": "2", + "numerator": "1", + "output_digest": "5c2c8fce03a50c47f5544df817965f84eb81942c5fba6cdf7c658513df63b64a", + "sample_stride": 1, + "samples": 5, + "seed_peak_bits": 0, + "state_bits": 28 + }, + { + "bits": 2, + "length": 7, + "max_squared_drift_ratio_log2": null, + "modulus": "3", + "numerator": "0", + "output_digest": "feee578192fae459090ae5adc46867ebce0742175247bd8fb3cfece8e12eae79", + "sample_stride": 1, + "samples": 7, + "seed_peak_bits": 0, + "state_bits": 30 + }, + { + "bits": 2, + "length": 7, + "max_squared_drift_ratio_log2": -23, + "modulus": "3", + "numerator": "1", + "output_digest": "bd117bbc3a0417640f85faa78005c64ecb274c3ebabaded5c59a591d4316ebd2", + "sample_stride": 1, + "samples": 7, + "seed_peak_bits": 6350, + "state_bits": 30 + }, + { + "bits": 2, + "length": 7, + "max_squared_drift_ratio_log2": -23, + "modulus": "3", + "numerator": "2", + "output_digest": "6b51645aa4ee8adc1a99b33151e91cf8303f751a7616ebdc1f6c79650de24e98", + "sample_stride": 1, + "samples": 7, + "seed_peak_bits": 6350, + "state_bits": 30 + }, + { + "bits": 2, + "length": 9, + "max_squared_drift_ratio_log2": null, + "modulus": "4", + "numerator": "0", + "output_digest": "feb7a39de8ed2a289234a545c3d716b4f1ef3241257a881da9a89947d91dd0a9", + "sample_stride": 1, + "samples": 9, + "seed_peak_bits": 0, + "state_bits": 30 + }, + { + "bits": 2, + "length": 9, + "max_squared_drift_ratio_log2": -202, + "modulus": "4", + "numerator": "1", + "output_digest": "3426cc72e8962ab4bb1a7e94d7bb58475b0efcb67add25f3b011862f863a21ef", + "sample_stride": 1, + "samples": 9, + "seed_peak_bits": 0, + "state_bits": 30 + }, + { + "bits": 2, + "length": 9, + "max_squared_drift_ratio_log2": -201, + "modulus": "4", + "numerator": "2", + "output_digest": "82f70db2051b0b8dde6ebd7f55164b45672b91203a8afd84863c51ffd766739d", + "sample_stride": 1, + "samples": 9, + "seed_peak_bits": 0, + "state_bits": 30 + }, + { + "bits": 2, + "length": 9, + "max_squared_drift_ratio_log2": -198, + "modulus": "4", + "numerator": "3", + "output_digest": "e2189eefb7efede961db1dc583465886021e73857e1b2fecf95a6237882834ce", + "sample_stride": 1, + "samples": 9, + "seed_peak_bits": 0, + "state_bits": 30 + }, + { + "bits": 2, + "length": 11, + "max_squared_drift_ratio_log2": null, + "modulus": "5", + "numerator": "0", + "output_digest": "1c0c33d09f730bcdd2a71a7daddc3d6435299115bd78a3f5aee9e372207a845f", + "sample_stride": 1, + "samples": 11, + "seed_peak_bits": 0, + "state_bits": 30 + }, + { + "bits": 2, + "length": 11, + "max_squared_drift_ratio_log2": -20, + "modulus": "5", + "numerator": "1", + "output_digest": "5fa9e57929211257ef1c3afbfdc1c0fbb2aeb524ecd1cb7d817f23f2de2fea25", + "sample_stride": 1, + "samples": 11, + "seed_peak_bits": 5767, + "state_bits": 30 + }, + { + "bits": 2, + "length": 11, + "max_squared_drift_ratio_log2": -24, + "modulus": "5", + "numerator": "2", + "output_digest": "bec24431db361dd4e9179558867f804f51ae4b774a34eecb7e6b8a748241b066", + "sample_stride": 1, + "samples": 11, + "seed_peak_bits": 5767, + "state_bits": 30 + }, + { + "bits": 2, + "length": 11, + "max_squared_drift_ratio_log2": -24, + "modulus": "5", + "numerator": "3", + "output_digest": "c7cf46acde165c2f660b43d21ea08b5d50ea9c6f832c2eb9f2f0eb704ba2766a", + "sample_stride": 1, + "samples": 11, + "seed_peak_bits": 5767, + "state_bits": 30 + }, + { + "bits": 2, + "length": 11, + "max_squared_drift_ratio_log2": -20, + "modulus": "5", + "numerator": "4", + "output_digest": "e26a6fa5f55a3ed40904b4a1e3d2699ffef0d456f3df8a05ada4fc97b60ceaff", + "sample_stride": 1, + "samples": 11, + "seed_peak_bits": 5767, + "state_bits": 30 + }, + { + "bits": 2, + "length": 13, + "max_squared_drift_ratio_log2": null, + "modulus": "6", + "numerator": "0", + "output_digest": "3af2a9879c4d85df93038c16b2808c771aa8d64af544f0f1b072d34ba7f1df77", + "sample_stride": 1, + "samples": 13, + "seed_peak_bits": 0, + "state_bits": 30 + }, + { + "bits": 2, + "length": 13, + "max_squared_drift_ratio_log2": -23, + "modulus": "6", + "numerator": "1", + "output_digest": "bd2944d72c84d0072eeab7c9a18373627bc81ed3cd9d8c3066d6a62e589d8fff", + "sample_stride": 1, + "samples": 13, + "seed_peak_bits": 6350, + "state_bits": 30 + }, + { + "bits": 2, + "length": 13, + "max_squared_drift_ratio_log2": -23, + "modulus": "6", + "numerator": "2", + "output_digest": "167253d82487d7eb32f76a70df649568b5bd901486212121c9b9a0e838089ecf", + "sample_stride": 1, + "samples": 13, + "seed_peak_bits": 6350, + "state_bits": 30 + }, + { + "bits": 2, + "length": 13, + "max_squared_drift_ratio_log2": -201, + "modulus": "6", + "numerator": "3", + "output_digest": "498365cb59f1e30e8259af12c0077abdc022477c4237b9a3b071b94135e795f5", + "sample_stride": 1, + "samples": 13, + "seed_peak_bits": 0, + "state_bits": 30 + }, + { + "bits": 2, + "length": 13, + "max_squared_drift_ratio_log2": -23, + "modulus": "6", + "numerator": "4", + "output_digest": "576d6ccf89d2465522823a649042d386faf119c8816f3611f8350dc4437d9ab0", + "sample_stride": 1, + "samples": 13, + "seed_peak_bits": 6350, + "state_bits": 30 + }, + { + "bits": 2, + "length": 13, + "max_squared_drift_ratio_log2": -23, + "modulus": "6", + "numerator": "5", + "output_digest": "cb15509754a72ff28734502fa0cdda096e57ec605f9f72bdef866dc85345bdb0", + "sample_stride": 1, + "samples": 13, + "seed_peak_bits": 6350, + "state_bits": 30 + }, + { + "bits": 2, + "length": 15, + "max_squared_drift_ratio_log2": null, + "modulus": "7", + "numerator": "0", + "output_digest": "c083e34f39b224159508d44898f75084de1d10a0de9f35632579eb23c03b1504", + "sample_stride": 1, + "samples": 15, + "seed_peak_bits": 0, + "state_bits": 32 + }, + { + "bits": 2, + "length": 15, + "max_squared_drift_ratio_log2": -21, + "modulus": "7", + "numerator": "1", + "output_digest": "2266f64e9865e5a386deffa26d89a89efe729627dd76f624c6ce4e8bee46450d", + "sample_stride": 1, + "samples": 15, + "seed_peak_bits": 6730, + "state_bits": 32 + }, + { + "bits": 2, + "length": 15, + "max_squared_drift_ratio_log2": -23, + "modulus": "7", + "numerator": "2", + "output_digest": "974f1e8b67ea754c57225e226a4048594e70c436d1aa9cf61b9409692c81c3ce", + "sample_stride": 1, + "samples": 15, + "seed_peak_bits": 6815, + "state_bits": 32 + }, + { + "bits": 2, + "length": 15, + "max_squared_drift_ratio_log2": -21, + "modulus": "7", + "numerator": "3", + "output_digest": "56e493ee76f6eedc650bc6b2f0f07c6401e84346fc1c77d3726c20513d836ceb", + "sample_stride": 1, + "samples": 15, + "seed_peak_bits": 6815, + "state_bits": 32 + }, + { + "bits": 2, + "length": 15, + "max_squared_drift_ratio_log2": -21, + "modulus": "7", + "numerator": "4", + "output_digest": "a0b74dc4e52e4f3031600a9e7068830e8db5c9a74daa38af1d5e1f4d725dd951", + "sample_stride": 1, + "samples": 15, + "seed_peak_bits": 6815, + "state_bits": 32 + }, + { + "bits": 2, + "length": 15, + "max_squared_drift_ratio_log2": -23, + "modulus": "7", + "numerator": "5", + "output_digest": "0842569b0b6c695ef82891931cb26d69b70393b55b37c0f3f2997fb543a0091b", + "sample_stride": 1, + "samples": 15, + "seed_peak_bits": 6815, + "state_bits": 32 + }, + { + "bits": 2, + "length": 15, + "max_squared_drift_ratio_log2": -21, + "modulus": "7", + "numerator": "6", + "output_digest": "28a1ace41bef8749132ad13e1fd1425832744ba21762bc3b58edebac258abdec", + "sample_stride": 1, + "samples": 15, + "seed_peak_bits": 6730, + "state_bits": 32 + }, + { + "bits": 2, + "length": 17, + "max_squared_drift_ratio_log2": null, + "modulus": "8", + "numerator": "0", + "output_digest": "d7cb1b4bb16f99029f8f06ec08d366c76b07287d94ec16dac6524a575ea38102", + "sample_stride": 1, + "samples": 17, + "seed_peak_bits": 0, + "state_bits": 32 + }, + { + "bits": 2, + "length": 17, + "max_squared_drift_ratio_log2": -29, + "modulus": "8", + "numerator": "1", + "output_digest": "a0f30e7619f12b5d6bda1d17283c3d999f8448a2d61556f054000665586a125f", + "sample_stride": 1, + "samples": 17, + "seed_peak_bits": 6619, + "state_bits": 32 + }, + { + "bits": 2, + "length": 17, + "max_squared_drift_ratio_log2": -202, + "modulus": "8", + "numerator": "2", + "output_digest": "682d69a99f6ab79e8ae7f355d1a5b619cc9ec931ce3856fb2bfe023592f85a54", + "sample_stride": 1, + "samples": 17, + "seed_peak_bits": 0, + "state_bits": 32 + }, + { + "bits": 2, + "length": 17, + "max_squared_drift_ratio_log2": -29, + "modulus": "8", + "numerator": "3", + "output_digest": "14f8b2f48ee296be9f83d24f0ba59f24d22a1cdc98962b5b3f540c59bb4852f2", + "sample_stride": 1, + "samples": 17, + "seed_peak_bits": 6619, + "state_bits": 32 + }, + { + "bits": 2, + "length": 17, + "max_squared_drift_ratio_log2": -201, + "modulus": "8", + "numerator": "4", + "output_digest": "84adb7daa70669a6a28549b19ce2bbebf02af1c7ba3a090c5eae6082872b3ea8", + "sample_stride": 1, + "samples": 17, + "seed_peak_bits": 0, + "state_bits": 32 + }, + { + "bits": 2, + "length": 17, + "max_squared_drift_ratio_log2": -29, + "modulus": "8", + "numerator": "5", + "output_digest": "5cfd5d7e8df1d28c07ff43a40a1bb681c4ddadfe5bf9af892cf2f79523226118", + "sample_stride": 1, + "samples": 17, + "seed_peak_bits": 6619, + "state_bits": 32 + }, + { + "bits": 2, + "length": 17, + "max_squared_drift_ratio_log2": -198, + "modulus": "8", + "numerator": "6", + "output_digest": "d2ab24687a1f6186cdebc6189ebaf4468e195a5c635dbed20ad0d310514d3925", + "sample_stride": 1, + "samples": 17, + "seed_peak_bits": 0, + "state_bits": 32 + }, + { + "bits": 2, + "length": 17, + "max_squared_drift_ratio_log2": -29, + "modulus": "8", + "numerator": "7", + "output_digest": "fd5807bf4bb409ad36e10ad40b7fd191e87779fa93f918951a7b82b610827148", + "sample_stride": 1, + "samples": 17, + "seed_peak_bits": 6619, + "state_bits": 32 + }, + { + "bits": 2, + "length": 19, + "max_squared_drift_ratio_log2": null, + "modulus": "9", + "numerator": "0", + "output_digest": "ccac1ca9acb94b91457c13f69c23e2862d15c2e280938fce8d15080a24e8d6ae", + "sample_stride": 1, + "samples": 19, + "seed_peak_bits": 0, + "state_bits": 32 + }, + { + "bits": 2, + "length": 19, + "max_squared_drift_ratio_log2": -22, + "modulus": "9", + "numerator": "1", + "output_digest": "726b36dabc891ff12d84290233b32b2ee144a63de41dec968c38cbbd98b52837", + "sample_stride": 1, + "samples": 19, + "seed_peak_bits": 6768, + "state_bits": 32 + }, + { + "bits": 2, + "length": 19, + "max_squared_drift_ratio_log2": -23, + "modulus": "9", + "numerator": "2", + "output_digest": "ffa827ad4b5cb014ed09aa12fd3065401d2f812ed6832fe774c969486d7c886f", + "sample_stride": 1, + "samples": 19, + "seed_peak_bits": 6768, + "state_bits": 32 + }, + { + "bits": 2, + "length": 19, + "max_squared_drift_ratio_log2": -24, + "modulus": "9", + "numerator": "3", + "output_digest": "71911bd88dcd52e563c67081e0337da5b2b1b0d5bf7467a9d7f28ec475743d08", + "sample_stride": 1, + "samples": 19, + "seed_peak_bits": 6593, + "state_bits": 32 + }, + { + "bits": 2, + "length": 19, + "max_squared_drift_ratio_log2": -20, + "modulus": "9", + "numerator": "4", + "output_digest": "22050542881f0469729205c1d47f0777cc0f820cb20ab8d3431e4fd4b88c7a66", + "sample_stride": 1, + "samples": 19, + "seed_peak_bits": 6740, + "state_bits": 32 + }, + { + "bits": 2, + "length": 19, + "max_squared_drift_ratio_log2": -20, + "modulus": "9", + "numerator": "5", + "output_digest": "783f74b70a5f6b0abc8824c2b15f6acb197e98969cf07149b7831097e9887e29", + "sample_stride": 1, + "samples": 19, + "seed_peak_bits": 6740, + "state_bits": 32 + }, + { + "bits": 2, + "length": 19, + "max_squared_drift_ratio_log2": -24, + "modulus": "9", + "numerator": "6", + "output_digest": "a9793b6f303ffc11264f71c10bd8d6de759c3cc825f5fa4b18d425da55974473", + "sample_stride": 1, + "samples": 19, + "seed_peak_bits": 6593, + "state_bits": 32 + }, + { + "bits": 2, + "length": 19, + "max_squared_drift_ratio_log2": -24, + "modulus": "9", + "numerator": "7", + "output_digest": "4654b81a75c70a58102dcdb81f514a66658887e5fd0a2f84fca58a2db3263181", + "sample_stride": 1, + "samples": 19, + "seed_peak_bits": 6768, + "state_bits": 32 + }, + { + "bits": 2, + "length": 19, + "max_squared_drift_ratio_log2": -22, + "modulus": "9", + "numerator": "8", + "output_digest": "9d38a15cc56ca42184f9c8dd33bf0f3d403b012736b04d613620ce621726e9d5", + "sample_stride": 1, + "samples": 19, + "seed_peak_bits": 6768, + "state_bits": 32 + }, + { + "bits": 24, + "length": 3, + "max_squared_drift_ratio_log2": null, + "modulus": "1", + "numerator": "0", + "output_digest": "4c36e53f87d9bafe6fda9826af4bc19e550196a557a8b4fa9d940a6d383a8f21", + "sample_stride": 1, + "samples": 3, + "seed_peak_bits": 0, + "state_bits": 50 + }, + { + "bits": 24, + "length": 5, + "max_squared_drift_ratio_log2": null, + "modulus": "2", + "numerator": "0", + "output_digest": "18494cd0341be801ba4227e281b866afdb493eb8df2294a59a9c8326bf0ccb6c", + "sample_stride": 1, + "samples": 5, + "seed_peak_bits": 0, + "state_bits": 50 + }, + { + "bits": 24, + "length": 5, + "max_squared_drift_ratio_log2": -201, + "modulus": "2", + "numerator": "1", + "output_digest": "5654b9211a6fa988f0d8931c6477ada119cd273021c649345f1ba37bfb1f16fc", + "sample_stride": 1, + "samples": 5, + "seed_peak_bits": 0, + "state_bits": 50 + }, + { + "bits": 24, + "length": 7, + "max_squared_drift_ratio_log2": null, + "modulus": "3", + "numerator": "0", + "output_digest": "fa1196b752531064b0048aacd62a68b2b50c53a7ba62add752e29b881abd7bd2", + "sample_stride": 1, + "samples": 7, + "seed_peak_bits": 0, + "state_bits": 52 + }, + { + "bits": 24, + "length": 7, + "max_squared_drift_ratio_log2": -24, + "modulus": "3", + "numerator": "1", + "output_digest": "96fe23b6159e5fd1a8ef0a3f477c390736fcbcbd041f68d6a5f790456645c15a", + "sample_stride": 1, + "samples": 7, + "seed_peak_bits": 12085, + "state_bits": 52 + }, + { + "bits": 24, + "length": 7, + "max_squared_drift_ratio_log2": -24, + "modulus": "3", + "numerator": "2", + "output_digest": "2d79790ec97c35a2128712c119894393b132e15c090f92b8937d26b6d4806d0f", + "sample_stride": 1, + "samples": 7, + "seed_peak_bits": 12085, + "state_bits": 52 + }, + { + "bits": 24, + "length": 9, + "max_squared_drift_ratio_log2": null, + "modulus": "4", + "numerator": "0", + "output_digest": "4ad9798e005fe636996e2a410fc9872532f5358b364753f08bc316b2b4d6bc7b", + "sample_stride": 1, + "samples": 9, + "seed_peak_bits": 0, + "state_bits": 52 + }, + { + "bits": 24, + "length": 9, + "max_squared_drift_ratio_log2": -202, + "modulus": "4", + "numerator": "1", + "output_digest": "dac4346d4498dc45c606279ce79e22575960ed9e72066210e3af408204f06abc", + "sample_stride": 1, + "samples": 9, + "seed_peak_bits": 0, + "state_bits": 52 + }, + { + "bits": 24, + "length": 9, + "max_squared_drift_ratio_log2": -201, + "modulus": "4", + "numerator": "2", + "output_digest": "7919e87f49b9c8cfb4b92df85d50b15e10fdba8b1a8f59cec32421f3dc9cd03a", + "sample_stride": 1, + "samples": 9, + "seed_peak_bits": 0, + "state_bits": 52 + }, + { + "bits": 24, + "length": 9, + "max_squared_drift_ratio_log2": -198, + "modulus": "4", + "numerator": "3", + "output_digest": "9461635c87a884b700a7accad4b7e30e125eaa7b4a893b942fd923f57e588f00", + "sample_stride": 1, + "samples": 9, + "seed_peak_bits": 0, + "state_bits": 52 + }, + { + "bits": 24, + "length": 11, + "max_squared_drift_ratio_log2": null, + "modulus": "5", + "numerator": "0", + "output_digest": "d7720b0628de294608435ccfdd0950397cafded1bff1ed46ddae1ba3de718bb3", + "sample_stride": 1, + "samples": 11, + "seed_peak_bits": 0, + "state_bits": 52 + }, + { + "bits": 24, + "length": 11, + "max_squared_drift_ratio_log2": -23, + "modulus": "5", + "numerator": "1", + "output_digest": "cbc2bab9e583c60eb99e004f22bb65e01498bf7332af4309e48bb5d798127ee6", + "sample_stride": 1, + "samples": 11, + "seed_peak_bits": 12271, + "state_bits": 52 + }, + { + "bits": 24, + "length": 11, + "max_squared_drift_ratio_log2": -22, + "modulus": "5", + "numerator": "2", + "output_digest": "a622acf40c13f6a3ff05531e84d261e1f1889695e38e94f8dfafaa5d921b2c55", + "sample_stride": 1, + "samples": 11, + "seed_peak_bits": 12271, + "state_bits": 52 + }, + { + "bits": 24, + "length": 11, + "max_squared_drift_ratio_log2": -22, + "modulus": "5", + "numerator": "3", + "output_digest": "9dae8cdb9ecefae4c9937ef582950aecf1f9d84222f9dfb85e15d21f7f0d88e3", + "sample_stride": 1, + "samples": 11, + "seed_peak_bits": 12271, + "state_bits": 52 + }, + { + "bits": 24, + "length": 11, + "max_squared_drift_ratio_log2": -23, + "modulus": "5", + "numerator": "4", + "output_digest": "d07cc4ce2edd85ba95e52416c21687533784c8521ffd08acfa27b1ebdac938de", + "sample_stride": 1, + "samples": 11, + "seed_peak_bits": 12271, + "state_bits": 52 + }, + { + "bits": 24, + "length": 13, + "max_squared_drift_ratio_log2": null, + "modulus": "6", + "numerator": "0", + "output_digest": "c3956986df0e5040317391f16deb419e20e3646a9e5d16e8e5ddba6add70cad7", + "sample_stride": 1, + "samples": 13, + "seed_peak_bits": 0, + "state_bits": 52 + }, + { + "bits": 24, + "length": 13, + "max_squared_drift_ratio_log2": -24, + "modulus": "6", + "numerator": "1", + "output_digest": "aa73c2d947849093550b78f6842bfec280ce7081ef13b0aee1e6471b275babf4", + "sample_stride": 1, + "samples": 13, + "seed_peak_bits": 12085, + "state_bits": 52 + }, + { + "bits": 24, + "length": 13, + "max_squared_drift_ratio_log2": -24, + "modulus": "6", + "numerator": "2", + "output_digest": "2f26218f875176c57b78aaddccf5670b3034ef9845c458b7f89466fb92b0d589", + "sample_stride": 1, + "samples": 13, + "seed_peak_bits": 12085, + "state_bits": 52 + }, + { + "bits": 24, + "length": 13, + "max_squared_drift_ratio_log2": -201, + "modulus": "6", + "numerator": "3", + "output_digest": "02a7e76bd528302dda68f8f32f521b1837ee265aa249af52bfece45abba9c695", + "sample_stride": 1, + "samples": 13, + "seed_peak_bits": 0, + "state_bits": 52 + }, + { + "bits": 24, + "length": 13, + "max_squared_drift_ratio_log2": -24, + "modulus": "6", + "numerator": "4", + "output_digest": "4f7ece1ea1ac8bd4ce232c67555d5575c2aaf1b21d820ac251acc09370ae5a07", + "sample_stride": 1, + "samples": 13, + "seed_peak_bits": 12085, + "state_bits": 52 + }, + { + "bits": 24, + "length": 13, + "max_squared_drift_ratio_log2": -24, + "modulus": "6", + "numerator": "5", + "output_digest": "5a7b4941908450352a0d5c87f5fc57cfaecc76f00a32b02153ad2719688bffc3", + "sample_stride": 1, + "samples": 13, + "seed_peak_bits": 12085, + "state_bits": 52 + }, + { + "bits": 24, + "length": 15, + "max_squared_drift_ratio_log2": null, + "modulus": "7", + "numerator": "0", + "output_digest": "229a423b9eccb1a3309ec56ac8f337792a3acb1e5b3838f158fbc13f64ec1be8", + "sample_stride": 1, + "samples": 15, + "seed_peak_bits": 0, + "state_bits": 54 + }, + { + "bits": 24, + "length": 15, + "max_squared_drift_ratio_log2": -23, + "modulus": "7", + "numerator": "1", + "output_digest": "6a1e6f852f471a27c3b58430ccc87bb2b636edccb996d05c5468c6b7f34b7639", + "sample_stride": 1, + "samples": 15, + "seed_peak_bits": 12846, + "state_bits": 54 + }, + { + "bits": 24, + "length": 15, + "max_squared_drift_ratio_log2": -21, + "modulus": "7", + "numerator": "2", + "output_digest": "5751caacea72d58dde39a5b1be16fc067076bc87751376bfa8e10ed9e7cc0e64", + "sample_stride": 1, + "samples": 15, + "seed_peak_bits": 12911, + "state_bits": 54 + }, + { + "bits": 24, + "length": 15, + "max_squared_drift_ratio_log2": -22, + "modulus": "7", + "numerator": "3", + "output_digest": "467b681d42dab00aa4b2606d0a230dafc1955d84bae75ead130802da0f0d8894", + "sample_stride": 1, + "samples": 15, + "seed_peak_bits": 12911, + "state_bits": 54 + }, + { + "bits": 24, + "length": 15, + "max_squared_drift_ratio_log2": -22, + "modulus": "7", + "numerator": "4", + "output_digest": "cc0e941f32c9d7388dcc2d33bb5a41975371b77029519062c822838e333453fe", + "sample_stride": 1, + "samples": 15, + "seed_peak_bits": 12911, + "state_bits": 54 + }, + { + "bits": 24, + "length": 15, + "max_squared_drift_ratio_log2": -21, + "modulus": "7", + "numerator": "5", + "output_digest": "f8f6974e9e27ce3a63581ae1c55a7467f6703e486af56df01b25eef198b5dfc8", + "sample_stride": 1, + "samples": 15, + "seed_peak_bits": 12911, + "state_bits": 54 + }, + { + "bits": 24, + "length": 15, + "max_squared_drift_ratio_log2": -23, + "modulus": "7", + "numerator": "6", + "output_digest": "dc65c9921de3dfd62442311b1816953c41600a63f1972d4b0d9c747432a22177", + "sample_stride": 1, + "samples": 15, + "seed_peak_bits": 12846, + "state_bits": 54 + }, + { + "bits": 24, + "length": 17, + "max_squared_drift_ratio_log2": null, + "modulus": "8", + "numerator": "0", + "output_digest": "457d7d393bc99d32b26abc0f398b424e40d9f09b8015ac367511e3a9c2d187a6", + "sample_stride": 1, + "samples": 17, + "seed_peak_bits": 0, + "state_bits": 54 + }, + { + "bits": 24, + "length": 17, + "max_squared_drift_ratio_log2": -24, + "modulus": "8", + "numerator": "1", + "output_digest": "be9471d919467f5a165ab4bf7610a860d3d77197865691526a93f172ffc0df6b", + "sample_stride": 1, + "samples": 17, + "seed_peak_bits": 12849, + "state_bits": 54 + }, + { + "bits": 24, + "length": 17, + "max_squared_drift_ratio_log2": -202, + "modulus": "8", + "numerator": "2", + "output_digest": "836ede0a2d80c5b6f1dd30b2cc968a1003c687aa303e44d1c7ed4e72deab3358", + "sample_stride": 1, + "samples": 17, + "seed_peak_bits": 0, + "state_bits": 54 + }, + { + "bits": 24, + "length": 17, + "max_squared_drift_ratio_log2": -24, + "modulus": "8", + "numerator": "3", + "output_digest": "5d7c9f7e4269c88546a0916efe7e4bd74357b9ec89a833acccf4ed70737ee26f", + "sample_stride": 1, + "samples": 17, + "seed_peak_bits": 12849, + "state_bits": 54 + }, + { + "bits": 24, + "length": 17, + "max_squared_drift_ratio_log2": -201, + "modulus": "8", + "numerator": "4", + "output_digest": "b67f45ede26eadf29264f6a2cb7ff5d4cba9fe4d093d63071ce6adc1fa4b1243", + "sample_stride": 1, + "samples": 17, + "seed_peak_bits": 0, + "state_bits": 54 + }, + { + "bits": 24, + "length": 17, + "max_squared_drift_ratio_log2": -24, + "modulus": "8", + "numerator": "5", + "output_digest": "545aff19ba054eea9aa0bd3b44b1e1cd58a951f450e840bec267b507ffdb038d", + "sample_stride": 1, + "samples": 17, + "seed_peak_bits": 12849, + "state_bits": 54 + }, + { + "bits": 24, + "length": 17, + "max_squared_drift_ratio_log2": -198, + "modulus": "8", + "numerator": "6", + "output_digest": "869dc83c5ad11398e8e6bd2a411f945d36b93b83c6068d16f3784e2d7c8debc6", + "sample_stride": 1, + "samples": 17, + "seed_peak_bits": 0, + "state_bits": 54 + }, + { + "bits": 24, + "length": 17, + "max_squared_drift_ratio_log2": -24, + "modulus": "8", + "numerator": "7", + "output_digest": "a874fa1d238169ec9743e426b95ca0c7c5f3ecdc415b0d70ca89063fa1782b6e", + "sample_stride": 1, + "samples": 17, + "seed_peak_bits": 12849, + "state_bits": 54 + }, + { + "bits": 24, + "length": 19, + "max_squared_drift_ratio_log2": null, + "modulus": "9", + "numerator": "0", + "output_digest": "ecfe92585edcd0030c0f16a0f6afa54857a02a73300c416879c9691296f3a295", + "sample_stride": 1, + "samples": 19, + "seed_peak_bits": 0, + "state_bits": 54 + }, + { + "bits": 24, + "length": 19, + "max_squared_drift_ratio_log2": -21, + "modulus": "9", + "numerator": "1", + "output_digest": "55b8e9b8fb46462d8b3d350232516359bda301f66711017aeca90eba8cd4b487", + "sample_stride": 1, + "samples": 19, + "seed_peak_bits": 12911, + "state_bits": 54 + }, + { + "bits": 24, + "length": 19, + "max_squared_drift_ratio_log2": -22, + "modulus": "9", + "numerator": "2", + "output_digest": "fee5894bd4f65269cb8aa8d4018f02b14051977dff629f6cce6f677989a8d7cc", + "sample_stride": 1, + "samples": 19, + "seed_peak_bits": 12911, + "state_bits": 54 + }, + { + "bits": 24, + "length": 19, + "max_squared_drift_ratio_log2": -26, + "modulus": "9", + "numerator": "3", + "output_digest": "8ea3bfb5a106425a3ae1f4f001e30d6dbec598cd2db5784ecb42ccb9a73272b9", + "sample_stride": 1, + "samples": 19, + "seed_peak_bits": 12812, + "state_bits": 54 + }, + { + "bits": 24, + "length": 19, + "max_squared_drift_ratio_log2": -22, + "modulus": "9", + "numerator": "4", + "output_digest": "3706ddbb5b6da82bf88389dcfcd05c907c9ec29f9390aacafdc026d4bd6c5d48", + "sample_stride": 1, + "samples": 19, + "seed_peak_bits": 12882, + "state_bits": 54 + }, + { + "bits": 24, + "length": 19, + "max_squared_drift_ratio_log2": -22, + "modulus": "9", + "numerator": "5", + "output_digest": "ca0b66af502ca3907933116636695a7d5d9ccb9e9867e9a4342d09625acc274a", + "sample_stride": 1, + "samples": 19, + "seed_peak_bits": 12882, + "state_bits": 54 + }, + { + "bits": 24, + "length": 19, + "max_squared_drift_ratio_log2": -26, + "modulus": "9", + "numerator": "6", + "output_digest": "96975d81999568a85849c0b45e66b1b0b2c4e0fca0be0b3c8c0553aa74017691", + "sample_stride": 1, + "samples": 19, + "seed_peak_bits": 12812, + "state_bits": 54 + }, + { + "bits": 24, + "length": 19, + "max_squared_drift_ratio_log2": -23, + "modulus": "9", + "numerator": "7", + "output_digest": "1e2c90133b465df7a1907e12eb1507f1b205138f24f324e139a42d66d53ed127", + "sample_stride": 1, + "samples": 19, + "seed_peak_bits": 12911, + "state_bits": 54 + }, + { + "bits": 24, + "length": 19, + "max_squared_drift_ratio_log2": -21, + "modulus": "9", + "numerator": "8", + "output_digest": "3e00026f94368013c9d08658ea4e98610a69000f4d7f9e7b69e0f70f15407015", + "sample_stride": 1, + "samples": 19, + "seed_peak_bits": 12911, + "state_bits": 54 + }, + { + "bits": 64, + "length": 33, + "max_squared_drift_ratio_log2": -22, + "modulus": "16", + "numerator": "-3", + "output_digest": "d8382c68ca9234436eb2f687eb65a351fd61d2f0c47ba3bc690c2ac9ea54595e", + "sample_stride": 1, + "samples": 33, + "seed_peak_bits": 30178, + "state_bits": 96 + }, + { + "bits": 64, + "length": 41, + "max_squared_drift_ratio_log2": -21, + "modulus": "17", + "numerator": "5", + "output_digest": "25a6cb4bb6be3f9800dcacfd2bcbc2561462d6d51082a705fe0591c2170eb785", + "sample_stride": 1, + "samples": 41, + "seed_peak_bits": 30229, + "state_bits": 96 + }, + { + "bits": 128, + "length": 65, + "max_squared_drift_ratio_log2": -20, + "modulus": "31", + "numerator": "-17", + "output_digest": "72d0b883ff1a0129b00b2ac26b69900450d5cffb4d5879322bc62d71d680239b", + "sample_stride": 1, + "samples": 65, + "seed_peak_bits": 71414, + "state_bits": 162 + }, + { + "bits": 96, + "length": 9, + "max_squared_drift_ratio_log2": -83, + "modulus": "1461501637330902918203684832716283019655932542983", + "numerator": "1", + "output_digest": "8bec2d4ab545c4237c7f95dd0b4e50866f8321e0dae2b8515d567bc535d453c7", + "sample_stride": 1, + "samples": 9, + "seed_peak_bits": 3427, + "state_bits": 124 + }, + { + "bits": 96, + "length": 9, + "max_squared_drift_ratio_log2": -37, + "modulus": "1461501637330902918203684832716283019655932542983", + "numerator": "1532495540865888858358347027150309183618739122183602185", + "output_digest": "8bec2d4ab545c4237c7f95dd0b4e50866f8321e0dae2b8515d567bc535d453c7", + "sample_stride": 1, + "samples": 9, + "seed_peak_bits": 45526, + "state_bits": 124 + }, + { + "bits": 96, + "length": 11, + "max_squared_drift_ratio_log2": -27, + "modulus": "1461501637330902918203684832716283019655932542976", + "numerator": "182687704666362864775460604089535377456991567871", + "output_digest": "5385c4c3e5bc728cc97998fc0dac03efb458eaea594b09698c6f3f6c2ab202db", + "sample_stride": 1, + "samples": 11, + "seed_peak_bits": 45661, + "state_bits": 124 + }, + { + "bits": 96, + "length": 11, + "max_squared_drift_ratio_log2": -27, + "modulus": "1461501637330902918203684832716283019655932542976", + "numerator": "182687704666362864775460604089535377456991567873", + "output_digest": "5385c4c3e5bc728cc97998fc0dac03efb458eaea594b09698c6f3f6c2ab202db", + "sample_stride": 1, + "samples": 11, + "seed_peak_bits": 45661, + "state_bits": 124 + }, + { + "bits": 24, + "length": 16384, + "max_squared_drift_ratio_log2": -27, + "modulus": "65537", + "numerator": "12345", + "output_digest": "c5e1a14c537554d8abe233e69e04cfab95111eaac4633d880f10867339dba274", + "sample_stride": 257, + "samples": 65, + "seed_peak_bits": 19565, + "state_bits": 74 + } + ], + "schema": 1, + "source_pins": 9, + "words": [ + { + "bits": 52, + "borrowed_arrays": 0, + "calls": 1, + "count": 1, + "dirt": "7", + "divisions": 24, + "length": 1, + "output_digest": "6b1338cd39a81566624207eb91d8ab2497bc6412d9d9538f511668d8b1bbc55d", + "products": 66, + "seed_terms": 0, + "sign": -1, + "stride": 0, + "tape_moves": 15 + }, + { + "bits": 52, + "borrowed_arrays": 0, + "calls": 1, + "count": 1, + "dirt": "7", + "divisions": 24, + "length": 1, + "output_digest": "e511dc9c3ca3e3b436127deaf6143ee8d07d432f9d262b13af5d8b943f92b933", + "products": 66, + "seed_terms": 0, + "sign": 1, + "stride": 0, + "tape_moves": 15 + }, + { + "bits": 79, + "borrowed_arrays": 2, + "calls": 10, + "count": 3, + "dirt": "7", + "divisions": 24, + "length": 1, + "output_digest": "9f991249c6cefe1c06fafa92d0cd865b1707af5a8ab95530de205d2a8b16c939", + "products": 66, + "seed_terms": 0, + "sign": -1, + "stride": 0, + "tape_moves": 128 + }, + { + "bits": 79, + "borrowed_arrays": 2, + "calls": 10, + "count": 3, + "dirt": "7", + "divisions": 24, + "length": 1, + "output_digest": "9273d9eb7540a6d8fd0b0e3696be08b9db2599105e1983ef042691bc0a6be377", + "products": 66, + "seed_terms": 0, + "sign": 1, + "stride": 0, + "tape_moves": 124 + }, + { + "bits": 52, + "borrowed_arrays": 0, + "calls": 1, + "count": 1, + "dirt": "-13", + "divisions": 120, + "length": 5, + "output_digest": "c6f93071b966bb49deeeea90b5d03eaf49d431d49a21efa4b7784f30a116c10b", + "products": 306, + "seed_terms": 4872, + "sign": -1, + "stride": 2, + "tape_moves": 1095 + }, + { + "bits": 52, + "borrowed_arrays": 0, + "calls": 1, + "count": 1, + "dirt": "-13", + "divisions": 120, + "length": 5, + "output_digest": "de0221cb455890a9224c40e9e6df02ca3564b4941cd4842c05d800c14a0421d4", + "products": 306, + "seed_terms": 4872, + "sign": 1, + "stride": 2, + "tape_moves": 1095 + }, + { + "bits": 80, + "borrowed_arrays": 2, + "calls": 10, + "count": 3, + "dirt": "-13", + "divisions": 120, + "length": 5, + "output_digest": "4f44ba568549c499a6ea04200b7a7d2be095773ed2c04e531df4a70a4dc894e4", + "products": 306, + "seed_terms": 6216, + "sign": -1, + "stride": 2, + "tape_moves": 10720 + }, + { + "bits": 80, + "borrowed_arrays": 2, + "calls": 10, + "count": 3, + "dirt": "-13", + "divisions": 120, + "length": 5, + "output_digest": "a134926e8e5025b8618366997c7754fb76aaedc5fd212a0343396da31f2d4510", + "products": 306, + "seed_terms": 6216, + "sign": 1, + "stride": 2, + "tape_moves": 10700 + }, + { + "bits": 88, + "borrowed_arrays": 0, + "calls": 1, + "count": 1, + "dirt": "1099511627776", + "divisions": 192, + "length": 8, + "output_digest": "f6ef7594d1ab74bc07de5f19cffda9d4cbc1d56decd676c9c2fff67b7ea3e47a", + "products": 486, + "seed_terms": 6696, + "sign": -1, + "stride": -3, + "tape_moves": 1816 + }, + { + "bits": 88, + "borrowed_arrays": 0, + "calls": 1, + "count": 1, + "dirt": "1099511627776", + "divisions": 192, + "length": 8, + "output_digest": "85656ff00c50d5def88a62068b8f7acd2849c9c0931c6061171d33befed309b8", + "products": 486, + "seed_terms": 6696, + "sign": 1, + "stride": -3, + "tape_moves": 1816 + }, + { + "bits": 115, + "borrowed_arrays": 2, + "calls": 10, + "count": 3, + "dirt": "1099511627776", + "divisions": 192, + "length": 8, + "output_digest": "10202e0da4332da48ac4db663b3006a43e6f31d4cc597cd24d2d90cd0d1285f9", + "products": 486, + "seed_terms": 7992, + "sign": -1, + "stride": -3, + "tape_moves": 17906 + }, + { + "bits": 115, + "borrowed_arrays": 2, + "calls": 10, + "count": 3, + "dirt": "1099511627776", + "divisions": 192, + "length": 8, + "output_digest": "e3e4fb48d6eef8244f68f8dd58d87e0c899af16fe5c9f7e2d1f3a49e95381e39", + "products": 486, + "seed_terms": 7992, + "sign": 1, + "stride": -3, + "tape_moves": 17874 + }, + { + "bits": 128, + "borrowed_arrays": 0, + "calls": 1, + "count": 1, + "dirt": "1208925819614629174706176", + "divisions": 408, + "length": 17, + "output_digest": "2f21573d8df58be1d58ecc653f505193fbf281e2db3ce00e000eeea6362c11b6", + "products": 1026, + "seed_terms": 8712, + "sign": -1, + "stride": 7, + "tape_moves": 5743 + }, + { + "bits": 128, + "borrowed_arrays": 0, + "calls": 1, + "count": 1, + "dirt": "1208925819614629174706176", + "divisions": 408, + "length": 17, + "output_digest": "b281d30d805b8ae363cafcf58bba0e9694d76e6072e9a1b7dc6c32b147190dd4", + "products": 1026, + "seed_terms": 8712, + "sign": 1, + "stride": 7, + "tape_moves": 5743 + }, + { + "bits": 155, + "borrowed_arrays": 2, + "calls": 10, + "count": 3, + "dirt": "1208925819614629174706176", + "divisions": 408, + "length": 17, + "output_digest": "1d1b80b30532c66fc68ba8f9f0e1556ea056879d94f527adac6612a18dfeac1e", + "products": 1026, + "seed_terms": 10008, + "sign": -1, + "stride": 7, + "tape_moves": 56576 + }, + { + "bits": 155, + "borrowed_arrays": 2, + "calls": 10, + "count": 3, + "dirt": "1208925819614629174706176", + "divisions": 408, + "length": 17, + "output_digest": "78282d9db45404942ecce203207fd4ad13ca26d4886415ba60109e505a912fef", + "products": 1026, + "seed_terms": 10008, + "sign": 1, + "stride": 7, + "tape_moves": 56508 + } + ] +} diff --git a/research/dirty-stride-fast-phases/NOTICE.md b/research/dirty-stride-fast-phases/NOTICE.md new file mode 100644 index 0000000000..858a4c6138 --- /dev/null +++ b/research/dirty-stride-fast-phases/NOTICE.md @@ -0,0 +1,32 @@ +# Attribution and license + +Prepared for Ali Romman (@aliromman), with substantial OpenAI Codex +assistance in research, proof, implementation, verification, and writing. +One assigned AI agent prepared this package. Parent-agent review is +separate from independent human review. New material is Apache-2.0. + +Ali Romman's PR1038 supplies the dirty quadratic-phase interface, three +rounded real shears, exact inverse, dirty-target error proof, independent +interval oracle, and test-data and tape-observer patterns used here. +`chirp.py` derives its dirty wrappers and precision function from that +source. `DERIVATION.diff` records the change. The original files remain +byte-exact under `vendor/pr1038/`, with Git-blob and SHA256 pins. + +The retained PR985 solver credits Rohan Guliani's PR846 Euclidean +decomposition, Charles Bethin's PR764 suffix estimate, and Chafik +Boukhalfa's PR628 dirty kernel lineage. Chafik Boukhalfa's PR479 supplies +the retained chirp family and scalar cost context. Ross Richardson's +PR948 gives independent symbolic evidence for the predecessor phase +stream. The original notices remain in the package. + +The PR1054 constant-coefficient precision result and PR1025 analytic +root-height result were read for their scope. Neither is an executed +dependency or an assumed combined theorem here. The previous local +cyclic-stability work separates raw coefficient growth from necessary +precision; this package makes no contrary obstruction claim. + +Multiplicative quadratic-phase recurrences, Machin series, Taylor series, +and three-shear plane rotations are established methods. No literature +priority is claimed. This contribution supplies their explicit reversible +stream, uniform guard bound, exact cleanup, and paid integration with the +specific dirty solver. diff --git a/research/dirty-stride-fast-phases/OVERLAP.md b/research/dirty-stride-fast-phases/OVERLAP.md new file mode 100644 index 0000000000..1af0042370 --- /dev/null +++ b/research/dirty-stride-fast-phases/OVERLAP.md @@ -0,0 +1,37 @@ +# Source and overlap refresh + +The experimental base remains +`abdc6a7002abdbcd9e21bfc5e1a795784deacfdf`. Its selected PR942 value was +`1661960547713947/10^18`, with the complex supplier binding. The source +manifests do not pin the moving root selected-result, current-status, or +PR-ledger files. + +On 2026-10-13 at 18:49 UTC, read-only GitHub API checks confirmed these +open predecessor heads: + +| PR | Head | Relation | +| --- | --- | --- | +| 1038 | `99c71d385171480ecc2c43287aa4a04e1027fcd9` | Direct phase and dirty-wrapper dependency | +| 1054 | `6bc675be1554b3a3be771147f21f9dc2c6dae50a` | Constant-coefficient precision; separate theorem | +| 1085 | `82b99c7d12dac009d35636562df4c6c91e885409` | Ideal dirty-prefix row bound; no phase-stream implementation | + +The coordinator's open-PR snapshot was searched for dirty, chirp, and +stream titles. It includes PR479, PR948, PR985, PR1038, PR1054, and +PR1057. The final API query inspected the 50 open PRs most recently +updated and searched their titles for phase, chirp, dirty, and stream. +PR1085 was the only returned new relevant title. Its body was read. +Its phase-access obligation stays open. This is a bounded overlap check, +not a literature-priority claim or a search of every external result. + +Live main changed during the run. At commit +`ba82555c64e1879e4b7decce28186f208fcb8262`, the root selected record gives +PR1042, `kappa=1672188119205421/10^18`, with review +`docs/research/community-batch1047-review.md` and audited source head +`7e8ef087c02a77038a28a337a7ee3a77039cac0d`. This was read directly through +the GitHub contents API. It is a separate live-status note. The experiments +keep their exact base and dependency pins. They make no numeric frontier +claim against either selected record. + +There was no pinned sprint issue in the supplied coordination state. +The local claim is `CLAIM.md`. No public comment, push, or PR was made +by this worker. The main agent owns publication after review. diff --git a/research/dirty-stride-fast-phases/PROOF.md b/research/dirty-stride-fast-phases/PROOF.md new file mode 100644 index 0000000000..fae6799b2f --- /dev/null +++ b/research/dirty-stride-fast-phases/PROOF.md @@ -0,0 +1,338 @@ +# A reversible quadratic phase stream + +## Contract and source boundary + +Let `L` be a positive integer. Let `a` be an integer and let `m>0` be an +integer. For `0<=j=2`. The phase state grid is `e=2^-Q`, +where + + b = bit_length(L+1), + Q = W+2b+20. + +There are two complex state values. Their forward update and their reverse +update are exact inverse permutations on the grid. After a pass, the code +runs all state updates backward. It then clears both initial state values. +An inverse data pass generates the same coefficients in forward index order +and reverses each data rotation. It does not invert a rounded complex +product, normalize a complex state, or infer a phase from the data. + +`chirp.py` retains the PR1038 leaf, factor-word, precision, and data-layout +interfaces. It changes phase preparation, phase passes, and phase cleanup. +The source pin is `99c71d385171480ecc2c43287aa4a04e1027fcd9`. The exact +retained bytes are under `vendor/pr1038/`. This proof uses PR1038 Sections +4--7 for the grid rotation, dirty conjugation, and inherited PR985 word. +It does not assume PR1054's improved precision. It keeps the original +`ceil(log2 n)` borrowed arrays for `n` factors. + +## 1. Bounded rotation coefficients + +For an exact unit phase `v`, select a quarter-turn `q` so that the real +component of `i^-q v` is positive. Write that reduced phase as `(c,s)`. +Its three real shear coefficients are `(t,s,t)`, with `t=s/(1+c)`. +The corresponding matrix is exactly the reduced rotation. The product +of the three shears is followed by `i^q`. + +The grid map rounds the product in each shear to the data grid. It adds +or subtracts that result exactly. Its inverse reverses the instruction +order and each addition. A shear source does not change in that shear. +Thus both directions use the same rounded product, including ties. The +map is an exact grid permutation for every input magnitude. + +We use only `|t|<=1/2` and `|s|<=1`. Each shear matrix has Euclidean norm +at most 2. If both computed coefficients differ from their exact values +by at most `delta`, the difference of the matrix products is at most +`12 delta`. The three grid roundings contribute at most +`(4+2+1)u/2 < 4u`. Thus coefficients within `2^-(W+8)` give + + |S(v)-D v| <= 4(|v|+1)u, + |S^-1(v)-D^-1 v| <= 4(|v|+1)u. (1) + +Here the two uses of `v` denote an arbitrary data value, and `D` is the +unit phase. The same bounds apply to the reverse word. These are the +PR1038 bounds. They include the initial dirty value. + +## 2. Read a phase state without normalization + +Suppose that a dyadic state `v` is within Euclidean distance `E<=1/64` +of a unit phase `v0`. `chart` selects the largest component in absolute +value and rotates its sign to the positive real axis. Ties have a fixed +order. Let `(c,s)=i^-q v` and `(c0,s0)=i^-q v0`. Then + + c >= |s|, c >= (1-E)/sqrt(2), + c0 >= (2/3)(1-E)-E >= 41/64 > 3/5. + +The exact reduced tangent satisfies `|s0/(1+c0)|<=1/2`. Compute +`s/(1+c)`, round it to spacing `delta`, and restrict it to `[-1/2,1/2]`. +Round `s` to that spacing and restrict it to `[-1,1]`. Restriction to an +interval cannot increase the error from a point in that interval. +The two denominators are at least 1. Therefore + + |t_hat-t0| <= 2E+delta/2, + |s_hat-s0| <= E+delta/2. (2) + +The quadrant need not agree with the exact phase's nearest quadrant. +The lower bound on `c0` applies to the quadrant actually selected. Thus +the proof also covers changes of chart near diagonal rays. + +The chart adds coefficients to a clean destination. After the rotation, +the code recomputes and subtracts those coefficients before it changes +their phase source. Both generation and erasure pay their division. + +At state precision `Q`, use `delta=e/256`. The matrix error from (2) is +at most `4(2(2E+delta/2)+(E+delta/2)) <= 20E+e`. +For a state input of modulus at most 2, one rounded update therefore has +error at most `40E+6e` relative to multiplication by `v0`. + +## 3. Paid seed preparation and exact cleanup + +The code prepares coefficients for `exp(2 pi i a/m)` and `omega` only. +Let `C=Q+8` and `P=C+16`. For one phase residue `A=a mod m`, set + + q = floor(4A/m+1/2), + h = 4A-qm, + theta = pi h/(2m). + +Then `|theta|<=pi/4`. A zero `h` gives exact zero shear coefficients. +Otherwise evaluate the first `P` terms of both Machin arctangent series +as exact rationals. Include the weights 16 and -4. The alternating tail +bound gives error below `20*2^-P`. Round the exact sum once to `2^-P`. +Round its product with `h/(2m)` once more. The angle error is less than +`6*2^-P`, and the rounded angle has absolute value less than 1. + +Evaluate sine and cosine with `P` exact rational Taylor terms at this +angle. The omitted tails are less than `2^-P`, because +`(2P)!>2^P`. Each total error is less than `7*2^-P`. Both the exact and +computed denominators in `sine/(1+cosine)` exceed 1. The tangent error +before final rounding is less than `14*2^-P`. Round both coefficients +to `2^-C`. Since `P=C+16`, both errors are less than `2^-C`. Their +magnitudes satisfy the bounds in Section 1. + +The exact series are not dropped. Each Machin power update is an exact +integer multiplication by 25 or 57121. Exact division reverses it. +Subtract the arctangent terms in reverse order to clear the sum. For a +nonzero rounded angle `X`, the Taylor term update is multiplication by +`-X^2/((2k)(2k+1))` or `-X^2/((2k-1)(2k))`. It is invertible over the +rationals. Subtract the last term from its sum, invert the term update, +and continue backward. Clear the initial terms `X` and 1. If `X=0`, +skip these term updates. This branch also covers a nonzero residue whose +angle rounds to zero. Finally, recompute and subtract the rounded angle +and rounded pi, and reverse the two Machin sums. The code checks zero +at each cleanup boundary. + +The coefficient output is a reversible addition to its destination. +Preparation adds the coefficients; erasure subtracts the same values. +In both cases all exact series state is cleared. To prepare `r_hat_0`, +apply the first seed's grid rotation to 1, then erase those coefficients. +The fixed coefficients of `omega` remain for stream use. At the end, +prepare the first seed's coefficients again, reverse its rotation, clear +the resulting 1, and erase the first and second seed coefficients. + +This uses a fixed number of rational arithmetic registers. Their widths +are polynomial in `P`; they are not all `O(P)`-bit registers. This change +is explicit in the work and workspace bill below. + +Directly rounding a complex product does not provide the required inverse +on the full data grid. At angle `pi/4`, both `(u,0)` and `(2u,0)` round +to `(u,u)`. The three-shear map avoids this collision by retaining each +shear source until its controlled addition is complete. This observation +does not rule out a restricted or checkpointed rounded-product method. + +## 4. Uniform drift bound + +Let `E_r(j)=|r_hat_j-r_j|` and `E_z(j)=|z_hat_j-z_j|`. +The first seed rotation and (1) give `E_r(0)<=8e`. Set `z_hat_0=1`. +As long as the states have modulus at most 2, the fixed `omega` rotation +gives + + E_r(j+1) <= E_r(j)+12e, + E_r(j) <= 12(j+1)e. (3) + +Use `chart(r_hat_j,Q)` for the update of `z_hat_j`. Section 2 gives + + E_z(j+1) <= E_z(j)+40E_r(j)+6e. + +Sum this inequality, starting at zero. For all `0<=j<=L`, + + E_z(j) <= (240j(j+1)+6j)e <= 256(j+1)^2 e. (4) + +The chosen `Q` makes the last bound at most `2^-(W+12)`. It also makes +the bound in (3) less than `1/64`. Hence each state is within `1/64` +of a unit phase and has modulus less than 2. Induction proves these +conditions and (3)--(4) together. There is no circular stability assumption. + +At emission use `chart(z_hat_j,W)`. Substitute +`E<=u/4096` and `delta=u/256` in (2). The tangent error is at most +`5u/2048= 2, 'working precision') + d.need(arrays and len(arrays[0]) > 0, 'nonempty arrays') + d.need(len({id(a) for a in arrays}) == len(arrays), 'distinct arrays') + length = len(arrays[0]) + d.need(all(len(a) == length for a in arrays), 'equal array lengths') + bound = F(0) + for data in arrays: + for i in range(length): + meter.access(data, i) + meter.input_reads += 1 + value = data[i] + d.need(isinstance(value, tuple) and len(value) == 2 + and all(type(v) in (int, F) for v in value), 'rational complex data') + d.need(d.rnd(value, bits) == value, 'fixed grid data') + bound = max(bound, abs(value[0])+abs(value[1])) + meter.initial_bound = bound + return bound + + +def _word_budget(bits, beta, length, stride, base_error, operator_norm, bound): + d.need(d.admitted(beta, bits, d.depth(length, stride)), 'inherited precision admission') + unit = F(1, 1 << bits) + d.need(3*base_error*unit <= 1, 'conjugation precision admission') + return 32*(1+operator_norm+base_error)*(bound+1)*unit + + +def _phase_descriptor(numerator, modulus): + d.need(type(numerator) is int and type(modulus) is int and modulus > 0, + 'integer phase descriptor') + + +def _coefficient(value): + return isinstance(value, tuple) and len(value) == 2 and all(type(v) in (int, F) for v in value) + + +def lift(x, y, b, beta, bits, stride, numerator, modulus, + sign=1, multiplier=d.ONE, shift=0): + """Add D times the pinned cyclic inverse times D inverse to a dirty target.""" + _phase_descriptor(numerator, modulus) + d.need(type(stride) is int and type(shift) is int + and type(sign) is int and sign in (-1, 1), 'integer stride and sign') + d.need(type(beta) in (int, F) and 0 < beta <= F(1, 2), 'stable margin') + d.need(_coefficient(b) and d.norm2(b) <= (1-beta)**2, 'stable coefficient') + d.need(_coefficient(multiplier), 'constant multiplier') + d.need(len(x) > 0 and shift % len(x) in (0, stride % len(x)), 'supported output shift') + meter = Meter() + meter.layout = {id(x): (0, 0), id(y): (1, 0)} + bound = _validate_arrays((x, y), bits, meter) + multiplier_bound = max(F(1), abs(multiplier[0])+abs(multiplier[1])) + base_error = 68*multiplier_bound/beta**5 + operator_norm = multiplier_bound/beta + meter.error_bound = _word_budget(bits, beta, len(x), stride, base_error, operator_norm, bound) + phase = Phase.prepare(numerator, modulus, bits, meter, len(x)) + phase_pass(x, phase, True) + phase_pass(y, phase, True) + meter.rewind() + d.lift(x, y, b, beta, bits, stride, sign, multiplier, shift, meter) + phase_pass(y, phase, False) + phase_pass(x, phase, False) + meter.rewind() + phase.clear() + return meter + + +def compose(x, y, factors, scratch, tau, beta, bits, stride, numerator, modulus, sign=1): + """Conjugate a complete PR985 factor word. Restore all borrowed arrays.""" + _phase_descriptor(numerator, modulus) + d.need(type(sign) is int and sign in (-1, 1) and type(stride) is int, 'integer stride and sign') + d.need(type(beta) in (int, F) and 0 < beta <= F(1, 2), 'stable margin') + d.need(factors and len(scratch) == (len(factors)-1).bit_length(), 'exact borrowed-array budget') + d.need(_coefficient(tau) and d.norm2(tau) == 1, 'unit gauge') + scale_bound = F(1) + for factor in factors: + d.need(isinstance(factor, tuple) and len(factor) == 3, 'factor descriptor') + inside, b, scale = factor + d.need(type(inside) is bool and _coefficient(b) and d.norm2(b) <= (1-beta)**2 + and type(scale) in (int, F), 'stable factor descriptor') + scale_bound = max(scale_bound, abs(scale)) + scratch = list(scratch) + arrays = [x, y]+scratch + meter = Meter() + if scratch: + meter.layout = {id(a): (0, i*len(x)) for i, a in enumerate(arrays[:-1])} + meter.layout[id(scratch[-1])] = (1, 0) + else: + meter.layout = {id(x): (0, 0), id(y): (1, 0)} + bound = _validate_arrays(arrays, bits, meter) + calls = d.leaf_calls(len(factors)) + growth = 2+scale_bound/beta + base_error = calls*68*(scale_bound+1)/beta**5*growth**(calls-1) + operator_norm = growth**calls + meter.error_bound = _word_budget(bits, beta, len(x), stride, base_error, operator_norm, bound) + phase = Phase.prepare(numerator, modulus, bits, meter, len(x)) + phase_pass(x, phase, True) + phase_pass(y, phase, True) + meter.rewind() + d.compose(x, y, factors, scratch, tau, beta, bits, stride, sign, meter) + phase_pass(y, phase, False) + phase_pass(x, phase, False) + meter.rewind() + phase.clear() + return meter + + +def working_bits(target_bits, incoming_bits, bound, beta, length, stride, + multiplier_bound=1, count=1, scale_bound=None): + """Choose a conservative precision for the complete conjugated word.""" + d.need(type(target_bits) is int and target_bits >= 0 + and type(incoming_bits) is int and incoming_bits >= 0 + and type(bound) in (int, F) and type(beta) in (int, F) + and bound >= 0 and 0 < beta <= F(1, 2) + and type(length) is int and length > 0 and type(stride) is int + and type(multiplier_bound) in (int, F) and multiplier_bound >= 1, 'working budget domain') + if scale_bound is None: + d.need(count == 1, 'leaf count') + base_error = 68*multiplier_bound/beta**5 + operator_norm = multiplier_bound/beta + else: + d.need(type(scale_bound) in (int, F) and scale_bound >= 1, 'factor scale bound') + calls = d.leaf_calls(count) + growth = 2+F(scale_bound)/beta + base_error = calls*68*(scale_bound+1)/beta**5*growth**(calls-1) + operator_norm = growth**calls + return max(2, incoming_bits, + target_bits+1+d.ceil_log2(32*(1+operator_norm+base_error)*(bound+1)), + d.ceil_log2(3*base_error), + d.ceil_log2(48*(d.depth(length, stride)+1)/beta**3)) diff --git a/research/dirty-stride-fast-phases/phase_stream.py b/research/dirty-stride-fast-phases/phase_stream.py new file mode 100644 index 0000000000..2db9eb3a53 --- /dev/null +++ b/research/dirty-stride-fast-phases/phase_stream.py @@ -0,0 +1,289 @@ +"""Generate quadratic phases with reversible rounded state updates.""" +from dataclasses import dataclass +from fractions import Fraction as F +import importlib.util +from pathlib import Path +import sys + + +_SPEC = importlib.util.spec_from_file_location( + 'fast_phase_predecessor', Path(__file__).resolve().parent/'vendor/pr1038/chirp.py') +old = importlib.util.module_from_spec(_SPEC) +sys.modules[_SPEC.name] = old +_SPEC.loader.exec_module(old) +d = old.d +real_round = old.real_round +rotate = old.rotate +quarter_turn = old.quarter_turn + + +@dataclass +class Meter(old.Meter): + seed_terms: int = 0 + seed_products: int = 0 + seed_divisions: int = 0 + seed_cleanups: int = 0 + seed_peak_bits: int = 0 + stream_charts: int = 0 + stream_steps: int = 0 + stream_reverse_steps: int = 0 + stream_cleanups: int = 0 + preparation_cleanups: int = 0 + state_bits: int = 0 + state_peak_l1: F = F(0) + + +def _width(value, meter): + value = F(value) + meter.seed_peak_bits = max(meter.seed_peak_bits, + abs(value.numerator).bit_length(), + value.denominator.bit_length()) + + +def _machin_add(total, denominator, weight, terms, meter, inverse=False): + """Add an exact arctangent sum and clear its power register.""" + power = denominator + if inverse: + for _ in range(terms): + power *= denominator*denominator + meter.seed_products += 1 + indices = range(terms-1, -1, -1) if inverse else range(terms) + for k in indices: + if inverse: + power //= denominator*denominator + meter.seed_divisions += 1 + term = F(weight if k % 2 == 0 else -weight, (2*k+1)*power) + total = total-term if inverse else total+term + meter.seed_terms += 1 + meter.seed_products += 1 + meter.seed_divisions += 1 + _width(total, meter) + _width(power, meter) + if not inverse: + power *= denominator*denominator + meter.seed_products += 1 + if not inverse: + for _ in range(terms): + power //= denominator*denominator + meter.seed_divisions += 1 + power -= denominator + d.need(power == 0, 'Machin power cleanup') + return total + + +def seed_coefficients(numerator, modulus, coefficient_bits, meter, + destination=(0, F(0), F(0)), inverse=False): + """Add or remove seed coefficients and clear the exact series state.""" + sign = -1 if inverse else 1 + residue = numerator % modulus + quadrant = (8*residue+modulus)//(2*modulus) + remainder = 4*residue-quadrant*modulus + if remainder == 0: + meter.seed_cleanups += 1 + return destination[0]+sign*(quadrant % 4), destination[1], destination[2] + precision = coefficient_bits+16 + pi_sum = _machin_add(F(0), 5, 16, precision, meter) + pi_sum = _machin_add(pi_sum, 239, -4, precision, meter) + pi_grid = real_round(pi_sum, precision) + angle = real_round(pi_grid*remainder/(2*modulus), precision) + meter.seed_products += 1 + meter.seed_divisions += 1 + d.need(abs(angle) < 1, 'seed reduced angle') + sine = sine_term = angle + cosine = cosine_term = F(1) + if angle: + squared = angle*angle + meter.seed_products += 1 + for k in range(1, precision): + sine_term = -sine_term*squared/(2*k*(2*k+1)) + cosine_term = -cosine_term*squared/((2*k-1)*2*k) + sine += sine_term + cosine += cosine_term + meter.seed_terms += 2 + meter.seed_products += 2 + meter.seed_divisions += 2 + for value in (sine, cosine, sine_term, cosine_term): + _width(value, meter) + tangent = real_round(sine/(1+cosine), coefficient_bits) + sine_output = real_round(sine, coefficient_bits) + output = (destination[0]+sign*(quadrant % 4), + destination[1]+sign*tangent, destination[2]+sign*sine_output) + meter.seed_divisions += 1 + if angle: + for k in range(precision-1, 0, -1): + sine -= sine_term + cosine -= cosine_term + sine_term = -sine_term*(2*k*(2*k+1))/squared + cosine_term = -cosine_term*((2*k-1)*2*k)/squared + meter.seed_terms += 2 + meter.seed_products += 2 + meter.seed_divisions += 2 + squared -= angle*angle + meter.seed_products += 1 + d.need(squared == 0, 'Taylor square cleanup') + sine -= angle + sine_term -= angle + cosine -= 1 + cosine_term -= 1 + d.need(sine == sine_term == cosine == cosine_term == 0, 'Taylor cleanup') + angle -= real_round(pi_grid*remainder/(2*modulus), precision) + pi_grid -= real_round(pi_sum, precision) + meter.seed_products += 1 + meter.seed_divisions += 1 + pi_sum = _machin_add(pi_sum, 239, -4, precision, meter, True) + pi_sum = _machin_add(pi_sum, 5, 16, precision, meter, True) + d.need(angle == pi_grid == pi_sum == 0, 'seed cleanup') + d.need(abs(tangent) <= F(1, 2) and abs(sine_output) <= 1, + 'bounded seed coefficients') + meter.seed_cleanups += 1 + return output + + +def chart(value, bits, meter, destination=(0, F(0), F(0)), inverse=False): + """Add or remove chart coefficients without changing their source.""" + x, y = value + if x >= abs(y): + quadrant = 0 + elif y >= abs(x): + quadrant = 1 + elif -x >= abs(y): + quadrant = 2 + else: + quadrant = 3 + cosine, sine = quarter_turn(value, (-quadrant) % 4) + d.need(cosine >= abs(sine), 'phase chart domain') + tangent = real_round(sine/(1+cosine), bits+8) + sine = real_round(sine, bits+8) + tangent = min(F(1, 2), max(-F(1, 2), tangent)) + sine = min(F(1), max(-F(1), sine)) + meter.phase_divisions += 1 + meter.stream_charts += 1 + sign = -1 if inverse else 1 + return (destination[0]+sign*quadrant, destination[1]+sign*tangent, + destination[2]+sign*sine) + + +@dataclass +class State: + phase: object + value: tuple + ratio: tuple + index: int = 0 + + def coefficients(self): + d.need(0 <= self.index < self.phase.length, 'phase stream index') + self.phase.meter.phase_evaluations += 1 + return chart(self.value, self.phase.bits, self.phase.meter) + + def erase_coefficients(self, coefficients): + result = chart(self.value, self.phase.bits, self.phase.meter, coefficients, True) + d.need(result == (0, 0, 0), 'emitted coefficient cleanup') + return result + + def step(self, inverse=False): + p = self.phase + if inverse: + d.need(self.index > 0, 'phase reverse index') + self.ratio = rotate(self.ratio, p.omega, p.state_bits, True, p.meter) + coefficients = chart(self.ratio, p.state_bits, p.meter) + self.value = rotate(self.value, coefficients, p.state_bits, True, p.meter) + coefficients = chart(self.ratio, p.state_bits, p.meter, coefficients, True) + self.index -= 1 + p.meter.stream_reverse_steps += 1 + else: + d.need(self.index < p.length, 'phase forward index') + coefficients = chart(self.ratio, p.state_bits, p.meter) + self.value = rotate(self.value, coefficients, p.state_bits, False, p.meter) + coefficients = chart(self.ratio, p.state_bits, p.meter, coefficients, True) + self.ratio = rotate(self.ratio, p.omega, p.state_bits, False, p.meter) + self.index += 1 + p.meter.stream_steps += 1 + d.need(coefficients == (0, 0, 0), 'state coefficient cleanup') + p.meter.state_peak_l1 = max(p.meter.state_peak_l1, + abs(self.value[0])+abs(self.value[1]), + abs(self.ratio[0])+abs(self.ratio[1])) + + def clear(self): + """Reverse the state word and clear its two initial values.""" + d.need(self.phase.busy, 'open phase stream') + while self.index: + self.step(True) + self.value = d.sub(self.value, d.ONE) + self.ratio = d.sub(self.ratio, self.phase.seed) + d.need(self.value == self.ratio == d.ZERO, 'phase state cleanup') + self.phase.busy = False + self.phase.meter.stream_cleanups += 1 + + +@dataclass +class Phase: + numerator: int + modulus: int + bits: int + length: int + state_bits: int + seed: tuple + omega: tuple + meter: Meter + active: bool = True + busy: bool = False + + @classmethod + def prepare(cls, numerator, modulus, bits, meter, length): + d.need(type(numerator) is int and type(modulus) is int and modulus > 0, + 'integer phase descriptor') + d.need(type(bits) is int and bits >= 2 and type(length) is int and length > 0, + 'phase stream domain') + state_bits = bits+2*(length+1).bit_length()+20 + meter.state_bits = state_bits + meter.coefficient_bits = bits+8 + meter.internal_bits = state_bits+24 + meter.descriptor_bits += abs(numerator).bit_length()+modulus.bit_length()+2 + a = numerator % modulus + first = seed_coefficients(a, modulus, state_bits+8, meter) + seed = rotate(d.ONE, first, state_bits, False, meter) + first = seed_coefficients(a, modulus, state_bits+8, meter, first, True) + d.need(first == (0, 0, 0), 'first seed coefficient cleanup') + omega = seed_coefficients(2*a, modulus, state_bits+8, meter) + return cls(a, modulus, bits, length, state_bits, seed, omega, meter) + + def start(self): + d.need(self.active, 'active phase preparation') + d.need(not self.busy, 'one open phase stream') + self.busy = True + return State(self, d.ONE, self.seed) + + def clear(self): + """Recompute the two seeds and subtract them from their registers.""" + d.need(self.active, 'active phase preparation') + d.need(not self.busy, 'phase stream must be clear') + first = seed_coefficients(self.numerator, self.modulus, self.state_bits+8, self.meter) + self.seed = rotate(self.seed, first, self.state_bits, True, self.meter) + self.seed = d.sub(self.seed, d.ONE) + first = seed_coefficients(self.numerator, self.modulus, self.state_bits+8, + self.meter, first, True) + d.need(first == (0, 0, 0), 'last seed coefficient cleanup') + self.omega = seed_coefficients(2*self.numerator, self.modulus, self.state_bits+8, + self.meter, self.omega, True) + d.need(self.seed == d.ZERO and self.omega == (0, 0, 0), 'phase preparation cleanup') + self.active = False + self.meter.preparation_cleanups += 1 + + +def phase_pass(data, phase, inverse): + """Apply one data pass and then reverse its phase state.""" + d.need(len(data) == phase.length, 'phase data length') + meter = phase.meter + before = meter.tape_moves + state = phase.start() + for i in range(len(data)): + coefficients = state.coefficients() + meter.access(data, i) + value = rotate(data[i], coefficients, phase.bits, inverse, meter) + meter.access(data, i, True) + data[i] = value + coefficients = state.erase_coefficients(coefficients) + state.step() + state.clear() + meter.phase_passes += 1 + meter.phase_record_moves += meter.tape_moves-before diff --git a/research/dirty-stride-fast-phases/test_package.py b/research/dirty-stride-fast-phases/test_package.py new file mode 100644 index 0000000000..cc646040e6 --- /dev/null +++ b/research/dirty-stride-fast-phases/test_package.py @@ -0,0 +1,96 @@ +"""Reject semantic faults and replay a relocated read-only package.""" +import difflib +import json +from pathlib import Path +import shutil +import subprocess +import sys +import tempfile + +HERE = Path(__file__).resolve().parent + + +def need(condition, message): + if not condition: + raise ValueError(message) + + +def probe(folder, expression): + program = ( + "import importlib.util,sys; from pathlib import Path; " + "p=Path(sys.argv[1]); s=importlib.util.spec_from_file_location('fault_check',p/'verify.py'); " + "v=importlib.util.module_from_spec(s); sys.modules[s.name]=v; s.loader.exec_module(v); " + +expression) + return subprocess.run([sys.executable, '-I', '-O', '-B', '-c', program, str(folder)], + text=True, capture_output=True, timeout=600) + + +def main(): + old = (HERE/'vendor/pr1038/chirp.py').read_text() + new = (HERE/'chirp.py').read_text() + difference = ''.join(difflib.unified_diff(old.splitlines(True), new.splitlines(True), + fromfile='PR1038/chirp.py', tofile='chirp.py', n=0)) + need(difference == (HERE/'DERIVATION.diff').read_text(), 'exact wrapper derivation') + cases = [ + ('wrong_second_difference', 'phase_stream.py', + 'omega = seed_coefficients(2*a, modulus, state_bits+8, meter)', + 'omega = seed_coefficients(a, modulus, state_bits+8, meter)', + 'v.phase_case(5,17,12,33)', 'ratio state drift'), + ('wrong_emitted_state', 'phase_stream.py', + 'return chart(self.value, self.phase.bits, self.phase.meter)', + 'return chart(self.ratio, self.phase.bits, self.phase.meter)', + 'v.phase_case(5,17,12,33)', 'phase quadrant'), + ('stale_reverse_chart', 'phase_stream.py', + 'self.value = rotate(self.value, coefficients, p.state_bits, True, p.meter)', + 'self.value = rotate(self.value, chart(p.seed, p.state_bits, p.meter), p.state_bits, True, p.meter)', + 'v.phase_case(5,17,12,33)', 'phase state cleanup'), + ('missing_guard_bits', 'phase_stream.py', + 'state_bits = bits+2*(length+1).bit_length()+20', 'state_bits = bits', + 'v.phase_case(1,31,2,65)', 'phase coefficient accuracy'), + ('uncleared_series', 'phase_stream.py', 'sine -= sine_term', 'sine -= 0', + 'v.phase_case(5,17,12,33)', 'Taylor cleanup'), + ('wrong_dirty_source_inverse', 'chirp.py', + 'phase_pass(x, phase, False)', 'phase_pass(x, phase, True)', + 'v.word_case(5,2,-3,10,2**40)', 'dirty source restoration'), + ('wrong_seed_erase', 'phase_stream.py', + 'self.seed = rotate(self.seed, first, self.state_bits, True, self.meter)', + 'self.seed = rotate(self.seed, first, self.state_bits, False, self.meter)', + 'v.phase_case(5,17,12,33)', 'phase preparation cleanup'), + ] + with tempfile.TemporaryDirectory(prefix='dirty-fast-phase-controls-') as temp: + root = Path(temp) + for name, filename, before, after, expression, expected in cases: + folder = root/name + shutil.copytree(HERE, folder) + path = folder/filename + source = path.read_text() + need(before in source, 'missing fault site: '+name) + path.write_text(source.replace(before, after, 1)) + result = probe(folder, expression) + need(result.returncode != 0 and 'ValueError: '+expected+'\n' in result.stderr, + 'wrong fault rejection: '+name+'\n'+result.stdout+result.stderr) + print('REJECT '+name+': '+expected, flush=True) + shutil.rmtree(folder) + folder = root/'source_corruption' + shutil.copytree(HERE, folder) + source = folder/'vendor/pr1038/oracle.py' + source.write_bytes(source.read_bytes()+b'\n') + result = probe(folder, 'v.pins()') + need(result.returncode != 0 and 'ValueError: source SHA256: vendor/pr1038/oracle.py\n' + in result.stderr, 'source corruption rejection') + print('REJECT source_corruption: source SHA256', flush=True) + shutil.rmtree(folder) + relocated = root/'relocated' + shutil.copytree(HERE, relocated) + for path in relocated.rglob('*'): + if path.is_file(): + path.chmod(0o444) + result = subprocess.run([sys.executable, '-I', '-O', '-B', str(relocated/'verify.py')], + cwd=root, capture_output=True, text=True, timeout=600) + need(result.returncode == 0, 'relocated receipt replay\n'+result.stdout+result.stderr) + print(result.stdout.strip(), flush=True) + print('PASS_FAST_PHASE_PACKAGE controls=8 relocation=1 derivation=1') + + +if __name__ == '__main__': + main() diff --git a/research/dirty-stride-fast-phases/vendor/pr1038/NOTICE.md b/research/dirty-stride-fast-phases/vendor/pr1038/NOTICE.md new file mode 100644 index 0000000000..ddb3b67bc9 --- /dev/null +++ b/research/dirty-stride-fast-phases/vendor/pr1038/NOTICE.md @@ -0,0 +1,38 @@ +# Attribution and license + +New material is prepared for Ali Romman (@aliromman), with substantial +OpenAI Codex assistance in research, derivation, implementation, checks, +and writing. One assigned AI agent prepared this package. Parent-agent +review is separate from independent human review. + +The new material uses Apache-2.0. Retained files keep their original bytes, +licenses, and notices. `SOURCE.json` records exact Git-blob and SHA256 pins. + +- Ali Romman's PR985 supplies the dirty arbitrary-stride solver, two-tape + transport, source-prefix closure, complete factor word, and inherited + conservative precision proof. The local code reads each validation value + once into a scalar. This makes the literal accesses agree with the meter. + The two edits are shown in `DERIVATION.diff`; the vendor copy is unchanged. +- Rohan Guliani's PR846 supplies the Euclidean decomposition. Charles + Bethin's PR764 supplies the shared-suffix estimate used in that lineage. +- Chafik Boukhalfa's PR479 supplies the retained chirp family and the exact + single-axis finite-difference rule. Its multi-axis stream and faster + scalar cost are not implemented or claimed by this package. +- Ross Richardson's PR948 checks PR479's state transitions by separate + symbolic and integer methods. It was read for scope and prior evidence; + no PR948 code is executed here. +- Chafik Boukhalfa's PR628 supplies the factor fixtures, dirty commutator, + and kernel certificates retained by PR985. The new radix-four check uses + its original root bytes and a separate exact polynomial/dense calculation. +- The retained manuscript supplies the chirp identity and transform context. + The PR621/PR628 kernel lineage credits OpenAI's exact-DFT preprint, + classical Newton and q-Pochhammer factorization, Alexey Kuznetsov, and + the methods listed in the original retained notices. + +Quadratic finite differences, the Machin identity, Taylor error estimates, +and three-shear plane rotations are established methods. No literature +priority is claimed. The contribution is their paid reversible phase +interface around the specific dirty solver and complete factor word. + +No new kappa, selected supplier, measured speedup, full transform, scalar +bit-reversible compiler, or full multiplication theorem is claimed. diff --git a/research/dirty-stride-fast-phases/vendor/pr1038/PROOF.md b/research/dirty-stride-fast-phases/vendor/pr1038/PROOF.md new file mode 100644 index 0000000000..fdbf277a03 --- /dev/null +++ b/research/dirty-stride-fast-phases/vendor/pr1038/PROOF.md @@ -0,0 +1,410 @@ +# Reversible quadratic phases for a dirty cyclic inverse + +## 1. Contract and the changed operator + +Let `Z^w x` have entry `x_((i-w) mod L)`. Let + + D_i = exp(2 pi i a i^2 / m), 0 <= i < L, + +where `a` is an integer and `m` is a positive integer. All phase residues +are reduced modulo `m`. No array of phase values is an input or a temporary. +The arrays are distinct arrays of complex dyadic numbers with grid spacing +`u=2^-W`. The source, target, and borrowed arrays can all contain arbitrary +nonzero data. The implementation measures an upper bound `B` on their +initial complex moduli by a paid scan of their component absolute sums. +Use the maximum complex modulus as the norm of an array, and use its +induced operator norm. + +The main leaf, with no output shift or multiplier, approximates + + y <- y + (I - b D Z^w D^-1)^-1 x, + +where `0 < beta <= 1/2` and `|b| <= 1-beta`. The source is restored exactly. +More generally, the optional constant multiplier `v` and allowed output +shift `s` give the increment + + D v Z^s (I-b Z^w)^-1 D^-1 x, s=0 or s=w modulo L. + +A negative call is the exact inverse instruction word. It approximates the +negative increment. It is not implemented by negating an approximate answer. +The error can depend on the initial dirty target. The precision bound below +includes that dependence. + +The complete factor interface implements `D F D^-1`, where `F` is the +complete PR985 factor word. This needs four phase passes for the whole word, +independent of its leaf count. The borrowed arrays keep the PR985 budget +`h=ceil(log2 n)` for `n` factors. This interface requires exactly `h` arrays. +It restores every borrowed array exactly. Only the source and target pass +through the phase maps; borrowed data can stay arbitrary. + +This changes the site coefficients. Put `U=D Z^w D^-1`. Its entry is + + (Ux)_i = exp(2 pi i a (i^2-j^2)/m) x_j, + j=(i-w) mod L. + +For the retained example `L=8, a=-3, m=16, w=3`, the phase residues on these +edges are `11, 9, 7, 5, 3, 1, 15, 13`. They are not one constant gauge. +The verifier separates its response from the unphased inverse by more than +`1/8` in a real component for a unit source and `b=1/2`. + +## 2. Relation to a retained transform + +PR479 and the retained manuscript use + + a_j = exp(-pi i sum_k s_k j_k^2/t_k). + +The new interface covers each single-axis family exactly: set `L=t`, +`a=-s`, and `m=2t`. For even `t`, this is the retained periodic chirp +`exp(-pi i s j^2/t)`. The sign, including its conjugate, is available by +changing the sign of the integer numerator. Arbitrary positive moduli and +integer numerators also give a well-defined diagonal on a finite array; +periodicity outside that array is not assumed. + +This is an executable dirty realization for that phase family. It is not a +complete transform. It does not implement the multi-axis mixed-radix stream, +arbitrary supplied phase arrays, variable coefficient magnitudes, or a +general variable-coefficient Euclidean solver. It does not replace PR479's +fast scalar interface or its asymptotic scalar row. The slower elementary +phase evaluator used here is fully charged in Section 7. + +## 3. Phase generation and coefficient accuracy + +Each array pass starts with `E=0` and `delta=a mod m`. After it uses index `j`, +it updates + + E <- (E+delta) mod m, + delta <- (delta+2a) mod m. + +Induction gives `E=a j^2 mod m` and `delta=a(2j+1) mod m` before each use. +This is the single-axis finite-difference rule from PR479. The implementation +restarts the rule on every pass. It never reads a stored coefficient table. + +Set `C=W+8`, `P=C+2 bit_length(C+32)+20`, and `e=2^-P`. At each index let + + q = floor(4E/m + 1/2), + h = 4E-q m, + theta = pi h/(2m). + +Then `|h| <= m/2`, `|theta| <= pi/4`, and +`exp(2 pi i E/m) = i^q exp(i theta)`. The code keeps `q modulo 4` after it +forms `h`. A zero `h` has exact zero shear coefficients. + +The remaining cases use the following explicit approximation. These are +error bounds for the executed routine, not assumptions supplied by its +admission guard. + +1. Use `pi = 16 atan(1/5) - 4 atan(1/239)`. Each arctangent uses its first + `P` alternating terms. Round each term separately to spacing `e`. Its + rounding error is at most `e/2`. Its omitted tail is less than `e` by + the alternating-series bound. Both Machin weights are included. Thus + the error in `pi_approx` is at most `(10P+20)e`. +2. Round `pi_approx*h/(2m)` to the same spacing. The resulting angle `X` + differs from `theta` by at most `(3P+6)e`. For the admitted `P` this + also proves `|X|<1`. +3. Start the sine term at `X` and the cosine term at `1`. For `k=1,...,P-1`, + multiply by `-X^2/((2k)(2k+1))` for sine and + `-X^2/((2k-1)(2k))` for cosine, rounding each new term to spacing `e`. + These recurrence factors have modulus at most `1/2`. If the preceding + term error is at most `e`, the next is at most `e/2+e/2=e`. Summing the + `P-1` term errors and the Taylor tail costs at most `Pe`. The tail is + less than `e` since `|X|<1` and `(2P)! > 2^P`. Both sine and cosine + are 1-Lipschitz on the real line. The total error against their values + at `theta` is therefore at most `E0=(4P+8)e`. +4. The exact cosine is at least `sqrt(2)/2`; `E0<1/4`. Thus the computed + `1+cosine` is greater than `1`, as is the exact denominator. Subtract + the two quotients to bound the error in `sine/(1+cosine)` by `2E0`. + The final internal division rounding adds at most `e/2`. The bound + `(8P+17)e` covers this result. +5. Round sine and tangent once more to spacing `2^-C`. The program checks + `(8P+17)2^-P <= 2^(-C-2)`. This holds for the stated `P`; it is also a + fail-closed integer/rational guard. To check all admitted widths, put + `b=bit_length(C+32)`. Then `C>=10`, `b>=6`, and `P<3*2^b`. Thus + `8P+17<2^(b+5)<=2^(2b+18)=2^(P-C-2)`. Each final coefficient is within + `3*2^-C/4 < 2^-C` of its exact value. In particular the rounded + coefficients satisfy `|tangent|<=1/2` and `|sine|<=1`. + +All polynomial sums are sums of dyadic values and are exact. No uncharged +trigonometric library call or supplied approximate root is used. The +independent oracle computes outward rational intervals for the full angle +`2 pi E/m`. It does not use the quadrant reduction, shears, recurrence of +residues, or production term rounding. + +## 4. An exactly reversible grid rotation + +For exact `t=tan(theta/2)` and `s=sin(theta)`, the three real shears + + x1 = x - t y, + y1 = y + s x1, + x2 = x1 - t y1 + +give the rotation by `theta`. This follows from +`1-st=cos(theta)` and `-2t+s t^2=-sin(theta)`. Append the exact quarter-turn +`i^q`. In the implementation each product in a shear is rounded to spacing +`u`, with ties toward positive infinity. Each subtraction or addition is +exact on the grid. + +The inverse first reverses the quarter-turn, then executes + + x1 = x2 + rnd(t y1), + y = y1 - rnd(s x1), + x = x1 + rnd(t y). + +Every rounded value in this inverse is the same rounded value that the +forward instruction used. Thus the inverse restores both components +exactly for every grid input, with no sign-symmetry assumption about the +rounding rule. Different array sites are disjoint, so both directions can +scan in increasing physical index order. + +Let `S_i` be this grid permutation and let `D_i` be the exact unit phase. +A shear with coefficient modulus at most `1` has Euclidean operator norm +at most `2`. The product of the other two shears has norm at most `4`. +The difference caused by the three coefficient errors is at most +`12*2^-C |z|`. The three data-rounding errors contribute at most +`(4+2+1)u/2`. Quarter-turns have norm `1`. The same reasoning applies to the +reverse word. Hence, in both directions, + + |S_i(z)-D_i z| <= 4(|z|+1)u, + |S_i^-1(z)-D_i^-1 z| <= 4(|z|+1)u. (1) + +The sharper coefficient contribution is not needed. The program uses the +same immutable descriptor, precision and deterministic evaluator in both +directions. No coefficient approximation is recovered from modified data. + +## 5. Conjugation, restoration, and the dirty-value error + +Execute, in this order: + + S^-1 on x; S^-1 on y; + the complete pinned dirty word A, with the requested sign; + S on y; S on x. + +The middle word restores its source and every borrowed array exactly. The +outer maps are exact grid inverses. Consequently the whole construction +restores the original source and all borrowed arrays exactly. Replacing +`A` by its complete inverse gives the exact inverse of this whole word. + +Suppose the ideal increment of the middle word is `F x`, with operator norm +at most `K`, and its execution error at initial array bound `D0` is at most + + H(D0+1)u. + +Both signs have the same bound. Require `W>=2`, the inherited PR985 admission, +and `3Hu<=1`. From (1), the transformed arrays have norm at most +`2(B+1)`; the unchanged borrowed arrays also satisfy this bound. Thus the +middle error is at most `3H(B+1)u <= B+1`. If `V` is the target just after +the middle word, then + + |V| <= (3+2K)(B+1). + +The first target transformation contributes at most `4(B+1)u` to final +error. The first source transformation contributes at most `4K(B+1)u`. +The middle word contributes at most `3H(B+1)u`. The final target rotation +contributes at most `(16+8K)(B+1)u`. Their sum is at most + + 32(1+K+H)(B+1)u. (2) + +This proves the stated error against the ideal phased operator, including +large dirty targets. It does not assume that the rounded increment is +independent of those targets. It does not tighten PR985's precision result. + +For a leaf, let `M>=max(1,|v|)`. The program uses the component absolute sum +of the constant multiplier as one such bound. From the pinned leaf theorem, +we can use + + H=68 M/beta^5, K=M/beta. + +For a complete word with `n` factors, let `r>=1` bound every factor scale. +The validated inner multiplier has modulus at most `r`; the outer multiplier +also has modulus at most `r` since `|b|<1`. Put + + N(1)=1, + N(n)=2N(floor(n/2))+2N(ceil(n/2)), + G=2+r/beta, + H=N*68(r+1)*G^(N-1)/beta^5, + K=G^N. + +PR985 supplies the middle execution bound with this `H`. Each ideal factor +has norm at most `r/beta`; `N>=n`, so `K` bounds the product. Exactly `N` +leaves and both signs of every dirty commutator are executed. The generic +kernel tests do not claim that arbitrary input factors are a certified +approximation to a particular retained polynomial. + +To pay execution error to `2^(-p-1)`, `working_bits` takes the maximum of +incoming grid precision, `2`, + + p+1+ceil(log2(32(1+K+H)(B+1))), + ceil(log2(3H)), + ceil(log2(48(depth(L,w)+1)/beta^3)). + +A retained factor certificate must separately pay its approximation error, +scaled by the source norm. Unit diagonal conjugation preserves that norm +bound. The radix-four test does this explicitly in Section 8. + +## 6. Data tapes, magnitude, and storage + +The local `dirty_stride.py` has two small changes from the pinned source. +Each input-validation loop reads `data[i]` once into a scalar before it +compares the rounded value. The predecessor expression read it twice but +charged one access. `DERIVATION.diff` records both changes. No arithmetic +or transport instruction in the inverse solver is changed. A separate array +observer checks literal reads and writes against the reported meter and +also gives a physical-movement lower bound. The original source stays in +`vendor/pr985/` with exact byte pins. + +Use the same two physical data tapes as PR985. A leaf puts one array on each +tape. A complete word puts its last borrowed array on tape B and its source, +target and other borrowed arrays consecutively on tape A. The inherited +routing swaps a leaf target onto B when necessary and restores that carrier. +No tape head is assigned separately to each array. + +Every array read and write, including phase scans and input validation, +charges the distance from the previous physical head position plus one. +The wrapper returns both heads to their origins before the middle word and +after the final phase passes. The middle word also pays its own returns. +The phase generator uses constant scalar registers beside the parked data +head. Reading or writing a complex record costs its complete record width. + +For a leaf, input validation costs at most `4L` record moves. The four phase +scans cost at most `16L` together. Their head returns cost at most `4L`. +Adding this to PR985's `4096L` bound gives the conservative bound `4160L`. +For a complete word, input validation, phase scans, travel between the +source and target slots, and returns cost at most `64(h+2)L`. Add this to +PR985's `8192(N+1)(h+2)L` bound. The result is at most + + 8256(N+1)(h+2)L. + +No control counter or coefficient tape is counted as a free data tape; their +bit costs are stated in Section 7. + +A leaf uses no borrowed data array. A word uses exactly `h` arbitrary +borrowed arrays. The implementation holds only a constant number of clean +arithmetic registers in phase preparation, rotation, or one leaf. A phase +object contains scalar descriptors and one scalar pi approximation. There +is no list or cache of arithmetic values indexed by `L` or `W`. The array +list in the word contains references to existing arrays only. The `n` +factor descriptors are immutable input records, with `O(nS)` input bits +for descriptor width `S`. Their sequential access is paid below. This +input storage is distinct from clean arithmetic workspace. + +The inherited leaf index stack occupies `O(log^2 L)` metadata bits. The word +also uses `O(log n (log n+log L))` metadata bits. The phase stream adds +`O(log m+log L+log P)` control bits and a constant number of `O(P)` scalar +values. Scalar arithmetic scratch is additional and is included below. +This is not a bit-reversible implementation of its scalar arithmetic. + +An explicit data-magnitude bound is + + 16384(M+1)(B+1)/beta^7 for a leaf, + 16384(r+1)(B+1)G^N/beta^7 for a word. + +PR985 bounds the middle values with transformed input bound `2(B+1)`. +The additional factor covers all three shear intermediates, whose norms +are at most eight times the input plus `4u`. The bound also covers the +untransformed initial arrays. The integer reserve is the ceiling of the +base-two logarithm of this bound, plus two sign/rounding bits. There is no +silent fixed-width overflow assumption. The finite replay measures array +writes and each transient real/imaginary pair inside the phase shears. It +does not measure the internal Machin integer widths; those use the +analytical bound in Section 7. + +## 7. Complete phase and bit cost + +Let `V=bit_length(m)`, and include the original numerator's full input width +when its descriptor is read and reduced. The phase stream uses `O(V)` bit +additions and reductions at each site. Its exact quadratic residue needs +no large square. Initial descriptor normalization and all precision-control +arithmetic are paid setup operations. + +One invocation prepares pi once, with `2P` rounded terms. It then evaluates +exactly `4L` phase coefficients. Each nontrivial evaluation has `2(P-1)` +Taylor terms. With the counters used by the reference, the phase layer has +at most + + 2P(4L+1) Taylor/Machin terms, + 4P(4L+1) counted scalar products, + 2P(4L+1) counted scalar divisions. + +The product total includes the three data shears per phase use. Integer +loop expressions and grid rounding are not extra unit-cost oracles: they +are included in the bit cost below. The scalar-operation count is +`O((L+1)P)`, not `O(L)`. The phase setup is not repeated at each inner leaf +of a complete factor word. + +The largest Machin power is below `239^(2P+1)`, so it has fewer than +`16P+8` bits. Every term is rounded before it is added to an accumulator. +The Taylor products and denominators have `O(P+V+log P)` bits. Include the +fractional precision, integer data reserve, original immutable coefficient +and phase-descriptor widths, and these phase temporaries in a scalar size +`T`. The exact rational budget constants `H` and `K` must also fit. If `S` +bounds the bit widths of each input scalar descriptor, their numerators and +denominators can have `O((N+1)S)` bits even when their magnitudes are small. +A constant multiple of `(N+1)S+P+V+2W+I+log L` suffices for `T`, where `I` +is the integer data reserve. This pays the rational control values as well +as the data arithmetic. It gives `O(T)` scalar storage, with a fixed number +of arithmetic work tapes. + +Let `C(T)` pay a scalar product or quotient, including rounding, on that +workspace. Schoolbook fixed-width arithmetic gives the sufficient choice +`C(T)=O(T^2)` on a fixed number of tapes. The exact-rational Python program +is a reference for the instructions and rounding. It is not a timing +implementation of that tape bound. + +For the complete `n`-factor word, a conservative paid bill is + + external factor production + + O((N L + (L+1)P + n+N+1) C(T)) + + O(N(h+2)L T + n N T) + + O((L+1)V + L log L). + +The API takes the factor descriptors as inputs; their production is an +inherited external cost and is not improved here. The `n+N+1` term pays +descriptor normalization, factor-margin validation, and the exact rational +precision-control computation, including binary powers and leaf counting. +The `n N T` term pays worst-case sequential factor-descriptor access for +every leaf. Validation of factor margins, source reads, all copies and +swaps, inverse calls, rounding, Euclidean children, and both head returns +are included. The elementary phase setup and term counters fit the displayed +phase cost. Oversized input descriptors and computing the precision bound +are included in setup; their lengths are not inferred from their magnitudes. +For a leaf use `N=1`, `h=0` and its constant descriptor cost. + +This construction does not inherit PR479's fast exponential complexity. +It supplies an explicit, slower complete phase interface. Fast phase +preparation can replace it only with a separate workspace and error proof. +No new supplier constant or kappa follows from this bill. + +## 8. Independent checks and a retained kernel + +The checker uses an outward rational interval oracle for full-angle sine +and cosine. It computes residues by direct integer squares. It evaluates +cyclic inverses by finite geometric sums. Thus it does not repeat the +production finite differences, shear rotation, or Euclidean recursion. +All error comparisons include the oracle interval width. + +For the retained `r=4, k=1` polynomial, set `zeta=i`. Its coefficients are +`(1/r) product_(s=j+1)^(r-1)(1-zeta^s)`, hence + + Q(X)=1+(1+i)X/2+(1+i)X^2/4+X^3/4. + +The unchanged PR628 root fixture, retained by PR985, gives a rational +factor polynomial `Qhat`. The checker expands its roots exactly and bounds +`epsilon=sum_j |Q_j-Qhat_j|` by component absolute sums. It checks each +stable factor margin. Therefore `A=4/beta^3` bounds `Qhat(U)^-1` for every +unitary weighted cyclic shift `U` in this family. If `A epsilon<1`, the +resolvent identity gives + + ||Q(U)^-1-Qhat(U)^-1|| <= A^2 epsilon/(1-A epsilon). + +The checker adds this bound times the source norm to (2). It separately +uses exact rational Gaussian elimination on the matrix for `Q(tau Z^w)`, +with interval right-hand sides after the ideal chirp. Two executions use +nonconstant quadratic chirps, dirty borrowed arrays, and the 128-bit total +budget. The length-eight case is in the retained periodic phase family. +Both the dense ideal comparison and exact complete-word inverse +are checked. This validates these cases and their stated uniform analytic +bound. It is not a complete-transform or multiplication verification. + +The finite replay is evidence for the written statements. It is not a Lean +proof, an independent human review, or a proof that the inherited PR985, +all-size compiler, common-label, physical supplier, analytic-transfer, and +multiplication assumptions hold. diff --git a/research/dirty-stride-fast-phases/vendor/pr1038/SOURCE.json b/research/dirty-stride-fast-phases/vendor/pr1038/SOURCE.json new file mode 100644 index 0000000000..3fc9cd2d0a --- /dev/null +++ b/research/dirty-stride-fast-phases/vendor/pr1038/SOURCE.json @@ -0,0 +1,123 @@ +{ + "base": "8bddbee377e2f3ce1f1247196c4b290d55767d20", + "sources": [ + { + "pr": 985, + "commit": "9df42a99dd8497b6e7d46cd2507bb89071372640", + "path": "research/dirty-stride-source-reuse/dirty_stride.py", + "local": "vendor/pr985/dirty_stride.py", + "sha256": "129b2bc8cf89159aae0a810fad4703caefefccf1c5f4805dd3c9b5e8d967354b", + "blob": "f469930b8b7d79cbc05a884dc2e224e32db7fc4a" + }, + { + "pr": 985, + "commit": "9df42a99dd8497b6e7d46cd2507bb89071372640", + "path": "research/dirty-stride-source-reuse/PROOF.md", + "local": "vendor/pr985/PROOF.md", + "sha256": "af6560b048b20649aebd31d68815948a104d2c87efcdf556c85a10f48bbc0489", + "blob": "01a3da73470dc773a33948ee5f8cfc8848798411" + }, + { + "pr": 985, + "commit": "9df42a99dd8497b6e7d46cd2507bb89071372640", + "path": "research/dirty-stride-source-reuse/NOTICE.md", + "local": "vendor/pr985/NOTICE.md", + "sha256": "759fd9ada0c70e657a0b903089b32c84239b05b1f214cb01e81e9ca61b0d87f4", + "blob": "a96d0c7fea342d1e6990b4d16d731038ba983511" + }, + { + "pr": 985, + "commit": "9df42a99dd8497b6e7d46cd2507bb89071372640", + "path": "research/dirty-stride-source-reuse/README.md", + "local": "vendor/pr985/README.md", + "sha256": "96527e5f166ab8223e5c652975f32b4ba1b6b1d7fd41667a46bf515571fdaf13", + "blob": "451ac2bac0cc263d03f3a489fc5dabd653e9a4c5" + }, + { + "pr": 985, + "commit": "9df42a99dd8497b6e7d46cd2507bb89071372640", + "path": "research/dirty-stride-source-reuse/LICENSE", + "local": "vendor/pr985/LICENSE", + "sha256": "c71d239df91726fc519c6eb72d318ec65820627232b2f796219e87dcf35d0ab4", + "blob": "261eeb9e9f8b2b4b0d119366dda99c6fd7d35c64" + }, + { + "pr": 479, + "commit": "40dd01058439cdffae3e3b41c4fe52edaccf5f91", + "path": "research/streamed-chirp-phases/README.md", + "local": "vendor/pr479/README.md", + "sha256": "ce0750869369141e7f1e70f8f5872a0b8d0a72b33bd0aad163798f64ee0ed442", + "blob": "b75ad5c75b6ab221e362771f2a788e54700806a8" + }, + { + "pr": 479, + "commit": "40dd01058439cdffae3e3b41c4fe52edaccf5f91", + "path": "research/streamed-chirp-phases/PROOF.md", + "local": "vendor/pr479/PROOF.md", + "sha256": "0c4438e09ec34ec8208dff90851fbc05bee0cbbefcac4d90208cfc63facc0207", + "blob": "fd05f84d603a509db594122c723d19cdf1cb029a" + }, + { + "pr": 479, + "commit": "40dd01058439cdffae3e3b41c4fe52edaccf5f91", + "path": "research/streamed-chirp-phases/SOURCE.json", + "local": "vendor/pr479/SOURCE.json", + "sha256": "a90653f2ab1a1ac32880e2b3818c13b8b7a6c3972b67f0f769ecb6641d8170ab", + "blob": "53657bfc2985f8cf4038347e5f577238ee2cac46" + }, + { + "pr": 628, + "retained_by_pr": 985, + "commit": "9df42a99dd8497b6e7d46cd2507bb89071372640", + "path": "research/dirty-stride-source-reuse/vendor/pr628/FACTORS.json", + "local": "vendor/pr628/FACTORS.json", + "sha256": "643d60ac24d5e67f6c9305175bd1c99bf5e42edfef6d1995c051e0b98ba446a7", + "blob": "d3894bd2a1660229f1a7ed16c9d982689effc96f" + }, + { + "pr": 628, + "retained_by_pr": 985, + "commit": "9df42a99dd8497b6e7d46cd2507bb89071372640", + "path": "research/dirty-stride-source-reuse/vendor/pr628/NOTICE.md", + "local": "vendor/pr628/NOTICE.md", + "sha256": "abeacf1ad1491d5797cfcf1f3832799bbf3b24b99597229d7bc884059a51ec68", + "blob": "22fbeb8ef0a28dd6b3feea718492091319ceb96e" + }, + { + "pr": 628, + "retained_by_pr": 985, + "commit": "9df42a99dd8497b6e7d46cd2507bb89071372640", + "path": "research/dirty-stride-source-reuse/vendor/pr628/SOURCE.json", + "local": "vendor/pr628/SOURCE.json", + "sha256": "38a425e2d3a72a5e39b18b54c728e40407e121d481a1c6423edd5c18e9b68a04", + "blob": "59bddaea4d66725bf675a67b6d2cf85a061b1fae" + }, + { + "pr": 628, + "retained_by_pr": 985, + "commit": "9df42a99dd8497b6e7d46cd2507bb89071372640", + "path": "research/dirty-stride-source-reuse/vendor/pr628/check.py", + "local": "vendor/pr628/check.py", + "sha256": "a3e46b556894d71a1c83bbbac01c50c1a541285fc9422a3567d5eb1d1ee86e7d", + "blob": "5ef6dd67b6c6cf949fedf16c0cd7b34ca4ccafc2" + }, + { + "pr": 628, + "retained_by_pr": 985, + "commit": "9df42a99dd8497b6e7d46cd2507bb89071372640", + "path": "research/dirty-stride-source-reuse/vendor/pr628/interval_kernel.py", + "local": "vendor/pr628/interval_kernel.py", + "sha256": "9b7bb24640a9a7d3f83c7951dfef483ec933db4c50de6c0130b85f450eeddfbe", + "blob": "1182175a55cb871481de0f56ea9ae56e474705fe" + }, + { + "pr": 628, + "retained_by_pr": 985, + "commit": "9df42a99dd8497b6e7d46cd2507bb89071372640", + "path": "research/dirty-stride-source-reuse/vendor/pr628/scan.py", + "local": "vendor/pr628/scan.py", + "sha256": "1a9c01ca1a29a9be1e756858161202621db29d5d332a4266b4017a9502a50d7e", + "blob": "d620e493596f12578c1ce47f216c44c79b65de08" + } + ] +} diff --git a/research/dirty-stride-fast-phases/vendor/pr1038/chirp.py b/research/dirty-stride-fast-phases/vendor/pr1038/chirp.py new file mode 100644 index 0000000000..8a5c397148 --- /dev/null +++ b/research/dirty-stride-fast-phases/vendor/pr1038/chirp.py @@ -0,0 +1,295 @@ +"""Use reversible quadratic phases around the pinned dirty-array solver.""" +from dataclasses import dataclass +from fractions import Fraction as F +import importlib.util +from pathlib import Path +import sys + +_SPEC = importlib.util.spec_from_file_location( + 'dirty_phase_predecessor', Path(__file__).resolve().parent/'dirty_stride.py') +d = importlib.util.module_from_spec(_SPEC) +sys.modules[_SPEC.name] = d +_SPEC.loader.exec_module(d) + + +def real_round(value, bits): + scale = 1 << bits + value = value*scale + F(1, 2) + return F(value.numerator//value.denominator, scale) + + +@dataclass +class Meter(d.Meter): + phase_evaluations: int = 0 + phase_passes: int = 0 + phase_terms: int = 0 + phase_products: int = 0 + phase_divisions: int = 0 + phase_record_moves: int = 0 + phase_recurrence_steps: int = 0 + input_reads: int = 0 + descriptor_bits: int = 0 + coefficient_bits: int = 0 + internal_bits: int = 0 + error_bound: F = F(0) + initial_bound: F = F(0) + phase_peak_l1: F = F(0) + + +def pi_approx(bits, meter): + """Use the Machin identity and round each arctangent term.""" + value = F(0) + for denominator, scale in ((5, 16), (239, -4)): + power = denominator + for k in range(bits): + term = real_round(F(1, (2*k+1)*power), bits) + value += scale*term if k % 2 == 0 else -scale*term + power *= denominator*denominator + meter.phase_terms += 1 + meter.phase_divisions += 1 + meter.phase_products += 2 + return value + + +@dataclass(frozen=True) +class Phase: + numerator: int + modulus: int + bits: int + internal: int + pi: F + meter: Meter + + @classmethod + def prepare(cls, numerator, modulus, bits, meter): + d.need(type(numerator) is int and type(modulus) is int and modulus > 0, + 'integer phase descriptor') + d.need(type(bits) is int and bits >= 2, 'phase precision') + coefficient_bits = bits+8 + internal = coefficient_bits + 2*(coefficient_bits+32).bit_length()+20 + d.need((8*internal+17)*F(1, 1 << internal) + <= F(1, 1 << (coefficient_bits+2)), 'phase coefficient budget') + meter.descriptor_bits += abs(numerator).bit_length()+modulus.bit_length()+2 + meter.coefficient_bits = coefficient_bits + meter.internal_bits = internal + return cls(numerator % modulus, modulus, bits, internal, + pi_approx(internal, meter), meter) + + def coefficients(self, residue): + """Return a quarter-turn and two dyadic shear coefficients.""" + d.need(type(residue) is int and 0 <= residue < self.modulus, 'phase residue') + self.meter.phase_evaluations += 1 + quadrant = (8*residue+self.modulus)//(2*self.modulus) + remainder = 4*residue-quadrant*self.modulus + if remainder == 0: + return quadrant % 4, F(0), F(0) + angle = real_round(self.pi*remainder/(2*self.modulus), self.internal) + self.meter.phase_products += 1 + self.meter.phase_divisions += 1 + d.need(abs(angle) <= 1, 'reduced angle') + sine = sine_term = angle + cosine = cosine_term = F(1) + for k in range(1, self.internal): + sine_term = real_round(-sine_term*angle*angle/((2*k)*(2*k+1)), self.internal) + cosine_term = real_round(-cosine_term*angle*angle/((2*k-1)*(2*k)), self.internal) + sine += sine_term + cosine += cosine_term + self.meter.phase_terms += 2 + self.meter.phase_products += 4 + self.meter.phase_divisions += 2 + tangent = real_round(sine/(1+cosine), self.internal) + self.meter.phase_divisions += 1 + tangent = real_round(tangent, self.bits+8) + sine = real_round(sine, self.bits+8) + d.need(abs(tangent) <= F(1, 2) and abs(sine) <= 1, 'bounded phase shears') + return quadrant % 4, tangent, sine + + +def quarter_turn(value, quadrant): + real, imag = value + if quadrant == 0: + return value + if quadrant == 1: + return -imag, real + if quadrant == 2: + return -real, -imag + return imag, -real + + +def observe_phase(value, meter): + if meter is not None: + meter.phase_peak_l1 = max(meter.phase_peak_l1, abs(value[0])+abs(value[1])) + + +def rotate(value, coefficients, bits, inverse=False, meter=None): + """Apply three grid shears. The inverse uses the reverse instruction order.""" + quadrant, tangent, sine = coefficients + observe_phase(value, meter) + if inverse: + real, imag = quarter_turn(value, (-quadrant) % 4) + real += real_round(tangent*imag, bits) + observe_phase((real, imag), meter) + imag -= real_round(sine*real, bits) + observe_phase((real, imag), meter) + real += real_round(tangent*imag, bits) + observe_phase((real, imag), meter) + result = real, imag + else: + real, imag = value + real -= real_round(tangent*imag, bits) + observe_phase((real, imag), meter) + imag += real_round(sine*real, bits) + observe_phase((real, imag), meter) + real -= real_round(tangent*imag, bits) + observe_phase((real, imag), meter) + result = quarter_turn((real, imag), quadrant) + if meter is not None: + meter.phase_products += 3 + return result + + +def phase_pass(data, phase, inverse): + """Scan one array. Use quadratic finite differences and no phase table.""" + meter = phase.meter + before = meter.tape_moves + residue = 0 + difference = phase.numerator + for i in range(len(data)): + coefficients = phase.coefficients(residue) + meter.access(data, i) + value = rotate(data[i], coefficients, phase.bits, inverse, meter) + meter.access(data, i, True) + data[i] = value + residue = (residue+difference) % phase.modulus + difference = (difference+2*phase.numerator) % phase.modulus + meter.phase_recurrence_steps += 1 + meter.phase_passes += 1 + meter.phase_record_moves += meter.tape_moves-before + + +def _validate_arrays(arrays, bits, meter): + d.need(type(bits) is int and bits >= 2, 'working precision') + d.need(arrays and len(arrays[0]) > 0, 'nonempty arrays') + d.need(len({id(a) for a in arrays}) == len(arrays), 'distinct arrays') + length = len(arrays[0]) + d.need(all(len(a) == length for a in arrays), 'equal array lengths') + bound = F(0) + for data in arrays: + for i in range(length): + meter.access(data, i) + meter.input_reads += 1 + value = data[i] + d.need(isinstance(value, tuple) and len(value) == 2 + and all(type(v) in (int, F) for v in value), 'rational complex data') + d.need(d.rnd(value, bits) == value, 'fixed grid data') + bound = max(bound, abs(value[0])+abs(value[1])) + meter.initial_bound = bound + return bound + + +def _word_budget(bits, beta, length, stride, base_error, operator_norm, bound): + d.need(d.admitted(beta, bits, d.depth(length, stride)), 'inherited precision admission') + unit = F(1, 1 << bits) + d.need(3*base_error*unit <= 1, 'conjugation precision admission') + return 32*(1+operator_norm+base_error)*(bound+1)*unit + + +def _phase_descriptor(numerator, modulus): + d.need(type(numerator) is int and type(modulus) is int and modulus > 0, + 'integer phase descriptor') + + +def _coefficient(value): + return isinstance(value, tuple) and len(value) == 2 and all(type(v) in (int, F) for v in value) + + +def lift(x, y, b, beta, bits, stride, numerator, modulus, + sign=1, multiplier=d.ONE, shift=0): + """Add D times the pinned cyclic inverse times D inverse to a dirty target.""" + _phase_descriptor(numerator, modulus) + d.need(type(stride) is int and type(shift) is int + and type(sign) is int and sign in (-1, 1), 'integer stride and sign') + d.need(type(beta) in (int, F) and 0 < beta <= F(1, 2), 'stable margin') + d.need(_coefficient(b) and d.norm2(b) <= (1-beta)**2, 'stable coefficient') + d.need(_coefficient(multiplier), 'constant multiplier') + d.need(len(x) > 0 and shift % len(x) in (0, stride % len(x)), 'supported output shift') + meter = Meter() + meter.layout = {id(x): (0, 0), id(y): (1, 0)} + bound = _validate_arrays((x, y), bits, meter) + multiplier_bound = max(F(1), abs(multiplier[0])+abs(multiplier[1])) + base_error = 68*multiplier_bound/beta**5 + operator_norm = multiplier_bound/beta + meter.error_bound = _word_budget(bits, beta, len(x), stride, base_error, operator_norm, bound) + phase = Phase.prepare(numerator, modulus, bits, meter) + phase_pass(x, phase, True) + phase_pass(y, phase, True) + meter.rewind() + d.lift(x, y, b, beta, bits, stride, sign, multiplier, shift, meter) + phase_pass(y, phase, False) + phase_pass(x, phase, False) + meter.rewind() + return meter + + +def compose(x, y, factors, scratch, tau, beta, bits, stride, numerator, modulus, sign=1): + """Conjugate a complete PR985 factor word. Restore all borrowed arrays.""" + _phase_descriptor(numerator, modulus) + d.need(type(sign) is int and sign in (-1, 1) and type(stride) is int, 'integer stride and sign') + d.need(type(beta) in (int, F) and 0 < beta <= F(1, 2), 'stable margin') + d.need(factors and len(scratch) == (len(factors)-1).bit_length(), 'exact borrowed-array budget') + d.need(_coefficient(tau) and d.norm2(tau) == 1, 'unit gauge') + scale_bound = F(1) + for factor in factors: + d.need(isinstance(factor, tuple) and len(factor) == 3, 'factor descriptor') + inside, b, scale = factor + d.need(type(inside) is bool and _coefficient(b) and d.norm2(b) <= (1-beta)**2 + and type(scale) in (int, F), 'stable factor descriptor') + scale_bound = max(scale_bound, abs(scale)) + scratch = list(scratch) + arrays = [x, y]+scratch + meter = Meter() + if scratch: + meter.layout = {id(a): (0, i*len(x)) for i, a in enumerate(arrays[:-1])} + meter.layout[id(scratch[-1])] = (1, 0) + else: + meter.layout = {id(x): (0, 0), id(y): (1, 0)} + bound = _validate_arrays(arrays, bits, meter) + calls = d.leaf_calls(len(factors)) + growth = 2+scale_bound/beta + base_error = calls*68*(scale_bound+1)/beta**5*growth**(calls-1) + operator_norm = growth**calls + meter.error_bound = _word_budget(bits, beta, len(x), stride, base_error, operator_norm, bound) + phase = Phase.prepare(numerator, modulus, bits, meter) + phase_pass(x, phase, True) + phase_pass(y, phase, True) + meter.rewind() + d.compose(x, y, factors, scratch, tau, beta, bits, stride, sign, meter) + phase_pass(y, phase, False) + phase_pass(x, phase, False) + meter.rewind() + return meter + + +def working_bits(target_bits, incoming_bits, bound, beta, length, stride, + multiplier_bound=1, count=1, scale_bound=None): + """Choose a conservative precision for the complete conjugated word.""" + d.need(type(target_bits) is int and target_bits >= 0 + and type(incoming_bits) is int and incoming_bits >= 0 + and type(bound) in (int, F) and type(beta) in (int, F) + and bound >= 0 and 0 < beta <= F(1, 2) + and type(length) is int and length > 0 and type(stride) is int + and type(multiplier_bound) in (int, F) and multiplier_bound >= 1, 'working budget domain') + if scale_bound is None: + d.need(count == 1, 'leaf count') + base_error = 68*multiplier_bound/beta**5 + operator_norm = multiplier_bound/beta + else: + d.need(type(scale_bound) in (int, F) and scale_bound >= 1, 'factor scale bound') + calls = d.leaf_calls(count) + growth = 2+F(scale_bound)/beta + base_error = calls*68*(scale_bound+1)/beta**5*growth**(calls-1) + operator_norm = growth**calls + return max(2, incoming_bits, + target_bits+1+d.ceil_log2(32*(1+operator_norm+base_error)*(bound+1)), + d.ceil_log2(3*base_error), + d.ceil_log2(48*(d.depth(length, stride)+1)/beta**3)) diff --git a/research/dirty-stride-fast-phases/vendor/pr1038/dirty_stride.py b/research/dirty-stride-fast-phases/vendor/pr1038/dirty_stride.py new file mode 100644 index 0000000000..027108871d --- /dev/null +++ b/research/dirty-stride-fast-phases/vendor/pr1038/dirty_stride.py @@ -0,0 +1,409 @@ +"""Add a cyclic inverse to a target and restore the source exactly. + +This is a fixed-grid arithmetic model with two data tapes. Every data access +pays its head movement. See PROOF.md for the complete bill. +""" +from dataclasses import dataclass, field +from fractions import Fraction as F +from math import gcd + +ZERO = (F(0), F(0)) +ONE = (F(1), F(0)) + + +def need(ok, message): + if not ok: + raise ValueError(message) + + +def add(a, b): + return (a[0] + b[0], a[1] + b[1]) + + +def sub(a, b): + return (a[0] - b[0], a[1] - b[1]) + + +def mul(a, b): + return (a[0]*b[0] - a[1]*b[1], a[0]*b[1] + a[1]*b[0]) + + +def norm2(a): + return a[0]*a[0] + a[1]*a[1] + + +def div(a, b): + d = norm2(b) + need(d != 0, 'zero divisor') + v = mul(a, (b[0], -b[1])) + return (v[0]/d, v[1]/d) + + +def rnd(a, bits): + """Use nearest dyadic rounding, with ties toward positive infinity.""" + scale = 1 << bits + def component(v): + t = v*scale + F(1, 2) + return F(t.numerator // t.denominator, scale) + return (component(a[0]), component(a[1])) + + +def trunc(a, bits): + """Truncate each component toward zero. This cannot increase the norm.""" + scale = 1 << bits + return tuple(F(int(v*scale), scale) for v in a) + + +@dataclass +class Meter: + reads: int = 0 + writes: int = 0 + tape_moves: int = 0 + products: int = 0 + power_products: int = 0 + divisions: int = 0 + nodes: int = 0 + leaves: int = 0 + validation_reads: int = 0 + max_depth: int = 0 + heads: dict = field(default_factory=dict) + layout: dict = field(default_factory=dict) + + def access(self, data, index, write=False): + tape, start = self.layout[id(data)] + position = start+index + self.tape_moves += abs(position - self.heads.get(tape, 0)) + 1 + self.heads[tape] = position + if write: + self.writes += 1 + else: + self.reads += 1 + + def rewind(self): + for tape, position in self.heads.items(): + self.tape_moves += position + self.heads[tape] = 0 + + +@dataclass(frozen=True) +class View: + data: list + start: int + step: int + length: int + meter: Meter + + def index(self, i): + need(0 <= i < self.length, 'view index') + return (self.start + self.step*i) % len(self.data) + + def get(self, i): + j = self.index(i) + self.meter.access(self.data, j) + return self.data[j] + + def put(self, i, value): + j = self.index(i) + self.meter.access(self.data, j, True) + self.data[j] = value + + def part(self, first, length): + need(0 <= first and first + length <= self.length, 'view interval') + return View(self.data, self.start + self.step*first, self.step, length, self.meter) + + def reverse(self): + return View(self.data, self.start + self.step*(self.length-1), + -self.step, self.length, self.meter) + + +def depth(length, stride): + count = 0 + while True: + stride %= length + stride = min(stride, length-stride) + if stride == 0: + return count + length, stride = stride, (-length) % stride + count += 1 + + +def admitted(beta, bits, levels): + return (type(bits) is int and bits >= 0 and 0 < beta <= F(1, 2) + and 48*(levels+1)*F(1, 1 << bits) <= beta**3) + + +def power(b, exponent, bits, meter): + """Use a left-to-right binary product and scalar working storage.""" + need(exponent >= 1, 'positive exponent') + result = b + for k in range(exponent.bit_length()-2, -1, -1): + result = trunc(mul(result, result), bits) + meter.power_products += 1 + if (exponent >> k) & 1: + result = trunc(mul(result, b), bits) + meter.power_products += 1 + return result + + +@dataclass(frozen=True) +class Geometry: + root: tuple + divisor: int + cycle: int + inverse: int + bits: int + meter: Meter + + def exponents(self, length, stride, orientation): + if stride == 0: + return (self.cycle, self.cycle) + def exponent(displacement): + need(displacement % self.divisor == 0, 'cycle lattice') + return ((displacement // self.divisor)*self.inverse) % self.cycle or self.cycle + return (exponent(orientation*stride), exponent(orientation*(stride-length))) + + def coefficients(self, length, stride, orientation): + a, b = self.exponents(length, stride, orientation) + left = power(self.root, a, self.bits, self.meter) + right = left if a == b else power(self.root, b, self.bits, self.meter) + return (left, right) + + +def emit(target, i, value, multiplier, sign, bits): + if multiplier != ONE: + target.meter.products += 1 + value = rnd(mul(multiplier, value), bits) + old = target.get(i) + target.put(i, add(old, value) if sign == 1 else sub(old, value)) + + +def swap(source, target): + for i in range(source.length): + left, right = source.get(i), target.get(i) + source.put(i, right) + target.put(i, left) + + +def causal(source, carrier, stride, geometry, undo=False): + """Apply the causal lifts with block swaps; restore the dirty carrier.""" + a, _ = geometry.coefficients(source.length, stride, source.step) + n = source.length + if undo: + swap(source, carrier) + for first in range(((n-1)//stride)*stride, -1, -stride): + size = min(stride, n-first) + swap(source.part(first, size), carrier.part(first, size)) + if first: + for i in range(first+size-1, first-1, -1): + value = rnd(mul(a, carrier.get(i-stride)), geometry.bits) + geometry.meter.products += 1 + source.put(i, sub(source.get(i), value)) + else: + for first in range(0, n, stride): + size = min(stride, n-first) + if first: + for i in range(first, first+size): + value = rnd(mul(a, carrier.get(i-stride)), geometry.bits) + geometry.meter.products += 1 + source.put(i, add(source.get(i), value)) + swap(source.part(first, size), carrier.part(first, size)) + swap(source, carrier) + + +def expand(source, target, stride, geometry, multiplier, sign): + """Conjugate the seed injection by exact inverse causal grid maps.""" + causal(target, source, stride, geometry, undo=True) + _, b = geometry.coefficients(source.length, stride, source.step) + coefficient = mul(multiplier, b) + geometry.meter.products += 1 + for i in range(stride): + value = rnd(mul(coefficient, source.get(i)), geometry.bits) + geometry.meter.products += 1 + emit(target, i, value, ONE, sign, geometry.bits) + causal(target, source, stride, geometry) + + +def closure(source, target, stride, geometry, sign, level): + """Move the child target to the other tape and restore the carrier.""" + borrowed = source.part(0, stride) + carrier = target.part(0, stride) + ends = source.part(source.length-stride, stride) + swap(borrowed, carrier) + solve(ends, carrier, (-source.length) % stride, geometry, sign=sign, level=level+1) + swap(borrowed, carrier) + + +def diagonal(source, target, geometry, multiplier, sign): + coefficient = power(geometry.root, geometry.cycle, geometry.bits, geometry.meter) + denominator = sub(ONE, coefficient) + for i in range(source.length): + value = rnd(div(source.get(i), denominator), geometry.bits) + geometry.meter.divisions += 1 + emit(target, i, value, multiplier, sign, geometry.bits) + + +def solve(source, target, stride, geometry, multiplier=ONE, sign=1, level=0): + """Use integer metadata at each level; scalar temporaries stay in helpers.""" + length = source.length + stride %= length + if 2*stride > length: + source, target, stride = source.reverse(), target.reverse(), length-stride + geometry.meter.nodes += 1 + geometry.meter.max_depth = max(geometry.meter.max_depth, level) + if stride == 0: + diagonal(source, target, geometry, multiplier, sign) + return + causal(source, target, stride, geometry) + if sign == 1: + for i in range(length): + emit(target, i, source.get(i), multiplier, 1, geometry.bits) + closure(source, target, stride, geometry, 1, level) + expand(source, target, stride, geometry, multiplier, 1) + closure(source, target, stride, geometry, -1, level) + expand(source, target, stride, geometry, multiplier, -1) + else: + expand(source, target, stride, geometry, multiplier, 1) + closure(source, target, stride, geometry, 1, level) + expand(source, target, stride, geometry, multiplier, -1) + closure(source, target, stride, geometry, -1, level) + for i in range(length): + emit(target, i, source.get(i), multiplier, -1, geometry.bits) + causal(source, target, stride, geometry, undo=True) + + +def lift(x, y, b, beta, bits, stride, sign=1, multiplier=ONE, shift=0, meter=None): + """Add multiplier*(I-b Z^stride)^-1 x at output indices i+shift.""" + need(type(stride) is int and type(shift) is int and type(sign) is int + and sign in (-1, 1), 'integer stride and sign') + need(x is not y and len(x) == len(y) and len(x) > 0, 'distinct equal arrays') + need(shift % len(x) in (0, stride % len(x)), 'output shift must be zero or the leaf stride') + need(isinstance(b, tuple) and len(b) == 2 and norm2(b) <= (1-beta)**2, 'stable constant coefficient') + need(isinstance(multiplier, tuple) and len(multiplier) == 2, 'constant multiplier') + need(admitted(beta, bits, depth(len(x), stride)), 'precision admission') + meter = meter if meter is not None else Meter() + if not meter.layout: + meter.layout = {id(x): (0, 0), id(y): (1, 0)} + need(id(x) in meter.layout and id(y) in meter.layout + and {meter.layout[id(x)][0], meter.layout[id(y)][0]} == {0, 1}, 'two distinct data tapes') + for data in (x, y): + for i in range(len(data)): + meter.access(data, i) + meter.validation_reads += 1 + value = data[i] + need(rnd(value, bits) == value, 'fixed grid data') + divisor = gcd(len(x), stride) + cycle = len(x)//divisor + inverse = pow((stride//divisor) % cycle, -1, cycle) if cycle > 1 else 0 + geometry = Geometry(b, divisor, cycle, inverse, bits, meter) + solve(View(x, 0, 1, len(x), meter), View(y, shift, 1, len(y), meter), + stride, geometry, multiplier, sign) + meter.leaves += 1 + meter.rewind() + return meter + + +def leaf_bound(bound, beta, bits, multiplier_bound=1): + return (64*multiplier_bound*(bound+1)/beta**5 + 3)*F(1, 1 << bits) + + +def leaf_calls(count): + need(type(count) is int and count >= 1, 'positive factor count') + if count == 1: + return 1 + return 2*leaf_calls(count//2) + 2*leaf_calls(count-count//2) + + +def compose(x, y, factors, scratch, tau, beta, bits, stride, sign=1, meter=None): + """Use PR628's balanced dirty commutator with the new arbitrary-stride leaf.""" + need(factors and len(scratch) >= (len(factors)-1).bit_length(), 'dirty depth') + arrays = [x, y] + scratch + need(len({id(a) for a in arrays}) == len(arrays), 'all arrays distinct') + need(len(x) > 0 and all(len(a) == len(x) for a in arrays), 'array lengths') + need(isinstance(tau, tuple) and norm2(tau) == 1, 'immutable unit gauge') + need(type(sign) is int and sign in (-1, 1) and type(stride) is int, 'integer stride and sign') + need(admitted(beta, bits, depth(len(x), stride)), 'precision admission') + for factor in factors: + need(isinstance(factor, tuple) and len(factor) == 3, 'factor descriptor') + inside, b, scale = factor + need(type(inside) is bool and isinstance(b, tuple) and len(b) == 2 + and norm2(b) <= (1-beta)**2 and isinstance(scale, (int, F)), 'stable factor descriptor') + meter = meter if meter is not None else Meter() + if scratch: + layout = {id(a): (0, i*len(x)) for i, a in enumerate(arrays[:-1])} + layout[id(scratch[-1])] = (1, 0) + else: + layout = {id(x): (0, 0), id(y): (1, 0)} + need(not meter.layout or meter.layout == layout, 'kernel tape layout') + meter.layout = layout + for data in arrays: + for i in range(len(data)): + meter.access(data, i) + meter.validation_reads += 1 + value = data[i] + need(rnd(value, bits) == value, 'fixed grid data') + def run(src, dst, lo, hi, level, direction): + if hi-lo == 1: + inside, b, scale = factors[lo] + phase = (tau[0], -tau[1]) if inside else tau + coefficient = mul(b, phase) + m = phase if inside else (-b[0], -b[1]) + m = (m[0]*scale, m[1]*scale) + if meter.layout[id(src)][0] == meter.layout[id(dst)][0]: + spare = scratch[-1] + destination = View(dst, 0, 1, len(dst), meter) + carrier = View(spare, 0, 1, len(spare), meter) + swap(destination, carrier) + lift(src, spare, coefficient, beta, bits, -stride if inside else stride, + direction, m, -stride if inside else 0, meter) + swap(destination, carrier) + else: + lift(src, dst, coefficient, beta, bits, -stride if inside else stride, + direction, m, -stride if inside else 0, meter) + return + mid = (lo+hi)//2 + z = scratch[level] + if direction == 1: + run(src, z, lo, mid, level+1, 1) + run(z, dst, mid, hi, level+1, 1) + run(src, z, lo, mid, level+1, -1) + run(z, dst, mid, hi, level+1, -1) + else: + run(z, dst, mid, hi, level+1, 1) + run(src, z, lo, mid, level+1, 1) + run(z, dst, mid, hi, level+1, -1) + run(src, z, lo, mid, level+1, -1) + run(x, y, 0, len(factors), 0, sign) + meter.rewind() + return meter + + +def ceil_log2(value): + value = F(value) + need(value > 0, 'positive logarithm') + k = value.numerator.bit_length() - value.denominator.bit_length() + power_of_two = F(2)**k + return k if value <= power_of_two else k+1 + + +def execution_bound(count, r, beta, initial_bound, bits): + need(r >= 1 and initial_bound >= 0 and 0 < beta <= F(1, 2) + and type(bits) is int and bits >= 0, 'execution budget domain') + calls = leaf_calls(count) + a = 68*(r+1)/beta**5 + growth = 2+F(r)/beta + need(a*F(1, 1 << bits) <= 1, 'execution precision') + return calls*a*(initial_bound+1)*F(1, 1 << bits)*growth**(calls-1) + + +def working_bits(count, r, beta, initial_bound, target_bits, incoming_bits, length, stride): + need(r >= 1 and initial_bound >= 0 and 0 < beta <= F(1, 2) + and type(target_bits) is int and target_bits >= 0 + and type(incoming_bits) is int and incoming_bits >= 0 + and type(length) is int and length > 0 and type(stride) is int, 'working budget domain') + calls = leaf_calls(count) + a = 68*(r+1)/beta**5 + growth = 2+F(r)/beta + levels = depth(length, stride) + return max(incoming_bits, target_bits+1+ceil_log2(calls*a*(initial_bound+1)) + +(calls-1)*ceil_log2(growth), ceil_log2(a), + ceil_log2(48*(levels+1)/beta**3)) diff --git a/research/dirty-stride-fast-phases/vendor/pr1038/oracle.py b/research/dirty-stride-fast-phases/vendor/pr1038/oracle.py new file mode 100644 index 0000000000..2851b219ad --- /dev/null +++ b/research/dirty-stride-fast-phases/vendor/pr1038/oracle.py @@ -0,0 +1,188 @@ +"""Enclose the ideal phases and operators with rational intervals. + +This checker does not use the production phase reduction, grid rotation, +finite differences, or Euclidean inverse program. +""" +from fractions import Fraction as F +from functools import lru_cache +from math import gcd + + +def outward(lo, hi, bits): + scale = 1 << bits + return F((lo*scale).__floor__(), scale), F((hi*scale).__ceil__(), scale) + + +def point(value): + return F(value), F(value) + + +def plus(a, b, bits): + return outward(a[0]+b[0], a[1]+b[1], bits) + + +def neg(a): + return -a[1], -a[0] + + +def times(a, b, bits): + values = (a[0]*b[0], a[0]*b[1], a[1]*b[0], a[1]*b[1]) + return outward(min(values), max(values), bits) + + +def scale(a, value, bits): + return times(a, point(value), bits) + + +@lru_cache(maxsize=32) +def pi_interval(bits): + def arctangent(denominator): + total = F(0) + power = denominator + k = 0 + while True: + total += (-1 if k % 2 else 1)*F(1, (2*k+1)*power) + k += 1 + power *= denominator*denominator + remainder = F(1, (2*k+1)*power) + if remainder < F(1, 1 << (bits+8)): + other = total+(-1 if k % 2 else 1)*remainder + return min(total, other), max(total, other) + a, b = arctangent(5), arctangent(239) + return outward(16*a[0]-4*b[1], 16*a[1]-4*b[0], bits) + + +@lru_cache(maxsize=4096) +def phase(numerator, modulus, index, bits): + residue = numerator*index*index % modulus + if residue == 0: + return point(1), point(0) + angle = scale(pi_interval(bits), F(2*residue, modulus), bits) + squared = times(angle, angle, bits) + sine = sine_term = angle + cosine = cosine_term = point(1) + tail = F(7*7, 2) + k = 1 + unit = F(1, 1 << bits) + while True: + sine_term = scale(times(sine_term, squared, bits), -F(1, 2*k*(2*k+1)), bits) + cosine_term = scale(times(cosine_term, squared, bits), -F(1, (2*k-1)*2*k), bits) + sine = plus(sine, sine_term, bits) + cosine = plus(cosine, cosine_term, bits) + tail *= F(49, (2*k+1)*(2*k+2)) + k += 1 + if tail < unit and k >= 4: + sine_tail = tail*F(7, 2*k+1) + return (outward(cosine[0]-tail, cosine[1]+tail, bits), + outward(sine[0]-sine_tail, sine[1]+sine_tail, bits)) + + +def complex_point(value): + return point(value[0]), point(value[1]) + + +def complex_add(a, b, bits): + return plus(a[0], b[0], bits), plus(a[1], b[1], bits) + + +def complex_mul(a, b, bits): + return (plus(times(a[0], b[0], bits), neg(times(a[1], b[1], bits)), bits), + plus(times(a[0], b[1], bits), times(a[1], b[0], bits), bits)) + + +def rational_mul(a, b): + return a[0]*b[0]-a[1]*b[1], a[0]*b[1]+a[1]*b[0] + + +def inverse(value): + norm = value[0]*value[0]+value[1]*value[1] + return value[0]/norm, -value[1]/norm + + +def cyclic(source, coefficient, stride, bits): + length = len(source) + cycle = length//gcd(length, stride) + power = (F(1), F(0)) + for _ in range(cycle): + power = rational_mul(power, coefficient) + denominator_inverse = complex_point(inverse((1-power[0], -power[1]))) + result = [] + for i in range(length): + total = complex_point((0, 0)) + power = (F(1), F(0)) + for j in range(cycle): + total = complex_add(total, complex_mul(complex_point(power), source[(i-j*stride) % length], bits), bits) + power = rational_mul(power, coefficient) + result.append(complex_mul(total, denominator_inverse, bits)) + return result + + +def conjugate_input(source, numerator, modulus, bits): + return [complex_mul((phase(numerator, modulus, i, bits)[0], + neg(phase(numerator, modulus, i, bits)[1])), complex_point(value), bits) + for i, value in enumerate(source)] + + +def conjugate_output(source, numerator, modulus, bits): + return [complex_mul(phase(numerator, modulus, i, bits), value, bits) + for i, value in enumerate(source)] + + +def leaf(source, coefficient, stride, numerator, modulus, bits, multiplier=(F(1), F(0)), shift=0): + transformed = conjugate_input(source, numerator, modulus, bits) + solved = cyclic(transformed, coefficient, stride, bits) + shifted = [complex_mul(complex_point(multiplier), solved[(i-shift) % len(source)], bits) + for i in range(len(source))] + return conjugate_output(shifted, numerator, modulus, bits) + + +def kernel(source, factors, tau, stride, numerator, modulus, bits): + out = conjugate_input(source, numerator, modulus, bits) + for inside, coefficient, factor_scale in factors: + if inside: + gauge = tau[0], -tau[1] + b = rational_mul(coefficient, gauge) + multiplier = gauge[0]*factor_scale, gauge[1]*factor_scale + solved = cyclic(out, b, -stride, bits) + out = [complex_mul(complex_point(multiplier), solved[(i+stride) % len(source)], bits) + for i in range(len(source))] + else: + b = rational_mul(coefficient, tau) + multiplier = -coefficient[0]*factor_scale, -coefficient[1]*factor_scale + out = [complex_mul(complex_point(multiplier), value, bits) + for value in cyclic(out, b, stride, bits)] + return conjugate_output(out, numerator, modulus, bits) + + +def squared_distance_bound(actual, expected): + return sum(max(abs(actual[i]-expected[i][0]), abs(actual[i]-expected[i][1]))**2 for i in range(2)) + + +def polynomial(source, coefficients, stride, numerator, modulus, bits, tau): + """Use exact Gaussian elimination for the polynomial matrix.""" + length = len(source) + matrix = [[(F(0), F(0)) for _ in range(length)] for _ in range(length)] + power = F(1), F(0) + for j, coefficient in enumerate(coefficients): + weight = rational_mul(coefficient, power) + for i in range(length): + column = (i-stride*j) % length + old = matrix[i][column] + matrix[i][column] = old[0]+weight[0], old[1]+weight[1] + power = rational_mul(power, tau) + rhs = conjugate_input(source, numerator, modulus, bits) + for k in range(length): + pivot = next(i for i in range(k, length) if matrix[i][k] != (0, 0)) + matrix[k], matrix[pivot] = matrix[pivot], matrix[k] + rhs[k], rhs[pivot] = rhs[pivot], rhs[k] + divisor = inverse(matrix[k][k]) + matrix[k] = [rational_mul(divisor, v) for v in matrix[k]] + rhs[k] = complex_mul(complex_point(divisor), rhs[k], bits) + for i in range(length): + if i == k: + continue + weight = matrix[i][k] + subtrahend = [rational_mul(weight, v) for v in matrix[k]] + matrix[i] = [(v[0]-w[0], v[1]-w[1]) for v, w in zip(matrix[i], subtrahend)] + rhs[i] = complex_add(rhs[i], complex_mul(complex_point((-weight[0], -weight[1])), rhs[k], bits), bits) + return conjugate_output(rhs, numerator, modulus, bits) diff --git a/research/dirty-stride-fast-phases/vendor/pr1038/vendor/pr985/LICENSE b/research/dirty-stride-fast-phases/vendor/pr1038/vendor/pr985/LICENSE new file mode 100644 index 0000000000..261eeb9e9f --- /dev/null +++ b/research/dirty-stride-fast-phases/vendor/pr1038/vendor/pr985/LICENSE @@ -0,0 +1,201 @@ + Apache License + Version 2.0, January 2004 + http://www.apache.org/licenses/ + + TERMS AND CONDITIONS FOR USE, REPRODUCTION, AND DISTRIBUTION + + 1. Definitions. + + "License" shall mean the terms and conditions for use, reproduction, + and distribution as defined by Sections 1 through 9 of this document. + + "Licensor" shall mean the copyright owner or entity authorized by + the copyright owner that is granting the License. + + "Legal Entity" shall mean the union of the acting entity and all + other entities that control, are controlled by, or are under common + control with that entity. For the purposes of this definition, + "control" means (i) the power, direct or indirect, to cause the + direction or management of such entity, whether by contract or + otherwise, or (ii) ownership of fifty percent (50%) or more of the + outstanding shares, or (iii) beneficial ownership of such entity. + + "You" (or "Your") shall mean an individual or Legal Entity + exercising permissions granted by this License. + + "Source" form shall mean the preferred form for making modifications, + including but not limited to software source code, documentation + source, and configuration files. + + "Object" form shall mean any form resulting from mechanical + transformation or translation of a Source form, including but + not limited to compiled object code, generated documentation, + and conversions to other media types. + + "Work" shall mean the work of authorship, whether in Source or + Object form, made available under the License, as indicated by a + copyright notice that is included in or attached to the work + (an example is provided in the Appendix below). + + "Derivative Works" shall mean any work, whether in Source or Object + form, that is based on (or derived from) the Work and for which the + editorial revisions, annotations, elaborations, or other modifications + represent, as a whole, an original work of authorship. For the purposes + of this License, Derivative Works shall not include works that remain + separable from, or merely link (or bind by name) to the interfaces of, + the Work and Derivative Works thereof. + + "Contribution" shall mean any work of authorship, including + the original version of the Work and any modifications or additions + to that Work or Derivative Works thereof, that is intentionally + submitted to Licensor for inclusion in the Work by the copyright owner + or by an individual or Legal Entity authorized to submit on behalf of + the copyright owner. For the purposes of this definition, "submitted" + means any form of electronic, verbal, or written communication sent + to the Licensor or its representatives, including but not limited to + communication on electronic mailing lists, source code control systems, + and issue tracking systems that are managed by, or on behalf of, the + Licensor for the purpose of discussing and improving the Work, but + excluding communication that is conspicuously marked or otherwise + designated in writing by the copyright owner as "Not a Contribution." + + "Contributor" shall mean Licensor and any individual or Legal Entity + on behalf of whom a Contribution has been received by Licensor and + subsequently incorporated within the Work. + + 2. Grant of Copyright License. Subject to the terms and conditions of + this License, each Contributor hereby grants to You a perpetual, + worldwide, non-exclusive, no-charge, royalty-free, irrevocable + copyright license to reproduce, prepare Derivative Works of, + publicly display, publicly perform, sublicense, and distribute the + Work and such Derivative Works in Source or Object form. + + 3. Grant of Patent License. Subject to the terms and conditions of + this License, each Contributor hereby grants to You a perpetual, + worldwide, non-exclusive, no-charge, royalty-free, irrevocable + (except as stated in this section) patent license to make, have made, + use, offer to sell, sell, import, and otherwise transfer the Work, + where such license applies only to those patent claims licensable + by such Contributor that are necessarily infringed by their + Contribution(s) alone or by combination of their Contribution(s) + with the Work to which such Contribution(s) was submitted. If You + institute patent litigation against any entity (including a + cross-claim or counterclaim in a lawsuit) alleging that the Work + or a Contribution incorporated within the Work constitutes direct + or contributory patent infringement, then any patent licenses + granted to You under this License for that Work shall terminate + as of the date such litigation is filed. + + 4. Redistribution. You may reproduce and distribute copies of the + Work or Derivative Works thereof in any medium, with or without + modifications, and in Source or Object form, provided that You + meet the following conditions: + + (a) You must give any other recipients of the Work or + Derivative Works a copy of this License; and + + (b) You must cause any modified files to carry prominent notices + stating that You changed the files; and + + (c) You must retain, in the Source form of any Derivative Works + that You distribute, all copyright, patent, trademark, and + attribution notices from the Source form of the Work, + excluding those notices that do not pertain to any part of + the Derivative Works; and + + (d) If the Work includes a "NOTICE" text file as part of its + distribution, then any Derivative Works that You distribute must + include a readable copy of the attribution notices contained + within such NOTICE file, excluding those notices that do not + pertain to any part of the Derivative Works, in at least one + of the following places: within a NOTICE text file distributed + as part of the Derivative Works; within the Source form or + documentation, if provided along with the Derivative Works; or, + within a display generated by the Derivative Works, if and + wherever such third-party notices normally appear. The contents + of the NOTICE file are for informational purposes only and + do not modify the License. You may add Your own attribution + notices within Derivative Works that You distribute, alongside + or as an addendum to the NOTICE text from the Work, provided + that such additional attribution notices cannot be construed + as modifying the License. + + You may add Your own copyright statement to Your modifications and + may provide additional or different license terms and conditions + for use, reproduction, or distribution of Your modifications, or + for any such Derivative Works as a whole, provided Your use, + reproduction, and distribution of the Work otherwise complies with + the conditions stated in this License. + + 5. Submission of Contributions. Unless You explicitly state otherwise, + any Contribution intentionally submitted for inclusion in the Work + by You to the Licensor shall be under the terms and conditions of + this License, without any additional terms or conditions. + Notwithstanding the above, nothing herein shall supersede or modify + the terms of any separate license agreement you may have executed + with Licensor regarding such Contributions. + + 6. Trademarks. This License does not grant permission to use the trade + names, trademarks, service marks, or product names of the Licensor, + except as required for reasonable and customary use in describing the + origin of the Work and reproducing the content of the NOTICE file. + + 7. Disclaimer of Warranty. Unless required by applicable law or + agreed to in writing, Licensor provides the Work (and each + Contributor provides its Contributions) on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or + implied, including, without limitation, any warranties or conditions + of TITLE, NON-INFRINGEMENT, MERCHANTABILITY, or FITNESS FOR A + PARTICULAR PURPOSE. You are solely responsible for determining the + appropriateness of using or redistributing the Work and assume any + risks associated with Your exercise of permissions under this License. + + 8. Limitation of Liability. In no event and under no legal theory, + whether in tort (including negligence), contract, or otherwise, + unless required by applicable law (such as deliberate and grossly + negligent acts) or agreed to in writing, shall any Contributor be + liable to You for damages, including any direct, indirect, special, + incidental, or consequential damages of any character arising as a + result of this License or out of the use or inability to use the + Work (including but not limited to damages for loss of goodwill, + work stoppage, computer failure or malfunction, or any and all + other commercial damages or losses), even if such Contributor + has been advised of the possibility of such damages. + + 9. Accepting Warranty or Additional Liability. While redistributing + the Work or Derivative Works thereof, You may choose to offer, + and charge a fee for, acceptance of support, warranty, indemnity, + or other liability obligations and/or rights consistent with this + License. However, in accepting such obligations, You may act only + on Your own behalf and on Your sole responsibility, not on behalf + of any other Contributor, and only if You agree to indemnify, + defend, and hold each Contributor harmless for any liability + incurred by, or claims asserted against, such Contributor by reason + of your accepting any such warranty or additional liability. + + END OF TERMS AND CONDITIONS + + APPENDIX: How to apply the Apache License to your work. + + To apply the Apache License to your work, attach the following + boilerplate notice, with the fields enclosed by brackets "[]" + replaced with your own identifying information. (Don't include + the brackets!) The text should be enclosed in the appropriate + comment syntax for the file format. We also recommend that a + file or class name and description of purpose be included on the + same "printed page" as the copyright notice for easier + identification within third-party archives. + + Copyright [yyyy] [name of copyright owner] + + Licensed under the Apache License, Version 2.0 (the "License"); + you may not use this file except in compliance with the License. + You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + + Unless required by applicable law or agreed to in writing, software + distributed under the License is distributed on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + See the License for the specific language governing permissions and + limitations under the License. diff --git a/research/dirty-stride-fast-phases/vendor/pr1038/vendor/pr985/NOTICE.md b/research/dirty-stride-fast-phases/vendor/pr1038/vendor/pr985/NOTICE.md new file mode 100644 index 0000000000..a96d0c7fea --- /dev/null +++ b/research/dirty-stride-fast-phases/vendor/pr1038/vendor/pr985/NOTICE.md @@ -0,0 +1,39 @@ +# Attribution and license + +New material is prepared for Ali Romman (@aliromman), with substantial +OpenAI Codex assistance in the derivation, implementation, tests, and text. +One AI agent performed this work. The separate exact arithmetic paths are +not independent human review or a formal proof of the full theorem. + +The new files use Apache-2.0. The retained predecessor files keep their +original notices and source pins in `SOURCE.json`. + +- Rohan Guliani (@Rohansguliani), PR846, supplies the Euclidean chain and + closure decomposition and its linear-tape clean-storage solver. This + submission uses that decomposition and supplies a separate dirty transport + proof. It does not reuse the clean queue or its traffic constant. +- Chafik Boukhalfa (@chafreaky), PR628, supplies certified factor data, + paid coefficient and phase preparation, rounded integer lifts, and the + balanced dirty commutator. Its unchanged checker and dense solver check + the original kernels here. Its original dirty proof supplies the + conservative full-word perturbation argument. +- Charles Bethin (@cbethin), PR764, supplies the shared-suffix estimate used + in the weighted error argument, also used by PR846. No PR764 code is used. +- Ali Romman, PR621, supplies the uniform cyclic inverse bound and retained + kernel lineage. PR628's notice credits the pinned OpenAI exact-DFT + preprint, classical Newton and q-Pochhammer factorization, Alexey + Kuznetsov, and the cited preparation and filtering methods. Those credits + remain applicable. No priority claim is made for those methods. +- 0xkilgore's PR663 independently audits exact restoration and the binary + leaf count. PR673, PR689, PR730, PR744, and later echo work study other + tree and scratch contracts. The present source reuse acts inside a leaf. + +The overlap check also read Chafik Boukhalfa's PR862, PR905, PR929, and PR951, +Gareth Davidson's PR877, Wentao Li's PR805, and Chafik Boukhalfa's new +PR977. PR977 improves the precision of the clean solver. Ross Richardson's +PR982 borrows outputs for a different fanout interface. These results supply different +precision or storage results. They are not implemented, benchmarked, or +combined by this package. Their inspected heads are recorded in CLAIM.json. + +No literature-priority, new kappa, selected-supplier, production-speed, +bit-reversible compiler, or complete multiplication theorem claim is made. diff --git a/research/dirty-stride-fast-phases/vendor/pr1038/vendor/pr985/PROOF.md b/research/dirty-stride-fast-phases/vendor/pr1038/vendor/pr985/PROOF.md new file mode 100644 index 0000000000..01a3da7347 --- /dev/null +++ b/research/dirty-stride-fast-phases/vendor/pr1038/vendor/pr985/PROOF.md @@ -0,0 +1,428 @@ +# A dirty source for the arbitrary-stride closure + +## Contract and limits + +Let `Z_L^w x` have entry `x_((i-w) mod L)`. Let `0 < beta <= 1/2`, +`|b| <= 1-beta`, and `u=2^-W`. The source `x` and target `y` are distinct +length-`L` arrays on the complex dyadic grid of spacing `u`. They can contain +arbitrary data. The coefficient `b` is one immutable rational complex number. Reduce the +stride modulo `L`. The bounds below use a residue with `O(log L)` bits; +reading an oversized stride descriptor is additional input work. + +The function `lift` adds an approximation to + + m (I-b Z_L^w)^-1 x + +to a cyclic shift of `y`. It restores `x` exactly. It needs no other data +array, no reserved zero bits, and only a constant number of clean complex +arithmetic registers. The recursion keeps integer indices and array views. +It does not keep a clean array of coefficients or closure values. + +Reversible block swaps implement the data operations on two sequential +tapes. No blank data array or coefficient table is used. The source and +target occupy one tape each during a leaf. The complete kernel routes its +borrowed arrays over the same two tapes, with all movement charged below. +The optional output shift is zero or the leaf stride; these are the two +cases required by the kernel factors. The construction does not implement +the variable-coefficient case of PR846's Theorem 1. + +Let `K` be the number of non-base Euclidean levels. Assume + + 48 (K+1) u <= beta^3. + +If `|x_i| <= B`, the unscaled target increment differs from the inverse by +at most + + 64 (B+1) u / beta^5. + +A multiplier of modulus at most `M`, with `M>=1`, changes this bound to + + (64 M (B+1) / beta^5 + 3) u. + +The rounded increment can depend on the target's initial value within this +error bound. A negative call executes the exact inverse word. It restores +the completed positive call, including all source and target data. + +## Source reuse + +As in Rohan Guliani's PR846, normalize the stride and reverse the index views +when necessary. Then `0 < w <= L/2`. Split the source into chains with +indices `c, c+w, c+2w, ...`. Apply the causal integer lifts + + x_i += rnd(a_i x_(i-w)), i=w,...,L-1. + +Call the resulting source `z`. Each lift has an exact inverse on the grid: +subtract the same rounded product in reverse order. Add `z` to the target +before using any source entries as temporary storage. + +The last `w` source entries contain the closure right-hand side. The first +`w` source entries are disjoint from that tail. Use this prefix as a dirty +target for the recursive closure. It contains arbitrary values; it is not +cleared. Let `F` be the completed rounded closure call. Let `H` be the +linear chain expansion. A reversible rounded conjugation implements the +additive `H` call, as specified below. +Execute + + F(tail, prefix) + H_add(prefix, y) + F_inverse(tail, prefix) + H_add_inverse(prefix, y) + +Then undo the original causal lifts. There are two recursive calls. + +The recursive call restores its input tail. The expansion preserves the +whole source, including the changed prefix and the tail. Thus `F_inverse` +sees precisely the state left by `F` and restores the prefix exactly. +Target independence and rounded linearity are not required. The two +expansion errors are charged separately. At every depth the prefix and +tail are disjoint because `2w <= L`. This proves restoration by induction. + +The negative leaf reverses the complete middle word: first the positive +expansion, then `F`, then the inverse expansion, then `F_inverse`, then the +negative baseline addition. It is not the positive word with all signs +changed in place. + +In exact arithmetic the source lifts give the chain decomposition +`S x = z + H S_child z_tail`. The four instructions above add precisely the +second term. This is the usual dirty commutator, with its connector taken +from the source prefix. PR628 supplies the dirty commutator and exact grid +lift methods; PR846 supplies the chain/closure decomposition. + +## Reversible transport on two dirty tapes + +Write `C` for the causal map whose non-start coefficient is `a`. Split the +source and carrier views into consecutive blocks of at most `w` records. +The carrier has arbitrary contents. For each block in forward order: + +1. For each position after the first block, add the rounded product of `a` + and the preceding block's value on the carrier tape to the source. +2. Swap this whole source block with its corresponding carrier block. + +At this point the source tape holds the original carrier, and the carrier +tape holds `C_round(source)`. Swap the two complete views. The causal map +has been applied, and the carrier has been restored exactly. Each block +scan and each rewind costs its block length. There is no length-`w` clean +queue. The old carrier data remains in the other tape throughout. + +For the inverse, first swap the complete views. Visit the blocks in reverse +order. Swap the current block back, then subtract its rounded predecessor +products. The predecessor block is still on the carrier tape. Reverse the +order within each block as well. This is the exact inverse grid word. +Each map costs at most `64n` data-record moves on views of length `n`, +including read/write charges, block rewinds and full-view swaps. This loose +bound also covers a cyclic view that crosses the root tape boundary once. + +An expansion uses this causal map on the **target**, with the source as +its dirty carrier. Apply `C_round_inverse` to the target, add +`rnd(m*b_start*prefix_i)` to its first `w` entries, then apply `C_round`. +The source is restored by each transport pass and is unchanged at the +injection. The negative expansion changes only the injection sign. It is +the exact inverse of the positive expansion. + +The closure source and its dirty prefix initially lie on the same tape. +Swap the source prefix with the first `w` records of the other tape. Run +the child on the source tail and this moved prefix, then swap back. The +outer target's displaced data is never read by the child and is restored +exactly. Both child views are now on distinct tapes. Their physical spans +are contiguous. At the root an output shift of zero or the leaf stride +ensures that the moved prefix does not cross the tape boundary. All later +views are reversals or subintervals of these contiguous views. + +## Two coefficients per level + +No length-dependent clean coefficient table is required. Put + + g = gcd(L0,w0), M = L0/g, v = (w0/g)^-1 modulo M. + +For `M=1`, all cycles have length one and no modular inverse is required. +At a recursive level the remaining original indices form an interval, +with orientation `s` equal to `+1` or `-1`. Its length and stride are `n,d`. +The displacement to a predecessor is `s*d` when `i>=d`, and `s*(d-n)` +otherwise. The number of original edges is the unique integer in `[1,M]` +congruent to + + (displacement/g) v modulo M. + +The induced edge stops at the first remaining vertex. It cannot complete +an extra original cycle first. Thus these residues are the actual return +times, not only congruence classes. At a base level each vertex has returned +to itself, so its exponent is `M`. There are only two possible coefficients +per level, both powers of the original `b`. + +The checker independently forms each new exponent by summing the old chain +exponents. It compares all resulting entries with this modular formula. +That finite comparison supports the all-length argument above. + +The implementation recomputes these two powers before each pass. Scalar +temporaries are held only inside `causal`, `expand`, and `diagonal`; they +are discarded before a recursive call. The recursion retains only views, +integer metadata, and references to the immutable root descriptor. + +## Rounded coefficient preparation + +`power` uses binary powering from the most significant bit. After each +product it truncates both components toward zero. A truncation has complex +error at most `2u` and cannot increase the norm. Hence the value prepared +for exponent `e` has norm at most `q^e`, where `q=|b|`. + +Expand the repeated squaring into a binary multiplication tree with `e` +leaves. A local error propagates through sibling products whose combined +norm is at most `q^(e-t)`, where that node contains `t` leaves. At most two +nodes have more than `floor(e/2)` leaves. All other internal nodes have at +least `ceil(e/2)` leaves outside them. Therefore the error is at most + + 2u [2 + e q^ceil(e/2)] <= 2u [2 + 2/beta] <= 8u/beta. + +Here `k q^k <= sum_(j=0)^(k-1) q^j <= 1/beta`. Exponents one and two +satisfy the same bound directly. The tree counts repeated operands with +their multiplicity; reuse of a scalar result does not remove its error. + +Let a chain contain `t` prepared coefficients. Each coefficient and its +exact counterpart has norm at most `q`. A telescoping product comparison +gives error at most + + (8u/beta) t q^(t-1) <= 8u/beta^2. + +Compare this chain product `pi` with the freshly prepared child coefficient. +They approximate the same original power. Their difference is at most + + eta = 12u/beta^2, + +since `8u/beta^2 + 8u/beta <= 12u/beta^2`. Every prepared coefficient has +margin at least `beta`. These statements include complex coefficients. + +## Error recurrence + +Consider the exact solution of the system with the prepared coefficients +at each level. Suppose its right-hand side satisfies + + |x_i| <= A (1-|a_i|). + +The shared-suffix identity used by Charles Bethin in PR764 bounds the exact +zero-seed chain value by `A(1-|pi|)`. The rounded causal pass adds at most +`u/beta`. Because the child coefficient differs from `pi` by at most `eta`, +the next weighted bound is + + A_next <= (1 + 12u/beta^3) A + u/beta^2. + +Start with `A0=B/beta`. The admission condition gives +`K*12u/beta^3 <= 1/4`. Thus `(1+12u/beta^3)^K < 2`, and at every level + + A <= Astar = 2(B/beta+1). + +The inverse of a system with coefficients of modulus at most `1-beta` +has infinity norm at most `1/beta`. Its stronger weighted estimate is: +if `|v_i| <= C(1-|a_i|)`, its solution has norm at most `C`. Both follow +from the absolutely convergent Neumann series and telescoping weights. + +Compare the exact child solution with the exact Schur solution. Their +right-hand sides differ by at most `u/beta`, and their coefficients differ +by at most `eta`. The Schur solution is a restriction of the parent solution +and has norm at most `Astar`. Their solution difference is at most + + Dstar = u/beta^2 + 12 Astar u/beta^3. + +The rounded inverse causal map differs from `C^-1` by at most `u` per +entry: each non-start entry subtracts one rounded product. The forward map +differs from `C` by at most `u/beta`. The seed injection adds at most `u` +error. Since `||C||<=1/beta`, the conjugated expansion differs from its +linear additive map by at most `3u/beta`, independently of both dirty +arrays. There are two expansions. The source pass contributes another +`u/beta`. Each chain prefix contains at +least one coefficient, so it contracts a child error by at most `1-beta`. +For the unscaled leaf error, + + e_k <= (1-beta) e_(k+1) + Dstar + 7u/beta, + e_K <= u. + +Summing the geometric series gives + + e_0 <= u + u/beta^3 + 12 Astar u/beta^4 + 7u/beta^2 + <= 64(B+1)u/beta^5. + +The root multiplier scales the baseline and each seed injection. The +injection rounding was already charged above. The baseline adds at most +one further `u`; a diagonal root adds at most one. Taking `M>=1` gives the +stated conservative bound with an additional `3u` allowance. + +## Complete arithmetic and storage recurrence + +Let `n_k` be the level lengths after orientation changes. Then + + n_k >= 2 n_(k+1) + n_(k+2). + +Indeed the next length is the normalized remainder, which is at most the +ordinary Euclidean remainder. Let `S(n,d)` be the sum of lengths over the +full binary recursion tree. It satisfies + + S(n,0)=n, + S(n,d)=n+2 S(d,(-n) mod d) <= 5n+2d <= 6n. + +The potential `5n+2d` proves the inequality directly using +`n>=2d+d_next`. Each non-base visit pays for six complete causal maps, +two seed injections, the baseline output, four closure transport swaps, +and both complete children. Each base visit pays for all divisions and +target additions. Validation scans and final head returns are also paid. +Reversal and slicing create index views; their physical accesses still pay +for every movement. + +Each visit prepares at most sixteen powers, each with `O(log M)` products. +This cost is also `O(L)`, although the constant is not optimized. From +`n_k>=5n_(k+2)`, a tree of depth `K` has `O(2^K)` nodes and +`2^K=O(L^log_5(4))`. Since `log_5(4)<1`, +`2^K log L=O(L)`. The Euclidean and modular-inverse setup uses integer +arithmetic on `O(log L)` bits. The implementation retains `O(log L)` +integer metadata words, or `O(log^2 L)` bits, with no data-sized values +on the recursion stack. Index arithmetic can be charged conservatively +with sequential counter updates. A block counter costs `O(w)` bit changes +over a block, including its `O(log(w+1))` reset. Integer setup costs +`O(log^3 L)` per visit, using a conservative schoolbook division bound. +The bound on `2^K` also gives `2^K log^3 L=O(L)`, so the total control cost +is `O(L)` bit operations in +the stated tape implementation. The rational reference uses Python +indices; it is not a timing implementation of this control bound. + +An explicit loose count is at most `1280L` scalar products and divisions, +including power preparation and a root multiplier. To see this, let `J` +be the number of tree visits. The direct count is at most +`9S + 32J bit_length(M) + 3L`. Since `4^8 <= 5^7`, the depth bound gives +`J <= 4 L^(7/8)`. Also `log2 L <= 8 L^(1/8)` for `L>=1`. It follows that +`J bit_length(M) <= 36L`, which proves the stated count. Mutable data +accesses are at most `71S+2L < 512L`, including both validation scans. +These constants are arithmetic bounds, not measured runtime claims. + +Thus one leaf uses `O(L)` rounded scalar operations, two data arrays, and +`O(1)` clean arithmetic registers, apart from scalar arithmetic workspace. +If both incoming arrays have norm at most `D`, internal values have +magnitude bounded by a polynomial in `(D+1)/beta`, independent of `L`. +The weighted bound gives `H=3(D+1)/beta` for causal source values at every +level. Under `a*u<=1` below, a completed child changes a prefix of this +size by at most `H/beta+(H+1)/2`. Thus the changed seed has norm at most +`10(D+1)/beta^2`. Each expansion target adds at most `M` times that seed +norm plus `3u/beta`. Inside its causal conjugation, the inverse pass has +norm at most twice its input plus `u`; the injection adds its seed term; +the forward pass divides this combined bound by at most `beta` and adds +`u/beta`. Apply these inequalities to both expansion calls and the baseline. +The loose bound `512(M+1)(D+1)/beta^7` covers every intermediate value, +including the displaced carrier data. At each recursive call its initial +target is a source prefix or the completed change to that prefix; no +unbounded sequence of target magnifications is assumed. The implementation +is an exact-rational reference with explicit +rounding; it is not a fast fixed-precision arithmetic compiler. + +## Products and paid precision + +For a certified kernel `Q(U)=r^-1 product_j(U-rho_j)`, use the same factor +maps as PR628 and PR846, with `U=tau Z_L^w` and constant `|tau|=1`. +For an outer root use coefficient `tau/rho` and multiplier `-1/rho`. +For an inner root use coefficient `rho*conj(tau)`, stride `-w`, multiplier +`conj(tau)`, and output shift `-w`. Include `r` in one multiplier. + +The implementation uses PR628's balanced four-call dirty commutator. +It requires `h=ceil(log2(r-1))` distinct arbitrary dirty arrays. The new +leaf needs no additional array. For `n` factors, the complete leaf count is + + N(1)=1, N(n)=2N(floor(n/2))+2N(ceil(n/2)) <= 4n^2. + +A negative composition reverses the order and signs of the four child +calls. Ancestor inputs and targets cannot be used as deeper scratch. +All source and borrowed data restore exactly. The final rounded increment +of a composition can depend on the initial borrowed values, as in PR628. + +For a simple sufficient bound retain PR628's full-word perturbation proof. +Put `a=68(r+1)/beta^5`, `G=2+r/beta`, and let every initial array have norm +at most `D`. If `a*u<=1`, then the execution error is at most + + N a (D+1) u G^(N-1). + +Choose `W` at least the incoming grid precision, the admission budget, and + + p+1 + ceil(log2(N a (D+1))) + (N-1) ceil(log2 G). + +This pays the execution error to `2^(-p-1)`. The pinned factor certificate +pays the remaining error to `2^(-p-1)` after scaling by the source norm. +For general input norm `X`, coefficient preparation must use target +`p+ceil(log2(max(1,X)))`. The supplied fixtures have ample certified margin; +their scaled perturbation is checked exactly, not assumed. + +With PR628's inherited `log2(1/beta)=O(r^3)`, this conservative choice is +`O(p+incoming_bits+log(D+1)+r^5+log(log L+2))` fractional bits. Integer +magnitude bits have the same order. An explicit safe integer reserve is + + ceil(log2(D+1)) + N ceil(log2 G) + + ceil(log2(512(r+1)/beta^7)) + 3. + +This note does not claim a precision improvement. Gareth Davidson's PR877 +and Chafik Boukhalfa's PR905 give stronger prefix bounds for the predecessor +dispatcher. Their cubic or quartic conclusions are not imported without a +separate check of the new leaf contract. + +## Complete tape and preparation bill + +The meter maps every mutable array to one of **two** physical tapes and +an absolute tape offset. It charges the distance from the previous head +position plus one unit for each read or write. It counts transport swaps, +validation, source restoration, and explicit final head returns. It does +not give each borrowed array a free tape head. + +One non-base visit costs at most `640n` record moves, apart from its two +children. Six causal maps cost at most `384n`. Baseline and seed scans cost +at most `24n`. The closure swaps and travel into and out of the child spans +cost at most `32n`. These total less than `640n`. The loose allowance also +covers the root cyclic view and entry positioning. A base visit costs at +most `8n`. Since `S<=6L`, validation and head returns leave the complete +leaf below **`4096L` data-record moves**. This is a proof bound, not a claim +that its constant is sharp. + +For a full kernel with `h>=1`, place the last borrowed array on tape B. +Place the source, target, and other borrowed arrays consecutively on tape +A. If a leaf's two arrays are both on A, swap its target with the array on +B, execute the leaf, and swap back. The B array is distinct from both +operands in this case. Its contents never enter the leaf calculation and +are restored exactly. If a leaf already uses both tapes, no such routing +is needed. For `h=0`, put the two arrays on separate tapes directly. + +This routing adds `O((h+2)L)` moves per leaf, including the distance between +array slots and returns to the tape origins. The full bound used by the +checker is `8192(N+1)(h+2)L` moves. It includes the initial validation scan. +No extra data array is allocated for transport. + +Let `W_total` include the fractional bits, integer reserve, and the bit +lengths of the immutable coefficient, gauge, and margin descriptors. The generic +rational API does not bound descriptor length from coefficient magnitude. +Let `M(W_total)` include the chosen rounded scalar multiplication and +division cost. With this complete scalar word size, the resulting bill is + + preparation + O(N L) M(W_total) + O(N(h+2)L W_total). + +Locating a factor on a sequential read-only descriptor tape costs at most +`O(r W_total)` per leaf, hence `O(r^3 W_total)` over the word. The original +factor and phase preparation charges remain inherited from PR628; this +descriptor access is included, not treated as free. Temporary powers, +source reads, multipliers, divisions, copies, inverse calls, and children +are all charged. The two data tapes, scalar arithmetic tapes, and metadata +stack use a fixed number of tapes. Borrowed data storage remains `hL` +complex words. Scalar arithmetic workspace has the declared word size; +the leaf index stack has `O(log^2 L)` bits. The balanced composition and +array routing retain another `O(log r (log r + log L))` metadata bits. +Their control cost is absorbed by the routing and descriptor terms, since +`log L=O(L)` and `log^2 r=O(r)`. These are metadata bits, not arithmetic +registers or data arrays. + +This is a constant-coefficient kernel construction. Variable site phases, +complete-transform composition and normalization, common-label/Clifford +compilation, physical supplier hosting, and the multiplication reduction +remain open. There is no new kappa, complex b, or BIT c claim. Finite replay +is not a formal proof of the full theorem or independent human review. + +Chafik Boukhalfa's concurrent PR951 gives two grouped additive calls with +no dirty connector in a neighboring-cycle layout. It uses a larger clean +scalar workspace. This note instead uses constant clean arithmetic storage +per first-order leaf, the original logarithmic dirty array budget, and +arbitrary-stride transport. It makes no claim to improve PR951's arithmetic +count or cubic precision. + +Chafik Boukhalfa's PR977 removes the length admission and proves cubic +precision for PR846 in its original clean scratch model. This construction +retains its stated admission and conservative precision. Combining those +precision arguments with the new dirty leaf requires a separate proof. diff --git a/research/dirty-stride-fast-phases/verify.py b/research/dirty-stride-fast-phases/verify.py new file mode 100644 index 0000000000..a039052e94 --- /dev/null +++ b/research/dirty-stride-fast-phases/verify.py @@ -0,0 +1,307 @@ +"""Check phase accuracy, exact cleanup, dirty calls, and the paid work.""" +import argparse +from fractions import Fraction as F +import hashlib +import importlib.util +import json +from pathlib import Path +import sys + +HERE = Path(__file__).resolve().parent + + +def load(path, name): + spec = importlib.util.spec_from_file_location(name, path) + module = importlib.util.module_from_spec(spec) + sys.modules[name] = module + spec.loader.exec_module(module) + return module + + +c = load(HERE/'chirp.py', 'checked_fast_chirp') +s = c.stream +o = load(HERE/'vendor/pr1038/oracle.py', 'fast_phase_oracle') + + +def need(condition, message): + if not condition: + raise ValueError(message) + + +def canonical(value): + return json.dumps(value, sort_keys=True, separators=(',', ':')) + + +def digest(values): + return hashlib.sha256(canonical([[str(a), str(b)] for a, b in values]).encode()).hexdigest() + + +def pins(): + source = json.loads((HERE/'SOURCE.json').read_text()) + for path, row in source['files'].items(): + value = (HERE/path).read_bytes() + need(hashlib.sha256(value).hexdigest() == row['sha256'], 'source SHA256: '+path) + blob = hashlib.sha1(b'blob '+str(len(value)).encode()+b'\0'+value).hexdigest() + need(blob == row['git_blob'], 'source Git blob: '+path) + return len(source['files']) + + +def coefficient_check(value, coefficients, ideal, bits): + quadrant, tangent, sine = coefficients + real, imag = ideal + reduced = ((real, imag), (imag, o.neg(real)), + (o.neg(real), o.neg(imag)), (o.neg(imag), real))[quadrant] + need(reduced[0][0] >= F(3, 5), 'phase quadrant') + denominator = o.plus(o.point(1), reduced[0], bits+96) + quotients = [v/w for v in reduced[1] for w in denominator] + error_t = max(abs(tangent-min(quotients)), abs(tangent-max(quotients))) + error_s = max(abs(sine-reduced[1][0]), abs(sine-reduced[1][1])) + need(max(error_t, error_s) <= F(1, 1 << (bits+8)), 'phase coefficient accuracy') + bound = 4*(abs(value[0])+abs(value[1])+1)*F(1, 1 << bits) + outputs = [] + for inverse in (False, True): + result = s.rotate(value, coefficients, bits, inverse) + phase = (ideal[0], o.neg(ideal[1])) if inverse else ideal + expected = o.complex_mul(phase, o.complex_point(value), bits+96) + need(o.squared_distance_bound(result, expected) <= bound**2, 'grid rotation accuracy') + need(s.rotate(result, coefficients, bits, not inverse) == value, 'grid rotation inverse') + outputs.append(result) + return digest(outputs) + + +def phase_case(numerator, modulus, bits, length, sample_stride=1): + meter = s.Meter() + plan = s.Phase.prepare(numerator, modulus, bits, meter, length) + state = plan.start() + unit = F(1, 1 << plan.state_bits) + value = F(2**80)+F(1, 1 << bits), -F(2**61)+F(3, 1 << bits) + result = hashlib.sha256() + samples = 0 + maximum_ratio = F(0) + for j in range(length): + if j % sample_stride == 0 or j == length-1: + ideal = o.phase(numerator, modulus, j, plan.state_bits+96) + ratio = o.phase(numerator*(2*j+1), modulus, 1, plan.state_bits+96) + phase_bound = 256*(j+1)**2*unit + ratio_bound = 12*(j+1)*unit + actual_error = o.squared_distance_bound(state.value, ideal) + need(actual_error <= phase_bound**2, 'phase state drift') + need(o.squared_distance_bound(state.ratio, ratio) <= ratio_bound**2, 'ratio state drift') + maximum_ratio = max(maximum_ratio, actual_error/phase_bound**2) + coefficients = state.coefficients() + result.update(coefficient_check(value, coefficients, ideal, bits).encode()) + need(state.erase_coefficients(coefficients) == (0, 0, 0), 'emitted coefficient zero') + samples += 1 + state.step() + state.clear() + need(state.value == state.ratio == c.d.ZERO, 'phase state zero') + plan.clear() + need(plan.seed == c.d.ZERO and plan.omega == (0, 0, 0), 'phase seed zero') + need(meter.stream_steps == meter.stream_reverse_steps == length, 'paid reverse reconstruction') + need(meter.stream_charts == 4*length+2*samples, 'paid coefficient reconstruction') + need(meter.phase_products == 12*length+6, 'paid phase state rotations') + need(meter.seed_cleanups == 6 and meter.stream_cleanups == meter.preparation_cleanups == 1, + 'complete phase cleanup') + return dict(numerator=str(numerator), modulus=str(modulus), bits=bits, length=length, + sample_stride=sample_stride, samples=samples, state_bits=meter.state_bits, + max_squared_drift_ratio_log2=c.d.ceil_log2(maximum_ratio) if maximum_ratio else None, + seed_peak_bits=meter.seed_peak_bits, output_digest=result.hexdigest()) + + +def phases(): + rows = [] + for bits in (2, 24): + for modulus in range(1, 10): + for numerator in range(modulus): + rows.append(phase_case(numerator, modulus, bits, 2*modulus+1)) + for bits, numerator, modulus, length in ( + (64, -3, 16, 33), (64, 5, 17, 41), (128, -17, 31, 65), + (96, 1, 2**160+7, 9), (96, 2**180+9, 2**160+7, 9), + (96, 2**157-1, 2**160, 11), (96, 2**157+1, 2**160, 11)): + rows.append(phase_case(numerator, modulus, bits, length)) + rows.append(phase_case(12345, 65537, 24, 16384, 257)) + return rows + + +class Bus: + def __init__(self): + self.heads = {0: 0, 1: 0} + self.minimum_moves = self.reads = self.writes = 0 + + def access(self, tape, position, write): + self.minimum_moves += abs(position-self.heads[tape])+1 + self.heads[tape] = position + if write: + self.writes += 1 + else: + self.reads += 1 + + +class Bank: + def __init__(self, values, bus, tape, offset=0): + self.values = list(values) + self.bus, self.tape, self.offset = bus, tape, offset + + def __len__(self): + return len(self.values) + + def __iter__(self): + raise ValueError('array copies are not allowed') + + def __getitem__(self, index): + need(type(index) is int, 'array slices are not allowed') + self.bus.access(self.tape, self.offset+index, False) + return self.values[index] + + def __setitem__(self, index, value): + need(type(index) is int, 'array slices are not allowed') + self.bus.access(self.tape, self.offset+index, True) + self.values[index] = value + + +def data(length, bits, salt, dirt=0): + unit = F(1, 1 << bits) + return [(dirt+F((7*j+salt) % 19-9, 8)+unit, + -dirt+F((11*j+salt) % 17-8, 16)+3*unit) for j in range(length)] + + +def word_case(length, stride, numerator, modulus, dirt, count=1, sign=1): + beta = F(1, 4) + factors = [(j % 2 == 0, (F(1, 4), F((-1)**j, 8)), F(1)) for j in range(count)] + tau = F(3, 5), F(4, 5) + bits = c.working_bits(24, 4, 2*abs(dirt)+16, beta, length, stride, + count=count, scale_bound=1 if count > 1 else None) + height = (count-1).bit_length() + initial = [data(length, bits, 3+7*j, dirt if j else 0) for j in range(height+2)] + bus = Bus() + banks = [Bank(a, bus, 1 if j == len(initial)-1 else 0, + 0 if j == len(initial)-1 else j*length) for j, a in enumerate(initial)] + x, y, scratch = banks[0], banks[1], banks[2:] + def call(direction): + if count == 1: + return c.lift(x, y, (F(3, 10), F(2, 5)), beta, bits, stride, + numerator, modulus, direction) + return c.compose(x, y, factors, scratch, tau, beta, bits, stride, + numerator, modulus, direction) + meter = call(sign) + need(x.values == initial[0] and [a.values for a in scratch] == initial[2:], 'dirty source restoration') + need(meter.reads == bus.reads and meter.writes == bus.writes, 'literal data access bill') + need(meter.tape_moves >= bus.minimum_moves, 'physical data movement bill') + need(set(meter.heads) == {0, 1} and all(v == 0 for v in meter.heads.values()), 'returned data heads') + calls = c.d.leaf_calls(count) + bill = 4160*length if count == 1 else 8256*(calls+1)*(height+2)*length + need(meter.tape_moves <= bill, 'complete two-tape bound') + need(meter.phase_evaluations == meter.stream_steps == meter.stream_reverse_steps == 4*length, + 'four phase streams') + need(meter.stream_cleanups == 4 and meter.preparation_cleanups == 1, 'word phase cleanup') + need(meter.phase_products == 60*length+6 and meter.phase_divisions == 24*length, + 'linear phase arithmetic') + need(meter.error_bound <= F(1, 1 << 25), 'complete dirty precision') + if count == 1: + exact = o.leaf(initial[0], (F(3, 10), F(2, 5)), stride, numerator, modulus, bits+96) + else: + exact = o.kernel(initial[0], factors, tau, stride, numerator, modulus, bits+96) + error = F(0) + nonzero = False + for j in range(length): + increment = tuple(sign*(v-w) for v, w in zip(y.values[j], initial[1][j])) + error = max(error, o.squared_distance_bound(increment, exact[j])) + nonzero |= any(not lo <= v <= hi for v, (lo, hi) in zip(increment, exact[j])) + need(error <= meter.error_bound**2, 'dirty operator accuracy') + need(nonzero, 'nonzero rounded execution error') + output = digest(y.values) + call(-sign) + need([a.values for a in banks] == initial, 'complete dirty inverse') + return dict(length=length, stride=stride, count=count, sign=sign, dirt=str(dirt), bits=bits, + borrowed_arrays=height, calls=calls, tape_moves=meter.tape_moves, + seed_terms=meter.seed_terms, products=meter.phase_products, + divisions=meter.phase_divisions, output_digest=output) + + +def words(): + rows = [] + for length, stride, dirt in ((1, 0, 7), (5, 2, -13), (8, -3, 2**40), (17, 7, 2**80)): + for count in (1, 3): + for sign in (-1, 1): + rows.append(word_case(length, stride, -3, 2*length, dirt, count, sign)) + return rows + + +def costs(): + rows = [] + for length, bits in ((32, 16), (257, 48), (2048, 96)): + before = data(length, bits, 3) + a, b = list(before), list(before) + old_meter, new_meter = s.old.Meter(), s.Meter() + old_meter.layout = {id(a): (0, 0)} + new_meter.layout = {id(b): (0, 0)} + old_plan = s.old.Phase.prepare(-3, 2*length+1, bits, old_meter) + new_plan = s.Phase.prepare(-3, 2*length+1, bits, new_meter, length) + for inverse in (False, True): + s.old.phase_pass(a, old_plan, inverse) + s.phase_pass(b, new_plan, inverse) + new_plan.clear() + need(a == b == before, 'same-input baseline restoration') + need(old_meter.reads == new_meter.reads and old_meter.writes == new_meter.writes + and old_meter.tape_moves == new_meter.tape_moves, 'same-input tape equality') + need(new_meter.phase_products == 30*length+6 and new_meter.phase_divisions == 12*length, + 'linear stream work') + need(new_meter.seed_terms <= 48*new_meter.internal_bits, + 'seed work independent of stream length') + rows.append(dict(length=length, bits=bits, old_terms=old_meter.phase_terms, + old_products=old_meter.phase_products, old_divisions=old_meter.phase_divisions, + new_seed_terms=new_meter.seed_terms, new_seed_products=new_meter.seed_products, + new_seed_divisions=new_meter.seed_divisions, new_stream_products=new_meter.phase_products, + new_stream_divisions=new_meter.phase_divisions, new_seed_peak_bits=new_meter.seed_peak_bits, + new_state_bits=new_meter.state_bits, tape_moves=new_meter.tape_moves)) + return rows + + +def controls(): + rows = [] + def reject(label, call, expected): + try: + call() + except ValueError as error: + need(str(error) == expected, 'wrong rejection: '+label) + rows.append(label) + else: + raise ValueError('accepted invalid input: '+label) + for label, args, expected in ( + ('zero_modulus', (1, 0, 12, 8), 'integer phase descriptor'), + ('boolean_numerator', (True, 7, 12, 8), 'integer phase descriptor'), + ('zero_length', (1, 7, 12, 0), 'phase stream domain'), + ('boolean_precision', (1, 7, True, 8), 'phase stream domain')): + a, m, w, length = args + reject(label, lambda a=a, m=m, w=w, length=length: + s.Phase.prepare(a, m, w, s.Meter(), length), expected) + p = s.Phase.prepare(1, 8, 12, s.Meter(), 3) + state = p.start() + reject('second_stream', p.start, 'one open phase stream') + reject('early_seed_clear', p.clear, 'phase stream must be clear') + state.step() + state.clear() + p.clear() + reject('cleared_seed', p.start, 'active phase preparation') + return rows + + +def main(): + parser = argparse.ArgumentParser() + parser.add_argument('--write', type=Path) + args = parser.parse_args() + result = dict(schema=1, source_pins=pins(), phases=phases(), words=words(), costs=costs(), controls=controls()) + encoded = json.dumps(result, indent=2, sort_keys=True)+'\n' + if args.write: + with args.write.open('x') as output: + output.write(encoded) + else: + need(canonical(result) == canonical(json.loads((HERE/'EXPECTED.json').read_text())), 'full typed receipt') + print('PASS_DIRTY_STRIDE_FAST_PHASES phases='+str(len(result['phases'])) + +' samples='+str(sum(row['samples'] for row in result['phases'])) + +' words='+str(len(result['words']))+' costs='+str(len(result['costs'])) + +' input_controls='+str(len(result['controls']))) + + +if __name__ == '__main__': + main()