diff --git a/docs/holomorphic-mathematical-contract.md b/docs/holomorphic-mathematical-contract.md index b83e8f6..00b528e 100644 --- a/docs/holomorphic-mathematical-contract.md +++ b/docs/holomorphic-mathematical-contract.md @@ -4,6 +4,17 @@ This document fixes the mathematics of the live explorer independently of CPU th The renderer may approximate these formulas numerically. It must not redefine the mathematical object in order to fit a particular implementation. +## Status labels used by this contract + +To keep theorem, prototype, and taste from drifting together, statements in this repository have these roles: + +- **MATHEMATICAL CONTRACT** — invariants every implementation must preserve, including `f = R exp(q)`, entire `q`, divisor preservation, and the phase/log-modulus handoff. +- **CURRENT REFERENCE IMPLEMENTATION** — the small polynomial `q`, coefficient envelope, and host/shader evaluators used today. +- **HISTORICAL EXPERIMENT** — the unit-disc Bergman workers and the old three-workers-plus-coordinator organization. +- **CANDIDATE FUTURE MATHEMATICS** — the Bargmann-Fock coordinate system below, pending host comparison and visual judgment. +- **BACKEND OPTIMIZATION IDEA** — scheduling, precision, descriptor packing, and CPU/GPU execution strategies; these do not define admissible functions. +- **VISUAL / ART-DIRECTION CHOICE** — scale, anchor distribution, amplitude, phase/sign, lifetime, temporal correlation, overlap, active count, and cadence. + ## 1. The live object The explorer displays @@ -37,7 +48,9 @@ The current implementation uses a small polynomial `q_t`. That is a finite compu ## 2. Why `exp(q)` is structural rather than cosmetic -If two meromorphic functions on `C` have the same finite zeros and poles with the same multiplicities, their quotient has no zeros or poles. When that quotient is entire and nonvanishing, it has an entire logarithm because `C` is simply connected. Thus the quotient can be written as `exp(q)` for an entire `q`. +If two nonzero meromorphic functions on `C` have the same finite zeros and poles with the same multiplicities, their quotient extends across the cancelled divisor to an entire, nonvanishing function. Every nonvanishing entire function `H` has an entire logarithm: `H'/H` is entire, its integral on the simply connected plane is an entire function `g`, and `H exp(-g)` is a nonzero constant whose logarithm can be absorbed into `g`. Thus `H = exp(q)` for an entire `q`. + +The simply connected hypothesis matters. A nonvanishing holomorphic function on an arbitrary multiply connected domain need not have a single-valued holomorphic logarithm. This contract makes the claim on `C`, not on every domain. So @@ -72,7 +85,7 @@ The present polynomial prototype omits the constant term, equivalently fixing `q There is no mathematical reason for the renderer to evaluate the complex exponential explicitly. -For +Away from the explicit zeros and poles, for ```text f(z) = R(z) exp(q(z)) @@ -82,7 +95,7 @@ we have exactly ```text log|f(z)| = log|R(z)| + Re(q(z)) -phase(f(z)) = phase(R(z)) + Im(q(z)) +phase(f(z)) = phase(R(z)) + Im(q(z)) (mod 2 pi) ``` The mathematical renderer boundary is therefore @@ -114,7 +127,7 @@ A local derivative may also be used when the desired perturbation is stated in t A canonical direction does not require a visual-energy heuristic. -Choose a Hilbert space `A` of admissible holomorphic perturbations in which point evaluation is continuous. Let `K(z,a)` be its reproducing kernel. If we require a perturbation `phi` to satisfy +Choose a Hilbert space `A` of admissible holomorphic perturbations with a stated norm, gauge, and continuous point evaluation. Use the convention `h(a) = `, with the inner product linear in its first argument. If we require a perturbation `phi` to satisfy ```text phi(a) = 1 @@ -128,6 +141,21 @@ phi_a(z) = K(z,a) / K(a,a). This is the precise meaning of a canonical direction of least holomorphic disturbance for a prescribed local value change. +Indeed, `||K(.,a)||^2 = K(a,a)`. Every admissible `h` with `h(a) = 1` decomposes as + +```text +h = K(.,a) / K(a,a) + g, +g(a) = 0. +``` + +The reproducing identity makes `g` orthogonal to `K(.,a)`, so + +```text +||h||^2 = 1 / K(a,a) + ||g||^2. +``` + +The normalized representer is therefore the unique minimum-norm solution. “Canonical” means canonical relative to the declared Hilbert space, norm, local functional, and gauge; it does not mean independent of those modeling choices. + More generally, if the prescribed local datum is a derivative or another continuous linear functional, its Riesz/reproducing representer gives the corresponding unique minimum-norm direction after normalization. Thus the mathematical pipeline can be @@ -144,13 +172,25 @@ The user remains the visual arbiter of anchor selection, amplitude, timing, over ## 8. Historical Bergman-disk construction -The historical `local-holomorphic-perturbations` experiment used the unit-disc Bergman extremal +The historical `local-holomorphic-perturbations` experiment used the unit-disc Bergman space + +```text +A^2(D) = {h holomorphic on D : (1/pi) integral_D |h(z)|^2 dA(z) < infinity}. +``` + +For this normalized area convention its kernel is + +```text +K_D(z,a) = 1 / (1 - conjugate(a) z)^2, +``` + +and its value-normalized extremal is ```text phi_a(z) = (1 - |a|^2)^2 / (1 - conjugate(a) z)^2 ``` -which satisfies `phi_a(a) = 1` and is the minimum Bergman-norm holomorphic function on the disc with that value. +which satisfies `phi_a(a) = 1` and is the minimum Bergman-norm holomorphic function on the disc with that value. Using unnormalized area instead multiplies `K_D` by `1/pi` and leaves `K_D(z,a)/K_D(a,a)` unchanged. That correctly demonstrated the canonical-direction idea on a bounded disc. It is not, by itself, the whole-plane basis for this repository: for `a != 0` it has a pole at `1 / conjugate(a)` outside the unit disc, and `exp(phi_a)` would have an essential singularity there. @@ -160,13 +200,21 @@ Do not revive that hidden singularity merely because the kernel was useful in th If we want the same extremal construction using entire functions, the admissible Hilbert space must itself consist of entire functions. -A natural candidate is a Bargmann-Fock space with an explicit length scale `s`. Its reproducing kernel has the form +A natural candidate is the Bargmann-Fock space, for `s > 0`, with the explicit convention ```text -K_s(z,a) = exp(z conjugate(a) / s^2) +F_s^2 = {h entire : ||h||_s^2 < infinity}, +||h||_s^2 = (1 / (pi s^2)) integral_C |h(z)|^2 exp(-|z|^2 / s^2) dA(z). ``` -up to the chosen normalization convention. The value-normalized extremal is therefore +The functions `z^n / (s^n sqrt(n!))` form an orthonormal basis for this convention, so the reproducing kernel is exactly + +```text +K_s(z,a) = exp(z conjugate(a) / s^2), +K_s(a,a) = exp(|a|^2 / s^2). +``` + +The value-normalized extremal is therefore ```text phi_a(z) = exp((z conjugate(a) - |a|^2) / s^2) @@ -174,13 +222,41 @@ phi_a(z) = exp((z conjugate(a) - |a|^2) / s^2) which is entire and satisfies `phi_a(a) = 1`. +Its squared norm is `exp(-|a|^2 / s^2)`. This is minimum disturbance in the stated Gaussian-weighted Fock norm. It is not a proof of visual localization in an unweighted viewport: its modulus can grow in the direction selected by the anchor. + This is a mathematically clean whole-plane candidate for canonical local perturbations. The choice of function-space norm and the scale `s` are modeling choices, not universal aesthetic truths. They should be exposed to visual evaluation rather than smuggled in as GPU constants. -If the gauge removes constant motion, use the corresponding gauge-fixed subspace or a derivative constraint rather than silently reintroducing the constant mode. +The displayed formula is **ungauged**. In particular, `phi_a(0) = exp(-|a|^2 / s^2)`, so it does not belong to the present prototype's `q(0) = 0` gauge. If that gauge is retained, use the closed subspace + +```text +F_{s,0}^2 = {h in F_s^2 : h(0) = 0} +``` + +whose kernel is + +```text +K_{s,0}(z,a) = exp(z conjugate(a) / s^2) - 1. +``` + +For a nonzero anchor, its normalized value representer is + +```text +phi_{a,0}(z) + = (exp(z conjugate(a) / s^2) - 1) + / (exp(|a|^2 / s^2) - 1). +``` + +It satisfies both `phi_{a,0}(0) = 0` and `phi_{a,0}(a) = 1`. A value constraint `h(0) = 1` is incompatible with this gauge. An origin derivative constraint is compatible; the minimum-norm direction satisfying `h'(0) = 1` is `h(z) = z`. The exact whole-plane perturbation space is therefore a mathematical/design decision to settle before optimizing a GPU implementation. The invariant that its elements are entire is not optional. -## 10. What randomness means +## 10. Superposition + +Every finite sum of admitted entire perturbation descriptors is entire. Therefore several value or derivative disturbances may overlap in `q` before exponentiation without creating finite singularities or changing the explicit divisor. + +For an infinite series, termwise entireness alone is insufficient: locally uniform convergence on `C` is a sufficient condition for the sum to remain entire. The current implementation and the candidate host reference use finite superpositions. + +## 11. What randomness means Randomness may choose: @@ -194,7 +270,7 @@ Randomness does not certify holomorphy. The admitted function space and exact fo The number of CPU workers is not mathematics. Three workers plus a coordinator was a useful implementation shape for a four-thread phone, but another CPU, x86-64 implementation, or GPU backend may schedule the same mathematical descriptors differently. -## 11. Current polynomial prototype +## 12. Current polynomial prototype The current implementation uses @@ -221,7 +297,7 @@ So the bound is usefully interpreted as an amplitude/numerical envelope on the r The present 128-candidate, three-worker score is a provisional CPU exploration strategy. It is not the definition of the canonical holomorphic direction and must not become part of the mathematical semantics merely because it exists in working code. -## 12. Safe motion descriptors +## 13. Safe motion descriptors If an implementation publishes a segment @@ -236,7 +312,7 @@ If the accepted coefficient set is convex, such as the current `sum |c_k| <= B` Other numerical, amplitude, derivative, or application-specific bounds must be named separately. Do not call them holomorphy checks. -## 13. Separation of responsibilities +## 14. Separation of responsibilities The mathematical evolution engine owns: @@ -256,7 +332,7 @@ Wegert owns the reusable rendering preference boundary from complex value / phas The user owns the final visual judgment about what motion is worth keeping. -## 14. Acceptance before GPU optimization +## 15. Acceptance before GPU optimization Before backend-specific optimization, host/reference tests should establish at least: @@ -269,4 +345,13 @@ Before backend-specific optimization, host/reference tests should establish at l 7. safe segments remain inside every claimed convex bound; 8. the current CPU heuristic is clearly labeled as an approximation/scheduling strategy rather than the mathematical definition. -Only after those are stable should PowerVR, FP16/FP32, fragment/compute division, register pressure, or other backend details be allowed to influence implementation choices. \ No newline at end of file +Only after those are stable should PowerVR, FP16/FP32, fragment/compute division, register pressure, or other backend details be allowed to influence implementation choices. + +The dependency-free host oracle in [`reference/entire_representer.py`](../reference/entire_representer.py) evaluates the two Fock value conventions, derivative representers, finite descriptor sums, and the complete regular-point state `R(z) exp(q(z))`. It is a **CANDIDATE FUTURE MATHEMATICS** oracle, not Android runtime code and not approval of a visual parameter set. + +## 16. Mathematical references + +- N. Aronszajn, “Theory of Reproducing Kernels,” *Transactions of the American Mathematical Society* 68 (1950), 337–404, [doi:10.1090/S0002-9947-1950-0051437-7](https://doi.org/10.1090/S0002-9947-1950-0051437-7). +- V. Bargmann, “On a Hilbert Space of Analytic Functions and an Associated Integral Transform,” *Communications on Pure and Applied Mathematics* 14 (1961), 187–214, [doi:10.1002/cpa.3160140303](https://doi.org/10.1002/cpa.3160140303). + +The repository named `isomorphismes/Conway` is a planar-symmetry/orbifold project, not a source for these complex-analysis results; no Conway chapter citation is asserted here. diff --git a/reference/entire_representer.py b/reference/entire_representer.py new file mode 100644 index 0000000..b81cb2b --- /dev/null +++ b/reference/entire_representer.py @@ -0,0 +1,233 @@ +"""Host oracle for candidate whole-plane canonical perturbations. + +This module is deliberately independent of Android, GLSL, and Idric backend +types. It fixes one Bargmann--Fock convention so that proposed descriptors can +be checked before a backend optimizes them. It is candidate mathematics, not +the live motion policy and not a choice of visual parameters. + +The norm convention is + + ||h||_s^2 = 1/(pi s^2) integral_C |h(z)|^2 exp(-|z|^2/s^2) dA(z). + +All state evaluation is at regular points of the explicit meromorphic factor. +""" + +from __future__ import annotations + +from dataclasses import dataclass +import cmath +import math +from typing import Iterable, Literal + + +Gauge = Literal["none", "zero_at_origin"] +DatumKind = Literal["value", "derivative"] + + +def _scale_squared(scale: float) -> float: + if not math.isfinite(scale) or scale <= 0.0: + raise ValueError("the Bargmann--Fock scale must be finite and positive") + return scale * scale + + +def fock_kernel(z: complex, anchor: complex, scale: float) -> complex: + """Return K_s(z,a) = exp(z conjugate(a) / s^2).""" + + scale_squared = _scale_squared(scale) + return cmath.exp(z * anchor.conjugate() / scale_squared) + + +def fock_value_direction( + z: complex, + anchor: complex, + scale: float, + *, + gauge: Gauge, +) -> complex: + """Return the minimum-norm direction whose value at ``anchor`` is one. + + ``gauge="none"`` uses the full Fock space. ``zero_at_origin`` uses the + closed subspace h(0)=0 and therefore rejects a value constraint at zero. + """ + + scale_squared = _scale_squared(scale) + anchor_squared = abs(anchor) ** 2 + + if gauge == "none": + return cmath.exp( + (z * anchor.conjugate() - anchor_squared) / scale_squared + ) + if gauge == "zero_at_origin": + if anchor == 0.0: + raise ValueError("h(0)=1 is incompatible with the h(0)=0 gauge") + numerator = cmath.exp(z * anchor.conjugate() / scale_squared) - 1.0 + denominator = math.expm1(anchor_squared / scale_squared) + return numerator / denominator + raise ValueError(f"unknown gauge: {gauge}") + + +def fock_derivative_direction( + z: complex, + anchor: complex, + scale: float, +) -> complex: + """Return the minimum-norm direction whose derivative at ``anchor`` is one. + + The derivative representer is + + r_a(z) = (z/s^2) exp(z conjugate(a)/s^2). + + It already vanishes at the origin, so projection to the h(0)=0 subspace + does not change it. + """ + + scale_squared = _scale_squared(scale) + exponent = abs(anchor) ** 2 / scale_squared + representer = (z / scale_squared) * cmath.exp( + z * anchor.conjugate() / scale_squared + ) + representer_derivative_at_anchor = ( + (scale_squared + abs(anchor) ** 2) + / (scale_squared * scale_squared) + * math.exp(exponent) + ) + return representer / representer_derivative_at_anchor + + +def minimum_norm_squared( + *, + anchor: complex, + scale: float, + kind: DatumKind, + gauge: Gauge, +) -> float: + """Return the squared norm of the corresponding normalized representer.""" + + scale_squared = _scale_squared(scale) + exponent = abs(anchor) ** 2 / scale_squared + + if kind == "value": + if gauge == "none": + return math.exp(-exponent) + if gauge == "zero_at_origin": + if anchor == 0.0: + raise ValueError("h(0)=1 is incompatible with the h(0)=0 gauge") + return 1.0 / math.expm1(exponent) + raise ValueError(f"unknown gauge: {gauge}") + + if kind == "derivative": + if gauge not in ("none", "zero_at_origin"): + raise ValueError(f"unknown gauge: {gauge}") + return ( + scale_squared + * scale_squared + * math.exp(-exponent) + / (scale_squared + abs(anchor) ** 2) + ) + + raise ValueError(f"unknown local datum kind: {kind}") + + +@dataclass(frozen=True) +class LocalDatum: + """One prescribed local value or derivative contribution to q.""" + + anchor: complex + amplitude: complex + scale: float + kind: DatumKind + gauge: Gauge + + def direction(self, z: complex) -> complex: + if self.kind == "value": + return fock_value_direction( + z, self.anchor, self.scale, gauge=self.gauge + ) + if self.kind == "derivative": + if self.gauge not in ("none", "zero_at_origin"): + raise ValueError(f"unknown gauge: {self.gauge}") + return fock_derivative_direction(z, self.anchor, self.scale) + raise ValueError(f"unknown local datum kind: {self.kind}") + + def contribution(self, z: complex) -> complex: + return self.amplitude * self.direction(z) + + +def evaluate_q(z: complex, data: Iterable[LocalDatum]) -> complex: + """Evaluate a finite superposition of admitted entire descriptors.""" + + return sum((datum.contribution(z) for datum in data), start=0j) + + +@dataclass(frozen=True) +class DivisorPoint: + location: complex + multiplicity: int = 1 + + def __post_init__(self) -> None: + if not isinstance(self.multiplicity, int) or isinstance( + self.multiplicity, bool + ) or self.multiplicity <= 0: + raise ValueError("divisor multiplicity must be a positive integer") + + +def evaluate_r( + z: complex, + *, + zeros: Iterable[DivisorPoint] = (), + poles: Iterable[DivisorPoint] = (), + gain: complex = 1.0, +) -> complex: + """Evaluate the explicitly factored meromorphic part R at a regular point.""" + + if gain == 0.0: + raise ValueError("the meromorphic factor gain must be nonzero") + value = complex(gain) + for zero in zeros: + value *= (z - zero.location) ** zero.multiplicity + for pole in poles: + value /= (z - pole.location) ** pole.multiplicity + return value + + +@dataclass(frozen=True) +class RegularPointState: + q: complex + re_q: float + im_q: float + r: complex + f: complex + log_modulus: float + phase: float + + +def evaluate_regular_state( + z: complex, + data: Iterable[LocalDatum], + *, + zeros: Iterable[DivisorPoint] = (), + poles: Iterable[DivisorPoint] = (), + gain: complex = 1.0, +) -> RegularPointState: + """Evaluate q, R exp(q), and the exact phase/log-modulus handoff. + + ``phase`` is an unwrapped representative. Rendering consumes it modulo + 2*pi. A point on the explicit divisor is rejected because phase and finite + log modulus are not defined there. + """ + + q = evaluate_q(z, data) + r = evaluate_r(z, zeros=zeros, poles=poles, gain=gain) + if r == 0.0: + raise ValueError("phase/log-modulus evaluation requires a regular point") + + f = r * cmath.exp(q) + return RegularPointState( + q=q, + re_q=q.real, + im_q=q.imag, + r=r, + f=f, + log_modulus=math.log(abs(r)) + q.real, + phase=cmath.phase(r) + q.imag, + ) diff --git a/tests/test_entire_representer.py b/tests/test_entire_representer.py new file mode 100644 index 0000000..a047449 --- /dev/null +++ b/tests/test_entire_representer.py @@ -0,0 +1,207 @@ +from __future__ import annotations + +import cmath +import math +from pathlib import Path +import sys +import unittest + + +ROOT = Path(__file__).resolve().parents[1] +sys.path.insert(0, str(ROOT / "reference")) + +from entire_representer import ( # noqa: E402 + DivisorPoint, + LocalDatum, + evaluate_q, + evaluate_regular_state, + fock_derivative_direction, + fock_kernel, + fock_value_direction, + minimum_norm_squared, +) + + +def polynomial_fock_norm_squared(coefficients: list[complex], scale: float) -> float: + return sum( + abs(coefficient) ** 2 * scale ** (2 * degree) * math.factorial(degree) + for degree, coefficient in enumerate(coefficients) + ) + + +class EntireRepresenterTests(unittest.TestCase): + def test_kernel_and_ungauged_value_normalization(self) -> None: + anchor = 0.7 - 0.35j + scale = 1.8 + self.assertAlmostEqual( + fock_kernel(anchor, anchor, scale).real, + math.exp(abs(anchor) ** 2 / scale**2), + delta=1.0e-14, + ) + direction = fock_value_direction(anchor, anchor, scale, gauge="none") + self.assertAlmostEqual(direction.real, 1.0, delta=1.0e-14) + self.assertAlmostEqual(direction.imag, 0.0, delta=1.0e-14) + self.assertAlmostEqual( + minimum_norm_squared( + anchor=anchor, scale=scale, kind="value", gauge="none" + ), + math.exp(-abs(anchor) ** 2 / scale**2), + delta=1.0e-14, + ) + + def test_zero_at_origin_gauge_uses_projected_kernel(self) -> None: + anchor = -0.45 + 0.8j + scale = 2.1 + self.assertAlmostEqual( + abs( + fock_value_direction( + 0j, anchor, scale, gauge="zero_at_origin" + ) + ), + 0.0, + delta=1.0e-14, + ) + self.assertAlmostEqual( + abs( + fock_value_direction( + anchor, anchor, scale, gauge="zero_at_origin" + ) + - 1.0 + ), + 0.0, + delta=1.0e-14, + ) + + def test_origin_value_is_incompatible_with_zero_gauge(self) -> None: + with self.assertRaisesRegex(ValueError, "incompatible"): + fock_value_direction(0j, 0j, 1.0, gauge="zero_at_origin") + + def test_derivative_representer_has_prescribed_local_derivative(self) -> None: + anchor = 0.3 + 0.4j + scale = 1.7 + step = 1.0e-6 + derivative = ( + fock_derivative_direction(anchor + step, anchor, scale) + - fock_derivative_direction(anchor - step, anchor, scale) + ) / (2.0 * step) + self.assertAlmostEqual(derivative.real, 1.0, delta=2.0e-10) + self.assertAlmostEqual(derivative.imag, 0.0, delta=2.0e-10) + self.assertEqual(fock_derivative_direction(0j, anchor, scale), 0j) + self.assertEqual(fock_derivative_direction(2.0 - 3.0j, 0j, scale), 2.0 - 3.0j) + + def test_minimum_norm_bounds_hold_for_polynomial_competitors(self) -> None: + anchor = 0.6 - 0.2j + scale = 1.4 + + ungauged_minimum = minimum_norm_squared( + anchor=anchor, scale=scale, kind="value", gauge="none" + ) + for slope in (0j, 0.2 + 0.4j, -0.7 + 0.1j): + # p(z) = 1 + slope * (z - anchor), hence p(anchor) = 1. + coefficients = [1.0 - slope * anchor, slope] + self.assertGreaterEqual( + polynomial_fock_norm_squared(coefficients, scale) + 1.0e-14, + ungauged_minimum, + ) + + gauged_minimum = minimum_norm_squared( + anchor=anchor, scale=scale, kind="value", gauge="zero_at_origin" + ) + # p(z) = z/anchor satisfies p(0)=0 and p(anchor)=1. + gauged_competitor = polynomial_fock_norm_squared([0j, 1.0 / anchor], scale) + self.assertGreaterEqual(gauged_competitor + 1.0e-14, gauged_minimum) + + def test_local_data_and_finite_superposition(self) -> None: + value_datum = LocalDatum( + anchor=0.8 + 0.1j, + amplitude=0.03 - 0.05j, + scale=1.5, + kind="value", + gauge="none", + ) + self.assertAlmostEqual( + abs(evaluate_q(value_datum.anchor, [value_datum]) - value_datum.amplitude), + 0.0, + delta=1.0e-14, + ) + + derivative_datum = LocalDatum( + anchor=-0.25 + 0.5j, + amplitude=-0.02 + 0.04j, + scale=2.0, + kind="derivative", + gauge="zero_at_origin", + ) + step = 1.0e-6 + derivative = ( + evaluate_q(derivative_datum.anchor + step, [derivative_datum]) + - evaluate_q(derivative_datum.anchor - step, [derivative_datum]) + ) / (2.0 * step) + self.assertAlmostEqual( + abs(derivative - derivative_datum.amplitude), 0.0, delta=2.0e-11 + ) + self.assertEqual(evaluate_q(0j, [derivative_datum]), 0j) + + z = -0.4 + 0.9j + self.assertEqual( + evaluate_q(z, [value_datum, derivative_datum]), + value_datum.contribution(z) + derivative_datum.contribution(z), + ) + + def test_complete_state_matches_direct_complex_evaluation(self) -> None: + data = [ + LocalDatum( + anchor=0.5 - 0.25j, + amplitude=0.08 + 0.03j, + scale=1.6, + kind="value", + gauge="zero_at_origin", + ), + LocalDatum( + anchor=-0.2 + 0.4j, + amplitude=-0.01 + 0.02j, + scale=1.2, + kind="derivative", + gauge="zero_at_origin", + ), + ] + self.assertEqual(evaluate_q(0j, data), 0j) + z = 0.3 + 0.7j + state = evaluate_regular_state( + z, + data, + zeros=[DivisorPoint(-0.5 + 0.2j, 2)], + poles=[DivisorPoint(1.1 - 0.4j)], + gain=0.8 - 0.3j, + ) + + self.assertAlmostEqual(abs(state.f - state.r * cmath.exp(state.q)), 0.0, delta=1.0e-14) + self.assertAlmostEqual(state.log_modulus, math.log(abs(state.f)), delta=1.0e-14) + phase_error = cmath.phase(cmath.exp(1j * (cmath.phase(state.f) - state.phase))) + self.assertAlmostEqual(phase_error, 0.0, delta=1.0e-14) + self.assertEqual(state.re_q, state.q.real) + self.assertEqual(state.im_q, state.q.imag) + + def test_entire_factor_is_nonzero_on_reference_grid(self) -> None: + data = [ + LocalDatum(0.7 + 0.2j, 0.1j, 1.3, "value", "none"), + LocalDatum(-0.4 + 0.5j, 0.03, 1.8, "derivative", "zero_at_origin"), + ] + for real in range(-3, 4): + for imaginary in range(-3, 4): + q = evaluate_q(complex(real, imaginary) / 2.0, data) + self.assertNotEqual(cmath.exp(q), 0j) + + def test_scale_and_multiplicity_validation(self) -> None: + with self.assertRaisesRegex(ValueError, "positive"): + fock_kernel(0j, 0j, 0.0) + with self.assertRaisesRegex(ValueError, "multiplicity"): + DivisorPoint(0j, 0) + with self.assertRaisesRegex(ValueError, "multiplicity"): + DivisorPoint(0j, 1.5) # type: ignore[arg-type] + with self.assertRaisesRegex(ValueError, "gain"): + evaluate_regular_state(1j, [], gain=0j) + + +if __name__ == "__main__": + unittest.main()