diff --git a/sage/period_computation.sage b/sage/period_computation.sage index 9ff6458..56eb473 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, var +from sage.all import sage_eval, PolynomialRing, QQ, SR, log, matrix, prod, var load("sage/util.py") @@ -104,6 +104,145 @@ class PFOperatorAnsatz(PFOperator): super().__init__(operator_string=operator_string, no_variables=no_variables, extra_locals=extra_locals) +class PFIdeal: + """ + A class representing a Picard-Fuchs ideal, which is a collection of PFOperators. + """ + def __init__(self, operators: list): + if not all(op.no_variables == operators[0].no_variables for op in operators): + raise ValueError("All operators must have the same number of variables.") + self.no_variables = operators[0].no_variables + self.operators = operators + logger.info("Initialised PFIdeal with %d operator(s).", len(self.operators)) + + def add_operator(self, pf_operator: PFOperator): + self.operators.append(pf_operator) + logger.info("Added PFOperator to PFIdeal: %s", pf_operator.operator_string) + + def remove_operator(self, pf_operator: PFOperator): + self.operators.remove(pf_operator) + logger.info("Removed PFOperator from PFIdeal: %s", pf_operator.operator_string) + + 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. + + A term c*z^p*theta^q sends a_k*z^k (true exponent k+indicials) to + c*(k+indicials)^q*a_k*z^(k+p): it shifts index k up by p (p>=0), never down or sideways. + E.g. z0*theta0 sends a_k*z0^k to (k+rho0)*a_k*z0^(k+1). + + This is the multivariate Frobenius method: canonical series solutions of a regular + holonomic D-ideal via its indicial ideal. See M. Saito, B. Sturmfels, N. Takayama, + "Gröbner Deformations of Hypergeometric Differential Equations", Algorithms and + Computation in Mathematics vol. 6, Springer, 2000, chs. 2-3. + """ + if len(indicials) != self.no_variables: + raise ValueError( + "Indicials must have length no_variables=%d, got %d." % (self.no_variables, len(indicials)) + ) + indicials = [QQ(rho) for rho in indicials] + + operator_terms = [ + [ + (monomial.exponents()[0][:self.no_variables], monomial.exponents()[0][self.no_variables:], coeff) + for coeff, monomial in pf_operator.operator + ] + for pf_operator in self.operators + ] + + def eigenvalue(k, q): + # theta_i^q_i acts on z_i^(k_i + indicials[i]) as multiplication by + # (k_i + indicials[i])^q_i; theta^q's combined eigenvalue is the product over i. + value = QQ(1) + for i in range(self.no_variables): + if q[i]: + value *= (k[i] + indicials[i]) ** q[i] + return value + + # Process multi-indices in order of increasing total degree: since every operator + # monomial has p >= 0 (componentwise), the z^m coefficient of L(y) only ever + # depends on a_k for k <= m, so by this point every k < m has already been solved. + z_indices = sorted(_multi_indices(self.no_variables, order), key=sum) + + # solved[k] holds a_k written as a vector of coefficients over the `dimension` + # independent solutions found so far (a basis of the solution space up to k). + dimension = 0 + solved = {} + + for m in z_indices: + # For each operator, "coefficient of z^m in L(y) = 0" splits into a diagonal + # part (the p=0, theta-only monomials, whose unknown is a_m itself) plus a + # known part contributed by already-solved a_k with k = m - p, p > 0. + diagonals = [] + known_contributions = [] + for terms in operator_terms: + diagonal = QQ(0) + contribution = [QQ(0)] * dimension + for p, q, coeff in terms: + k = tuple(m[i] - p[i] for i in range(self.no_variables)) + if any(ki < 0 for ki in k): + continue # this monomial would need a_k for a negative multi-index k: no such term + ev = eigenvalue(k, q) + if ev == 0: + continue # theta^q kills z_i^(k_i + indicials[i]) here, so this monomial contributes nothing + if k == m: + diagonal += coeff * ev # p = 0: coefficient multiplying the still-unknown a_m + else: + k_vector = solved[k] # p > 0: a_k is already known, add its contribution + for i in range(dimension): + contribution[i] += coeff * ev * k_vector[i] + diagonals.append(diagonal) + known_contributions.append(contribution) + + # An operator with a nonzero diagonal lets us solve a_m = -(known part)/diagonal + # directly; this is exactly the indicial equation being nonzero at m + indicials. + active = next((r for r, d in enumerate(diagonals) if d != 0), None) + + if active is not None: + value = [-known_contributions[active][i] / diagonals[active] for i in range(dimension)] + # Every operator's equation at m must independently be satisfied by this + # same a_m; disagreement means the ideal is inconsistent with these indicials. + for r, d in enumerate(diagonals): + if any(d * value[i] + known_contributions[r][i] != 0 for i in range(dimension)): + raise ValueError( + "Inconsistent Picard-Fuchs ideal or indicial exponents %s at multidegree %s." + % (indicials, m) + ) + solved[m] = value + else: + # Resonance: every operator's indicial part vanishes at m + indicials, so a_m + # cannot be pinned down by this equation. If the already-known lower-degree + # data still forces a nonzero constraint here, satisfying it would require a + # log(z)-term solution, which this method (deliberately) does not compute. + if any(x != 0 for contribution in known_contributions for x in contribution): + raise ValueError( + "Resonance at multidegree %s for indicials %s would require a logarithmic " + "solution, which find_power_series_solution does not compute." % (m, indicials) + ) + # Otherwise a_m is genuinely free: it starts a new independent solution, so + # extend every previously solved coefficient with a 0 in this new direction. + dimension += 1 + for v in solved.values(): + v.append(QQ(0)) + solved[m] = [QQ(0)] * (dimension - 1) + [QQ(1)] + + log_index = tuple([0] * self.no_variables) + solutions = [ + Period( + no_variables=self.no_variables, + coefficients={log_index: {m: solved[m][i] for m in z_indices}}, + order=order, + indicials=indicials, + ) + for i in range(dimension) + ] + + logger.info( + "Found %d power series solution(s) for indicials %s at order %d.", len(solutions), indicials, order + ) + return solutions + + class Period: """ A class representing a period as a formal power series in z-variables and their logs. @@ -166,7 +305,18 @@ class Period: expression += monomial log_substitutions = {log_gens[i]: log(SR(z_gens[i])) for i in range(self.no_variables)} - return str(SR(expression).subs(log_substitutions)) + result = SR(expression).subs(log_substitutions) + + # A nonzero indicial ρ_i means the coefficients above are for z_i^k, but the + # actual solution is z_i^(ρ_i + k); make that explicit in the printed form. + indicial_prefactor = prod( + (SR(z_gens[i]) ** self.indicials[i] for i in range(self.no_variables) if self.indicials[i] != 0), + SR(1), + ) + if indicial_prefactor != 1: + result *= indicial_prefactor + + return str(result) def _max_z_degree(self) -> int: # Highest total z-degree (sum of the z-multi-index) among all coefficients. @@ -194,11 +344,16 @@ class Period: shape as self.coefficients. It uses the product rule theta_i(z^a log(z)^k) = a_i * z^a log(z)^k + k_i * z^a log(z)^(k - e_i). + + z_index entries are offsets from self.indicials: a stored z_index of a really means + z^(a + self.indicials[index]), so theta_i's eigenvalue is a_i + self.indicials[index] + rather than the bare a_i (self.indicials is all-zero unless explicitly given, in which + case this reduces to the ordinary power-series rule). """ result = {} for log_index, z_dict in coefficients.items(): for z_index, coeff in z_dict.items(): - a_i = z_index[index] + a_i = z_index[index] + self.indicials[index] if a_i != 0: inner = result.setdefault(log_index, {}) inner[z_index] = inner.get(z_index, 0) + a_i * coeff @@ -250,6 +405,7 @@ class Period: no_variables=self.no_variables, coefficients=result_coefficients, order=self.order, + indicials=self.indicials, ) def find_annihilating_operators(self, z_degree: int, theta_degree: int) -> list: @@ -277,8 +433,18 @@ class Period: return operators - def __init__(self, no_variables: int = 1, coefficients: dict = None, period_string: str = None, order: int = None): + def __init__( + self, + no_variables: int = 1, + coefficients: dict = None, + period_string: str = None, + order: int = None, + indicials: list = None, + ): self.no_variables = no_variables + # indicials[i] is the (rational) Frobenius exponent ρ_i of z_i: a stored + # z-index of a really represents z^(a + indicials[i]). + self.indicials = list(indicials) if indicials is not None else [0] * no_variables # Use coefficients or period_string to initialize the period if period_string is not None: @@ -311,7 +477,7 @@ class PeriodAnsatz(Period): A class for period Ansätze, characterised by number of variables and their z- and log-multi-degrees. """ - def __init__(self, no_variables: int, z_degree: int, log_degree: int): + def __init__(self, no_variables: int, z_degree: int, log_degree: int, indicials: list = None): self.no_variables = no_variables self.z_degree = z_degree self.log_degree = log_degree @@ -329,7 +495,7 @@ class PeriodAnsatz(Period): for log_index in log_indices } - super().__init__(no_variables=no_variables, coefficients=coefficients, order=z_degree) + super().__init__(no_variables=no_variables, coefficients=coefficients, order=z_degree, indicials=indicials) self.expansion_coefficients = self.coefficients logger.info( "Initialised PeriodAnsatz with %d unknown coefficient(s).", diff --git a/tests/test_period_computation.py b/tests/test_period_computation.py index b1ab3d4..efa942e 100644 --- a/tests/test_period_computation.py +++ b/tests/test_period_computation.py @@ -99,3 +99,75 @@ def test_simplify_factorises_theta_polynomial_per_z_monomial(): # The underlying (expanded) operator is unchanged - only the display string differs. assert simplified.operator == op.operator assert simplified.operator_string == "-5*(5*theta0 + 4)*(5*theta0 + 3)*(5*theta0 + 2)*(5*theta0 + 1)*z0 + theta0^4" + + +def test_find_power_series_solution_recovers_quintic_period(): + # The quintic's Picard-Fuchs operator has indicial equation theta0^4 = 0 at z0 = 0 + # (a quadruple root at 0), so its holomorphic (non-logarithmic) power series solution + # is found at indicial 0. Up to normalisation this is the classic quintic period + # 1 + 120*z0 + 113400*z0^2 + ... + op = PFOperator( + "theta0^4 - 5*z0*(5*theta0 + 1)*(5*theta0 + 2)*(5*theta0 + 3)*(5*theta0 + 4)", + no_variables=1, + ) + ideal = PFIdeal([op]) + + solutions = ideal.find_power_series_solution(indicials=[0], order=3) + + assert len(solutions) == 1 + assert solutions[0].coefficients == {(0,): {(0,): 1, (1,): 120, (2,): 113400, (3,): 168168000}} + + +def test_find_power_series_solution_handles_rational_indicial(): + # theta0 - 1/2 kills z0^(1/2 + k) only when k = 0, since its eigenvalue is 1/2 + k; + # so at indicial 1/2 there is exactly one solution (a constant multiple of sqrt(z0)), + # while at indicial 0 no power series solution exists at all. + op = PFOperator("theta0 - 1/2", no_variables=1) + ideal = PFIdeal([op]) + + solutions = ideal.find_power_series_solution(indicials=[1/2], order=3) + assert len(solutions) == 1 + assert solutions[0].coefficients == {(0,): {(0,): 1, (1,): 0, (2,): 0, (3,): 0}} + assert solutions[0].period_string == "sqrt(z0)" + + assert ideal.find_power_series_solution(indicials=[0], order=3) == [] + + +def test_ideal_finds_holomorphic_solution_and_recovers_first_operator(): + # Ideal generated by two operators (paper's z1, z2, theta1, theta2 <-> code's z0, z1, + # theta0, theta1): + # 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", + 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", + no_variables=2, + ) + ideal = PFIdeal([M1, M2]) + + solutions = ideal.find_power_series_solution(indicials=[0, QQ(1) / 2], order=6) + assert len(solutions) == 1 + + period = solutions[0] + assert period.apply_operator(M1).coefficients == {} + assert period.apply_operator(M2).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 + + # 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