PF ideals and their power series solutions

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).",

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