From 82765b201c46e9595baa8364007b2e8c397387f7 Mon Sep 17 00:00:00 2001 From: Julian Piribauer Date: Wed, 9 Sep 2026 21:53:58 +0200 Subject: [PATCH] Implementing moving origin of ideal --- sage/period_computation.sage | 147 ++++++++++++++++++++++++++++++- tests/test_period_computation.py | 55 +++++++----- 2 files changed, 177 insertions(+), 25 deletions(-) diff --git a/sage/period_computation.sage b/sage/period_computation.sage index 70d5a65..44b8530 100644 --- a/sage/period_computation.sage +++ b/sage/period_computation.sage @@ -1,6 +1,6 @@ import copy import logging -from sage.all import sage_eval, PolynomialRing, QQ, SR, log, matrix, prod, var +from sage.all import sage_eval, PolynomialRing, QQ, SR, function, log, matrix, prod, solve, var load("sage/util.py") @@ -71,6 +71,139 @@ class PFOperator: logger.debug("Simplified PFOperator to: %s", simplified.operator_string) return simplified + def change_coordinates(self, *new_coordinates: str) -> "PFOperator": + """ + Rewrites this operator in new coordinates w0, ..., w{n-1} := new_coordinates(z0, ..., + z{n-1}), where new_coordinates is given as n expressions in the *old* coordinates + z0, ..., z{n-1}; the inverse change of coordinates needed to do this is solved for + automatically (the two aren't independent - one determines the other). The new + origin w = 0 is the point z0, ..., z{n-1} at which new_coordinates vanishes. + + E.g. pf_operator.change_coordinates("1 - 5*z0", "1/z1") moves the origin to + z0 = 1/5, z1 = infinity. + """ + n = self.no_variables + if len(new_coordinates) != n: + raise ValueError( + "Expected %d new coordinate expression(s), got %d." % (n, len(new_coordinates)) + ) + + z_gens = [var('z%d' % i) for i in range(n)] + parse_locals = {str(z): z for z in z_gens} + try: + new_coords = [SR(sage_eval(expr, locals=parse_locals)) for expr in new_coordinates] + except Exception as e: + raise ValueError("Invalid coordinate expression: %s" % e) + + w_gens = [var('w%d' % i) for i in range(n)] + equations = [w_gens[i] == new_coords[i] for i in range(n)] + solutions = solve(equations, z_gens, solution_dict=True) + if not solutions: + raise ValueError("Could not invert the given change of coordinates %s." % (new_coordinates,)) + inverse_coords = [solutions[0][z] for z in z_gens] + + Y = function('Y')(*w_gens) + y = Y.subs({w_gens[i]: new_coords[i] for i in range(n)}) + + def apply_theta(expr, i): + return z_gens[i] * expr.diff(z_gens[i]) + + total = SR(0) + for coeff, monomial in self.operator: + exponents = monomial.exponents()[0] + z_exponents = exponents[:n] + theta_exponents = exponents[n:] + + term = y + for i, power in enumerate(theta_exponents): + for _ in range(power): + term = apply_theta(term, i) + + z_monomial = prod(z_gens[i] ** z_exponents[i] for i in range(n)) + total += SR(coeff) * z_monomial * term + + total = total.subs({z_gens[i]: inverse_coords[i] for i in range(n)}) + total = total.simplify_full().expand() + + total_theta_degree = max(sum(monomial.exponents()[0][n:]) for _, monomial in self.operator) + theta_ring = PolynomialRing(QQ, n, ['theta%d' % i for i in range(n)]) + theta_gens = theta_ring.gens() + + # Peel off the coefficient of each derivative order k, highest total order first (so + # a coefficient can never still contain an as-yet-unextracted higher derivative as a + # factor), and rewrite w^k*D^k Y as a falling-factorial polynomial in theta applied + # to Y - the Euler-operator identity theta*(theta-1)*...*(theta-k+1) = w^k*d^k/dw^k, + # taken variable by variable. Each term's coefficient is left as a general rational + # function of w for now; only once every term has been collected is the whole + # operator scaled by their common denominator (see below). + remainder = total + terms = [] + for k in sorted(_multi_indices(n, total_theta_degree), key=lambda k: -sum(k)): + if sum(k) == 0: + dterm = Y + else: + args = [arg for i in range(n) if k[i] for arg in (w_gens[i], k[i])] + dterm = Y.diff(*args) + + coeff = remainder.coefficient(dterm) + if coeff == 0: + continue + remainder -= coeff * dterm + + falling_factorial = theta_ring(1) + for i in range(n): + for j in range(k[i]): + falling_factorial *= (theta_gens[i] - j) + + w_monomial = prod(w_gens[i] ** k[i] for i in range(n)) + terms.append(((coeff / w_monomial).simplify_rational(), falling_factorial)) + + remainder = remainder.simplify_full() + if remainder != 0: + raise ValueError("Residual nonzero after change of coordinates: %s" % remainder) + if not terms: + raise ValueError("Change of coordinates produced the zero operator.") + + w_ring = PolynomialRing(QQ, n, ['w%d' % i for i in range(n)]) + common_denominator = w_ring(1) + for rational_coeff, _ in terms: + common_denominator = common_denominator.lcm(w_ring(rational_coeff.denominator())) + + result_terms = {} + for rational_coeff, falling_factorial in terms: + scaled = w_ring((rational_coeff * SR(common_denominator)).simplify_rational()) + for w_coeff, w_monomial in scaled: + p = w_monomial.exponents()[0] + for theta_coeff, theta_monomial in falling_factorial: + q = theta_monomial.exponents()[0] + key = (p, q) + result_terms[key] = result_terms.get(key, QQ(0)) + w_coeff * theta_coeff + + content = gcd([QQ(coeff) for coeff in result_terms.values() if coeff != 0]) + if content not in (0, 1): + result_terms = {key: coeff / content for key, coeff in result_terms.items()} + + def _format_monomial(coeff, p, q): + factors = [str(QQ(coeff))] + for i in range(n): + if p[i]: + factors.append("z%d^%d" % (i, p[i]) if p[i] != 1 else "z%d" % i) + for i in range(n): + if q[i]: + factors.append("theta%d^%d" % (i, q[i]) if q[i] != 1 else "theta%d" % i) + return "*".join(factors) + + terms_string = [ + _format_monomial(coeff, p, q) for (p, q), coeff in result_terms.items() if coeff != 0 + ] + if not terms_string: + raise ValueError("Change of coordinates produced the zero operator.") + operator_string = " + ".join(terms_string) + + new_operator = PFOperator(operator_string, no_variables=n) + logger.debug("Changed coordinates of PFOperator to: %s", new_operator.operator_string) + return new_operator + class PFOperatorAnsatz(PFOperator): """ A class for Picard-Fuchs operator Ansätze, characterised by number of variables and their z- and theta-multi-degrees. @@ -122,7 +255,17 @@ class PFIdeal: def remove_operator(self, pf_operator: PFOperator): self.operators.remove(pf_operator) logger.info("Removed PFOperator from PFIdeal: %s", pf_operator.operator_string) - + + def change_coordinates(self, *new_coordinates: str) -> "PFIdeal": + """ + Moves the origin of the whole system of differential equations to a new point, by + returning a new PFIdeal whose operators are each rewritten in new_coordinates; see + PFOperator.change_coordinates for what new_coordinates means. + """ + new_operators = [pf_operator.change_coordinates(*new_coordinates) for pf_operator in self.operators] + logger.info("Changed coordinates of PFIdeal with %d operator(s) to %s.", len(new_operators), new_coordinates) + return PFIdeal(new_operators) + def find_power_series_solution(self, indicials: list, order: int) -> list: """ Finds power series solutions (no logs) to this PFIdeal at given indicial exponents/order. diff --git a/tests/test_period_computation.py b/tests/test_period_computation.py index 81a2270..277a701 100644 --- a/tests/test_period_computation.py +++ b/tests/test_period_computation.py @@ -125,42 +125,39 @@ def test_find_power_series_solution_handles_rational_indicial(): def test_ideal_finds_power_series_solution_and_recovers_first_operator(): - # Picard--Fuchs ideal for P_{22211}[8]: - # M1 = theta2*(-2*theta1+2*theta2-1) + 2*(theta1-2*theta2-1)*(theta1-2*theta2)*z2 - # M2 = theta1^2*(2*(theta1-2*theta2)*z2-theta2) - 16*(2*theta1+1)*(4*theta1+1)*(4*theta1+3)*z1*z2 - M1 = PFOperator( - "theta1*(-2*theta0 + 2*theta1 - 1) + 2*(theta0 - 2*theta1 - 1)*(theta0 - 2*theta1)*z1", + # Picard--Fuchs ideal for P_{22211}[8] at the MUM point: + # D1 = Theta_x^2*(Theta_x - 2*Theta_y) - 4*x*(4*Theta_x + 3)*(4*Theta_x + 2)*(4*Theta_x + 1) + # D2 = Theta_y^2 - y*(2*Theta_y - Theta_x + 1)*(2*Theta_y - Theta_x) + # + # We test moving the ideal to the intersection with the locus of Strong Coupling. + D1 = PFOperator( + "theta0^2*(theta0 - 2*theta1) - 4*z0*(4*theta0 + 3)*(4*theta0 + 2)*(4*theta0 + 1)", no_variables=2, ) - M2 = PFOperator( - "theta0^2*(2*(theta0 - 2*theta1)*z1 - theta1) - 16*(2*theta0 + 1)*(4*theta0 + 1)*(4*theta0 + 3)*z0*z1", + D2 = PFOperator( + "theta1^2 - z1*(2*theta1 - theta0 + 1)*(2*theta1 - theta0)", no_variables=2, ) - ideal = PFIdeal([M1, M2]) + moved = PFIdeal([D1, D2]).change_coordinates("z0", "z1 - 1/4") + op1, op2 = moved.operators - solutions = ideal.find_power_series_solution(indicials=[0, QQ(1) / 2], order=6) + solutions = moved.find_power_series_solution(indicials=[0, QQ(1) / 2], order=14) assert len(solutions) == 1 period = solutions[0] - assert period.apply_operator(M1).coefficients == {} - assert period.apply_operator(M2).coefficients == {} + assert period.apply_operator(op1).coefficients == {} + assert period.apply_operator(op2).coefficients == {} # Normalisation convention: the free parameter at the leading (0, 0) coefficient is 1. assert period.coefficients[(0, 0)][(0, 0)] == 1 - assert period.coefficients[(0, 0)][(1, 1)] == -32 + assert period.coefficients[(0, 0)][(1, 0)] == 0 - # M1 lives entirely within z_degree <= 1, theta_degree <= 2 (a single bare factor of - # z2, quadratic in theta), so this is the natural Ansatz level to look for it at. recovered = period.find_annihilating_operators(z_degree=1, theta_degree=2) assert len(recovered) == 1 - - # recovered[0] should be a scalar multiple of M1 - compare via the coefficient of the - # bare theta2 (theta1 in code) monomial, which is nonzero in M1. - ratio = ( - recovered[0].operator.monomial_coefficient(M1.theta_gens[1]) - / M1.operator.monomial_coefficient(M1.theta_gens[1]) - ) - assert ratio != 0 - assert recovered[0].operator == ratio * M1.operator + assert recovered[0].operator == PFOperator( + "-2*z1*theta0^2 + 8*z1*theta0*theta1 - 8*z1*theta1^2 + 2*z1*theta0 - 4*z1*theta1 " + "+ 2*theta0*theta1 - 2*theta1^2 + theta1", + no_variables=2, + ).operator def test_single_operator_recovers_itself_from_its_holomorphic_period(): @@ -192,3 +189,15 @@ def test_single_operator_recovers_itself_from_its_holomorphic_period(): ) assert ratio != 0 assert recovered[0].operator == ratio * L.operator + + +def test_change_coordinates_handles_non_monomial_denominator(): + op = PFOperator("theta0^2 - z0", no_variables=1) + + result = op.change_coordinates("z0/(1-z0)") + + expected = PFOperator( + "(1+z0)^3*theta0^2 + z0*(1+z0)^2*theta0 - z0", + no_variables=1, + ) + assert result.operator == expected.operator -- 2.54.0