Implementing Frobenius basis for PFIdeals

2026-11-01 16:46:13 +01:00
4 changed files with 895 additions and 110 deletions
@@ -8,7 +8,7 @@ It allows for computation of discriminant loci and topological data of manifolds
 Supported initialisations are represented by the following examples.
 ```python
-elliptic_curve_D = ToricPolytopeProjectiveSpace([1, 2, 3], model_name="elliptic_curve_D")
+elliptic_curve_D = ToricPolytopeProjectiveSpace([3, 2, 1], model_name="elliptic_curve_D")
 CY3_quintic = ToricPolytopeProjectiveSpace([1, 1, 1, 1, 1], model_name="quintic")
 CY3_bicubic = ToricPolytopeCICY([[3, 3]], model_name="bi-cubic")
@@ -150,3 +150,42 @@ recovered[0].simplify().operator_string
 `recovered[0]` is, up to scale, `M1` &mdash; recovered purely from `M1`, `M2`'s shared power series
 solution. 
 ### Frobenius bases
 A `PFIdeal` also computes its complete Frobenius basis &mdash; all solutions, logarithmic ones included &mdash;
 at the origin of its current coordinates; for another point, move it to the origin with `change_coordinates`
 first. The exponents are the `indicials`, the roots of the system's indicial ideal, which is completed by linear
 algebra where the operators' z-free parts alone don't determine it. The optional `hodge_numbers` &mdash; the
 numbers of solutions of top log-degree 0, 1, ..., n at a MUM point, e.g. `[1, 2, 2, 1]` for a two-parameter
 threefold &mdash; certify that completion through their sum, the rank, and bound the powers of logarithms. For
 P_{2,2,2,1,1}[8] at its MUM point and at the point (z0, z1) = (0, 1/4) on its strong coupling divisor:
 ```python
 D1 = PFOperator(
    "theta0^2*(theta0 - 2*theta1) - 4*z0*(4*theta0 + 3)*(4*theta0 + 2)*(4*theta0 + 1)",
    no_variables=2,
 )
 D2 = PFOperator("theta1^2 - z1*(2*theta1 - theta0 + 1)*(2*theta1 - theta0)", no_variables=2)
 ideal = PFIdeal([D1, D2], hodge_numbers=[1, 2, 2, 1])
 strong_coupling = ideal.change_coordinates("z0", "z1 - 1/4")
 ideal.indicials()
 strong_coupling.indicials()
 basis = strong_coupling.frobenius_basis(order=2)
 basis[5].period_string
 ```
 ```term
 [((0, 0), 6)]
 [((0, 0), 4), ((0, 1/2), 2)]
 INFO:__main__:Completed the indicial ideal at z = 0 at multiplier degree 2 with 1 extra operator(s): colength 6.
 INFO:__main__:Exponent class of (0, 0): 4 solution(s) with leading terms ['1', 'log(z0)', 'log(z0)^2', 'log(z0)^3'].
 INFO:__main__:Exponent class of (0, 1/2): 2 solution(s) with leading terms ['sqrt(z1)', 'sqrt(z1)*(log(z0) + log(z1))'].
 INFO:__main__:Found a Frobenius basis of 6 solution(s) at order 2.
 '-1/225*(7200*z0*z1*log(z0) - 720*z1^2*log(z0) + 7200*z0*z1*log(z1) - 720*z1^2*log(z1) + 33600*z0*z1 - 732*z1^2 + 300*z1*log(z0) + 300*z1*log(z1) + 250*z1 - 225*log(z0) - 225*log(z1))*sqrt(z1)'
 ```
 Within each exponent class the basis is canonical: every solution has its own leading term with coefficient 1,
 at which the class' other solutions vanish. Coefficients are those of plain powers of log(z_i), without factors
 of 2&pi;i.
@@ -8,20 +8,35 @@ PFOperator
    simplify                      Group by z-monomials and factorise theta-polynomials
    change_coordinates            Rewrite the operator in new coordinates
    _operator_from_string         Parse a string into the operator polynomial
    _theta_parts                  Theta-polynomial multiplying each z-monomial
 PFOperatorAnsatz(PFOperator)
    __init__                      Ansatz of given theta- and z-multi-degrees
 PFIdeal
-    __init__                      Construct from a list of PFOperators
+    __init__                      Construct from a list of PFOperators and optional Hodge numbers
    add_operator                  Append a PFOperator
    remove_operator               Remove a PFOperator
    change_coordinates            Rewrite all operators in new coordinates
    indicial_ideal                Indicial ideal at z = 0, completed by linear algebra
    indicials                     Indicials (exponents) at z = 0 with multiplicities
    frobenius_basis               Basis of all solutions at z = 0, logarithmic ones included
    find_power_series_solution    Power series solutions at given indicials/order
    _operator_theta_parts         Theta-parts of all operators in a common theta-ring
    _complete_indicial_ideal      Indicial ideal and the extra operators needed to see it
    _solve_exponent_class         Solutions with given indicials, up to a degree in log(z)
 Helpers
    _z_order_key                  Graded order on z-monomials (initial forms, solving order)
    _left_multiply                Normal-ordered product z^a * theta^c * operator
    _indicial_roots               Roots of an indicial ideal with their multiplicities
    _leading_term                 Lowest term of a solution, for logging
 """
 import copy
 import logging
-from sage.all import sage_eval, PolynomialRing, QQ, SR, function, prod, solve, var
+from sage.all import (
    sage_eval, PolynomialRing, QQ, SR, ZZ, binomial, block_matrix, floor, function, log, matrix, prod, solve, var
 )
 load("sage/util.py")
@@ -225,6 +240,24 @@ class PFOperator:
        logger.debug("Changed coordinates of PFOperator to: %s", new_operator.operator_string)
        return new_operator
    def _theta_parts(self, theta_ring) -> dict:
        """
        Splits this operator as sum_p z^p * P_p(theta) and returns {p: P_p}, each P_p an element of
        theta_ring = QQ[theta0, ..., theta{n-1}]; z-monomials whose theta-polynomial vanishes are dropped.
        theta_ring has to be built as PolynomialRing(QQ, n, names): without the explicit n, a single
        variable gives a univariate ring, whose exponents() and monomial() work differently.
        """
        if self.ring.base_ring() != QQ:
            raise ValueError("Only operators with rational coefficients can be split into theta-parts.")
        theta_gens = theta_ring.gens()
        parts = {}
        for coeff, monomial in self.operator:
            exponents = monomial.exponents()[0]
            z_exponents = tuple(exponents[:self.no_variables])
            theta_monomial = prod(theta_gens[i] ** exponents[self.no_variables + i] for i in range(self.no_variables))
            parts[z_exponents] = parts.get(z_exponents, theta_ring(0)) + coeff * theta_monomial
        return {p: P for p, P in parts.items() if P != 0}
 class PFOperatorAnsatz(PFOperator):
    """
    A class for Picard-Fuchs operator Ansätze, characterised by number of variables and their z- and theta-multi-degrees.
@@ -261,12 +294,22 @@ class PFOperatorAnsatz(PFOperator):
 class PFIdeal:
    """
    A class representing a Picard-Fuchs ideal, which is a collection of PFOperators.
    hodge_numbers optionally lists the numbers of solutions of top log-degree 0, 1, ..., n at a point of
    maximal unipotent monodromy, i.e. h^{n,0}, ..., h^{0,n} of the (horizontal part of the) middle
    cohomology - e.g. [1, 2, 2, 1] for a two-parameter Calabi-Yau threefold. Their sum is the holonomic
    rank, the number of solutions at every point, and n bounds the degree in log(z) of every solution at
    every point (monodromy theorem: N^(n+1) = 0 for a weight-n variation of Hodge structure). Without them,
    indicials and frobenius_basis still work, just without the certificates and checks these allow.
    """
-    def __init__(self, operators: list):
+    def __init__(self, operators: list, hodge_numbers: list = None):
        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.")
        if hodge_numbers is not None and (not hodge_numbers or any(h not in ZZ or h < 1 for h in hodge_numbers)):
            raise ValueError("hodge_numbers must be a non-empty list of positive integers, got %s." % (hodge_numbers,))
        self.no_variables = operators[0].no_variables
        self.operators = operators
        self.hodge_numbers = list(hodge_numbers) if hodge_numbers is not None else None
        logger.info("Initialised PFIdeal with %d operator(s).", len(self.operators))
    def add_operator(self, pf_operator: PFOperator):
@@ -285,123 +328,484 @@ class PFIdeal:
        """
        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)
+        return PFIdeal(new_operators, hodge_numbers=self.hodge_numbers)
    def indicial_ideal(self, max_multiplier_degree: int = 4):
        """
        Returns the indicial ideal of this PFIdeal at z = 0, an ideal of QQ[theta0, ..., theta{n-1}] whose
        roots are the indicials (exponents) of the solutions there; see _complete_indicial_ideal.
        """
        indicial_ideal, _ = self._complete_indicial_ideal(max_multiplier_degree)
        return indicial_ideal
    def indicials(self, max_multiplier_degree: int = 4) -> list:
        """
        Returns the indicials (exponents) of this PFIdeal at z = 0, i.e. the roots of its indicial ideal, as
        a sorted list of (rho, multiplicity) pairs, the multiplicities adding up to the holonomic rank. E.g.
        the quintic has [((0,), 4)] at its MUM point and [((0,), 1), ((1,), 2), ((2,), 1)] at its conifold point.
        """
        return _indicial_roots(self.indicial_ideal(max_multiplier_degree))
    def frobenius_basis(self, order: int, max_log_degree: int = None, max_multiplier_degree: int = 4) -> list:
        """
        Computes a basis of all solutions of this PFIdeal at z = 0, logarithmic ones included, as Periods
        to the given order. For the basis at another point, move that point to the origin first with
        change_coordinates.
        The solutions are found one exponent class at a time: the indicials are grouped by their residue
        mod Z^n, and each class is solved once, from its componentwise smallest exponent rho_min (which
        becomes the Periods' indicials), so that solutions starting at rho_min + k with k > 0 - e.g. at the
        indicials 0, 1, 1, 2 of a conifold point - come out of the same solve. Their degree in log(z) is at
        most max_log_degree if given, and otherwise at most the class' total multiplicity minus 1 (the
        logarithms of the monodromies' unipotent parts act on the class' solutions as commuting nilpotent
        maps, so any product of that many of them vanishes) and at most n = len(hodge_numbers) - 1.
        Within each class the basis is canonical (see _solve_exponent_class): every solution has its own
        leading term with coefficient 1, at which the class' other solutions vanish. Coefficients are those
        of plain log(z_i)^b, without factors of 2*pi*i or factorials.
        This assumes the origin is a normal-crossing point of the discriminant, where every solution is a
        sum of terms z^rho * (power series) * (polynomial in log z). Elsewhere fewer solutions than the
        rank turn up, which is logged as a warning - as is a log structure at a MUM point that differs from
        hodge_numbers.
        """
        n = self.no_variables
        indicial_ideal, operator_parts = self._complete_indicial_ideal(max_multiplier_degree)
        theta_ring = indicial_ideal.ring()
        roots = _indicial_roots(indicial_ideal)
        classes = {}
        for rho, multiplicity in roots:
            classes.setdefault(tuple(r - floor(r) for r in rho), []).append((rho, multiplicity))
        exponent_classes = sorted(
            ((tuple(min(rho[i] for rho, _ in members) for i in range(n)), members) for members in classes.values()),
            key=lambda exponent_class: exponent_class[0],
        )
        cy_dimension = len(self.hodge_numbers) - 1 if self.hodge_numbers is not None else None
        basis = []
        for rho_min, members in exponent_classes:
            class_size = sum(multiplicity for _, multiplicity in members)
            # A solution starting at the indicial rho only turns up once the solve reaches z^(rho - rho_min).
            required_order = max(sum(rho[i] - rho_min[i] for i in range(n)) for rho, _ in members)
            if order < required_order:
                raise ValueError(
                    "Order %d is too low to reach the indicials %s; it must be at least %d."
                    % (order, [rho for rho, _ in members], required_order)
                )
            if max_log_degree is not None:
                log_degree = max_log_degree
            elif cy_dimension is not None:
                log_degree = min(class_size - 1, cy_dimension)
            else:
                log_degree = class_size - 1
            solutions = self._solve_exponent_class(operator_parts, theta_ring, rho_min, order, log_degree)
            logger.info(
                "Exponent class of %s: %d solution(s) with leading terms %s.",
                rho_min, len(solutions), [_leading_term(coefficients, rho_min) for coefficients in solutions],
            )
            if len(solutions) != class_size:
                logger.warning(
                    "Found %d solution(s) in the exponent class of %s, whose indicials have total multiplicity "
                    "%d: z = 0 may not be a normal-crossing point of the discriminant, or the degree in log(z), "
                    "capped at %d, is too low.",
                    len(solutions), rho_min, class_size, log_degree,
                )
            basis.extend(
                Period(no_variables=n, coefficients=coefficients, order=order, indicials=list(rho_min))
                for coefficients in solutions
            )
        if self.hodge_numbers is not None:
            if len(basis) != sum(self.hodge_numbers):
                logger.warning(
                    "Found %d solution(s), but the rank sum(hodge_numbers) is %d.", len(basis), sum(self.hodge_numbers)
                )
            # At a point of maximal unipotent monodromy - all indicials 0, with a solution of top log-degree
            # n - the numbers of solutions of each top log-degree are the Hodge numbers. A single exponent
            # class at 0 isn't enough: at (0, -1/27) of P_{333111}[12], on the MUM divisor z0 = 0, the
            # indicials are (0, 0) and (0, 1) and there is a log(z0)^4, but no MUM point.
            top_log_degrees = [max(sum(log_index) for log_index in period.coefficients) for period in basis]
            if (
                len(roots) == 1
                and all(r == 0 for r in roots[0][0])
                and max(top_log_degrees, default=-1) == cy_dimension
            ):
                log_structure = [top_log_degrees.count(degree) for degree in range(cy_dimension + 1)]
                if log_structure != self.hodge_numbers:
                    logger.warning(
                        "The numbers of solutions per top log-degree %s at this MUM point differ from "
                        "hodge_numbers %s.", log_structure, self.hodge_numbers,
                    )
        logger.info("Found a Frobenius basis of %d solution(s) at order %d.", len(basis), order)
        return basis
    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.
+        Finds the power series solutions (no logs) to this PFIdeal at given indicial exponents/order, i.e.
        all solutions of the form z^indicials * (power series in z); frobenius_basis finds all solutions,
        logarithmic ones included.
-        A term c*z^p*theta^q sends a_k*z^k (true exponent k+indicials) to
+        This is _solve_exponent_class with no room for logarithms, on the operators as given (without
-        c*(k+indicials)^q*a_k*z^(k+p): it shifts index k up by p (p>=0), never down or sideways.
+        completing their indicial ideal). A combination of solutions that could only be continued with a
-        E.g. z0*theta0 sends a_k*z0^k to (k+rho0)*a_k*z0^(k+1).
+        log(z)-term, at a resonance, is not a power series solution and is dropped.
        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]
        theta_ring, operator_parts = self._operator_theta_parts()
        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(
+            Period(no_variables=self.no_variables, coefficients=coefficients, order=order, indicials=indicials)
-                no_variables=self.no_variables,
+            for coefficients in self._solve_exponent_class(operator_parts, theta_ring, indicials, order, 0)
                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
    def _operator_theta_parts(self) -> tuple:
        # The theta-ring QQ[theta0, ..., theta{n-1}] and, in it, the theta-parts (see
        # PFOperator._theta_parts) of every nonzero operator of this PFIdeal.
        theta_ring = PolynomialRing(QQ, self.no_variables, ['theta%d' % i for i in range(self.no_variables)])
        operator_parts = [parts for parts in (op._theta_parts(theta_ring) for op in self.operators) if parts]
        if not operator_parts:
            raise ValueError("The PFIdeal contains no nonzero operator.")
        return theta_ring, operator_parts
    def _complete_indicial_ideal(self, max_multiplier_degree: int = 4) -> tuple:
        """
        Computes the indicial ideal of this PFIdeal at z = 0, together with the operators needed to see it.
        The indicial ideal consists of the theta-polynomials Q for which some operator in the (left) ideal
        has initial part z^e * Q(theta), i.e. this as its part at the lowest z-monomial w.r.t. _z_order_key.
        Its roots are the indicials (exponents) of the solutions at z = 0, and for a regular holonomic system,
        as every Picard-Fuchs system is, its colength - the number of roots counted with multiplicity - is
        the holonomic rank. See M. Saito, B. Sturmfels, N. Takayama, "Gröbner Deformations of
        Hypergeometric Differential Equations", Algorithms and Computation in Mathematics vol. 6,
        Springer, 2000, ch. 2.
        The initial parts of the generators alone can generate a strictly smaller ideal. E.g. for
        P_{22211}[8] moved to (z0, z1) = (0, 1/4) they are theta0^2*theta1 and theta1*(2*theta1 - 2*theta0 - 1),
        which vanish on the whole line theta1 = 0. The missing elements are initial parts of combinations of
        the generators in which the generators' own initial parts cancel (S-pairs), and are found here by
        linear algebra: the left multiples z^a * theta^c * L of the generators L with |a| + |c| <= d span a
        finite-dimensional space (a Macaulay matrix), whose echelon form w.r.t. the order on z-monomials
        exposes every initial part it contains. d grows until the result is zero-dimensional and complete:
        of colength rank = sum(hodge_numbers) if those are known, and otherwise unchanged from d - 1.
        Returns (indicial ideal, operator parts): the theta-parts (see PFOperator._theta_parts) of the
        generators, followed by those of the few extra operators whose initial parts were needed - solving
        with all of them sees exactly the right starting terms.
        """
        n = self.no_variables
        theta_ring, generator_parts = self._operator_theta_parts()
        rank = sum(self.hodge_numbers) if self.hodge_numbers is not None else None
        generator_initials = [parts[min(parts, key=_z_order_key)] for parts in generator_parts]
        # A single ordinary differential operator needs no completion: initial parts are multiplicative, so
        # its indicial ideal is generated by its own initial part (Fuchs' indicial polynomial).
        single_ode = n == 1 and len(generator_parts) == 1
        previous = None
        for degree in range(max_multiplier_degree + 1):
            # One row per left multiple, one column per (z-monomial, theta-monomial), the z-monomials in
            # z-order, so that each row of the echelon form starts at its initial z-monomial.
            rows = [
                _left_multiply(parts, multiplier[:n], multiplier[n:], theta_ring)
                for parts in generator_parts
                for multiplier in _multi_indices(2 * n, degree)
            ]
            columns = sorted(
                {(p, tuple(q)) for row in rows for p, P in row.items() for q in P.dict()},
                key=lambda column: (_z_order_key(column[0]), column[1]),
            )
            column_position = {column: j for j, column in enumerate(columns)}
            # Dense on purpose, here and in _solve_exponent_class: Sage's echelon form of *sparse* rational
            # matrices can be wrong (seen with Sage 9.7 for a 1 x 2 matrix met when solving the quintic at
            # its orbifold point).
            macaulay = matrix(QQ, len(rows), len(columns))
            for i, row in enumerate(rows):
                for p, P in row.items():
                    for q, coeff in P.dict().items():
                        macaulay[i, column_position[(p, tuple(q))]] = coeff
            echelon = macaulay.echelon_form()
            # Keep, greedily, every echelon row whose initial part isn't in the ideal yet, then drop again
            # those made redundant by rows kept later: every extra operator costs a full pass of the
            # solver, and very few are needed (one for P_{22211}[8] at (0, 1/4)).
            extra_parts, extra_initials = [], []
            current = theta_ring.ideal(generator_initials)
            for i, pivot in enumerate(echelon.pivots()):
                row = echelon.row(i).dict()
                initial_z_index = columns[pivot][0]
                initial = sum(
                    coeff * theta_ring.monomial(*columns[j][1])
                    for j, coeff in row.items()
                    if columns[j][0] == initial_z_index
                )
                if current.reduce(initial) == 0:
                    continue
                parts = {}
                for j, coeff in row.items():
                    p, q = columns[j]
                    parts[p] = parts.get(p, theta_ring(0)) + coeff * theta_ring.monomial(*q)
                extra_parts.append(parts)
                extra_initials.append(initial)
                current = theta_ring.ideal(generator_initials + extra_initials)
            for i in reversed(range(len(extra_parts))):
                others = generator_initials + extra_initials[:i] + extra_initials[i + 1:]
                if theta_ring.ideal(others).reduce(extra_initials[i]) == 0:
                    del extra_parts[i], extra_initials[i]
            indicial_ideal = theta_ring.ideal(generator_initials + extra_initials)
            finite = indicial_ideal.dimension() <= 0
            colength = indicial_ideal.vector_space_dimension() if finite else None
            logger.debug(
                "Indicial ideal at multiplier degree %d with %d extra operator(s): %s.", degree, len(extra_parts),
                "colength %d" % colength if finite else "not zero-dimensional",
            )
            if finite and rank is not None:
                if colength == rank:
                    break
                if colength < rank:
                    raise ValueError(
                        "The indicial ideal at z = 0 has colength %d, less than the rank sum(hodge_numbers) = %d: "
                        "an operator is not in the Picard-Fuchs ideal, or hodge_numbers is wrong." % (colength, rank)
                    )
            elif finite and (single_ode or indicial_ideal == previous):
                if not single_ode:
                    logger.warning(
                        "Taking the indicial ideal at z = 0 as complete since it didn't change from multiplier "
                        "degree %d to %d; give the PFIdeal hodge_numbers to certify this by its rank.",
                        degree - 1, degree,
                    )
                break
            previous = indicial_ideal
        else:
            raise ValueError(
                "Could not complete the indicial ideal at z = 0 up to multiplier degree %d: the operators may not "
                "generate the whole Picard-Fuchs ideal (whose rank then exceeds sum(hodge_numbers)), or a higher "
                "max_multiplier_degree is needed." % max_multiplier_degree
            )
        logger.info(
            "Completed the indicial ideal at z = 0 at multiplier degree %d with %d extra operator(s): colength %d.",
            degree, len(extra_parts), colength,
        )
        return indicial_ideal, generator_parts + extra_parts
    def _solve_exponent_class(
        self, operator_parts: list, theta_ring, indicials, order: int, max_log_degree: int
    ) -> list:
        """
        Finds all solutions y = sum_k z^(indicials + k) * g_k(log z) of the operators given by their
        theta-parts (see PFOperator._theta_parts), to the given order: k runs over the multi-indices >= 0,
        and every g_k is a polynomial in L_i = log(z_i) of total degree at most max_log_degree. Returns one
        coefficients dictionary {log_index: {z_index: coeff}} (see Period) per solution of a canonical basis.
        Since theta_i acts on z^(indicials + k) * g(L) as indicials_i + k_i + d/dL_i, a term c*z^p*theta^q
        sends it to c*z^(indicials + k + p) * (indicials + k + d/dL)^q g(L): it shifts k up by p (p >= 0),
        never down or sideways. So, as in Frobenius' method, the solutions can be found z-monomial by
        z-monomial - only with every unknown a polynomial g_k rather than a number, on whose coefficients a
        theta-polynomial P acts by the matrix of P(indicials + k + d/dL). Logarithms, indicials differing by
        integers (resonance) and several operators constraining each other are then all handled by the
        kernels and consistency conditions of these small matrices, with no derivatives with respect to the
        exponent and no case analysis.
        The basis is canonical: in the reduced echelon form of all solutions, each has its own leading
        coefficient 1 at which all others vanish, so in particular no solution's tail contains another one's
        leading term - as for the canonical series solutions of M. Saito, B. Sturmfels, N. Takayama,
        "Gröbner Deformations of Hypergeometric Differential Equations", Algorithms and Computation in
        Mathematics vol. 6, Springer, 2000, ch. 2.
        """
        n = self.no_variables
        theta_gens = theta_ring.gens()
        indicials = [QQ(rho) for rho in indicials]
        # Monomials L^b spanning the polynomials of total degree <= max_log_degree in L, lowest degree first
        # (1, L0, L1, L0^2, L0*L1, ...); this order decides which coefficient of a solution is normalised.
        log_indices = sorted(_multi_indices(n, max_log_degree), key=lambda b: (sum(b), tuple(-x for x in b)))
        log_position = {log_index: i for i, log_index in enumerate(log_indices)}
        size = len(log_indices)
        # By Taylor's formula, P(s + d/dL) sends L^b to the sum over b' <= b of
        # binom(b, b') * (d^(b - b')P/dtheta^(b - b'))(s) * L^b'. With this pattern of entries and every
        # theta-polynomial's derivatives computed once, each matrix below costs only a few evaluations.
        action_pattern = [
            (
                log_position[lower],
                log_position[upper],
                tuple(upper[i] - lower[i] for i in range(n)),
                prod(binomial(upper[i], lower[i]) for i in range(n)),
            )
            for upper in log_indices
            for lower in log_indices
            if all(lower[i] <= upper[i] for i in range(n))
        ]
        def taylor_derivatives(P):
            derivatives = {}
            for alpha in _multi_indices(n, max_log_degree):
                derivative = P
                for i in range(n):
                    if alpha[i]:
                        derivative = derivative.derivative(theta_gens[i], alpha[i])
                if derivative != 0:
                    derivatives[alpha] = derivative
            return derivatives
        def action_matrix(derivatives, shift):
            values = {alpha: derivative(*shift) for alpha, derivative in derivatives.items()}
            entries = {}
            for row, column, alpha, binomial_factor in action_pattern:
                if values.get(alpha, 0) != 0:
                    entries[(row, column)] = binomial_factor * values[alpha]
            # matrix() would make a sparse matrix from a dictionary; see the comment on dense matrices in
            # _complete_indicial_ideal.
            return matrix(QQ, size, size, entries, sparse=False)
        operator_derivatives = [{p: taylor_derivatives(P) for p, P in parts.items()} for parts in operator_parts]
        supports = [sorted(parts, key=_z_order_key) for parts in operator_parts]
        z_indices = sorted(_multi_indices(n, order), key=_z_order_key)
        # solved[k] holds g_k as a size x dimension matrix, whose column j holds g_k's coefficients
        # (w.r.t. log_indices) in the j-th of the `dimension` independent solutions found so far.
        dimension = 0
        solved = {}
        for k in z_indices:
            # Collect the equations whose latest unknown is g_k. The coefficient of z^(indicials + m) in L(y)
            # involves g_(m - p) for all p <= m (componentwise) in L's support, the latest of them in z-order
            # being g_(m - p') for the smallest such p'. Taking m = k + p' for those p' that are the smallest
            # element p <= k + p' of the support thus visits every equation exactly once - also for operators
            # whose lowest z-monomial isn't 1, like the extra ones from _complete_indicial_ideal.
            diagonal_blocks, known_blocks = [], []
            for derivatives, support in zip(operator_derivatives, supports):
                for position, pivot_exponent in enumerate(support):
                    m = tuple(k[i] + pivot_exponent[i] for i in range(n))
                    if any(all(p[i] <= m[i] for i in range(n)) for p in support[:position]):
                        continue
                    diagonal_blocks.append(
                        action_matrix(derivatives[pivot_exponent], [indicials[i] + k[i] for i in range(n)])
                    )
                    known = matrix(QQ, size, dimension)
                    if dimension:
                        for p in support[position + 1:]:
                            if all(p[i] <= m[i] for i in range(n)):
                                j = tuple(m[i] - p[i] for i in range(n))  # earlier than k in z-order: solved
                                known += action_matrix(derivatives[p], [indicials[i] + j[i] for i in range(n)]) * solved[j]
                    known_blocks.append(known)
            diagonal = block_matrix(QQ, [[block] for block in diagonal_blocks], subdivide=False)
            known = block_matrix(QQ, [[block] for block in known_blocks], subdivide=False)
            # Row-reduce [diagonal | known]: rows with a pivot in the diagonal part give g_k in terms of the
            # solutions so far, rows without are linear conditions on those solutions, and the diagonal
            # part's free columns start new solutions at z^(indicials + k), where the indicial ideal vanishes.
            reduced = diagonal.augment(known).rref()
            conditions = [row[size:] for row in reduced.rows() if row[:size] == 0 and row[size:] != 0]
            if conditions:
                # Some combinations of the solutions so far can't be continued past this point - e.g. as they
                # would need a higher power of log(z) than max_log_degree. Keep those that can, by rewriting
                # every g_j in a basis of the combinations satisfying all conditions.
                restriction = matrix(QQ, conditions).right_kernel_matrix().transpose()
                for j in solved:
                    solved[j] = solved[j] * restriction
                dimension = restriction.ncols()
                reduced = diagonal.augment(known * restriction).rref()
            pivots = [column for column in reduced.pivots() if column < size]
            free = [column for column in range(size) if column not in pivots]
            # The free coefficients are 0 in the continuations of the solutions so far, so that no tail
            # contains the leading term of a solution starting here; each new solution has one of them 1.
            particular = matrix(QQ, size, dimension)
            new_solutions = matrix(QQ, size, len(free))
            for row_index, column in enumerate(pivots):
                particular.set_row(column, -reduced.row(row_index)[size:])
                for j, free_column in enumerate(free):
                    new_solutions[column, j] = -reduced[row_index, free_column]
            for j, free_column in enumerate(free):
                new_solutions[free_column, j] = 1
            if free:
                for j in solved:
                    solved[j] = solved[j].augment(matrix(QQ, size, len(free)))
            solved[k] = particular.augment(new_solutions)
            dimension += len(free)
        if dimension == 0:
            return []
        # Canonical basis: the reduced echelon form of all solutions' coefficients, ordered by z-monomial (in
        # z-order) and then log-monomial. Each solution then has its own leading coefficient 1 at which all
        # others vanish, independently of the order in which the solutions turned up above.
        coefficient_labels = [(k, log_index) for k in z_indices for log_index in log_indices]
        echelon = block_matrix(QQ, [[solved[k]] for k in z_indices], subdivide=False).transpose().rref()
        solutions = []
        for solution_vector in echelon.rows()[:dimension]:
            by_log_index = {}
            for (k, log_index), coeff in zip(coefficient_labels, solution_vector):
                by_log_index.setdefault(log_index, {})[k] = coeff
            # Every z-index up to order is kept, zeros included, for each log-monomial that occurs.
            solutions.append({
                log_index: z_dict for log_index, z_dict in by_log_index.items()
                if any(coeff != 0 for coeff in z_dict.values())
            })
        return solutions
 def _z_order_key(z_index) -> tuple:
    # The order on z-monomials used throughout: by total degree, ties broken lexicographically. Being graded,
    # it is compatible with truncating at a total degree; being a monomial order (p < q implies p + r < q + r),
    # it acts like a generic weight vector, so an operator's initial (lowest) part is a single z-monomial
    # times a theta-polynomial, and solving in this order only ever needs earlier unknowns.
    return (sum(z_index),) + tuple(z_index)
 def _left_multiply(parts: dict, z_shift, theta_power, theta_ring) -> dict:
    # Normal-ordered product (z^a * theta^c) * (sum_p z^p * P_p(theta)) = sum_p z^(a + p) * (theta + p)^c * P_p(theta)
    # of theta-parts (see PFOperator._theta_parts), by theta_i * z^p = z^p * (theta_i + p_i).
    theta_gens = theta_ring.gens()
    n = len(theta_gens)
    product = {}
    for p, P in parts.items():
        shifted_power = prod((theta_gens[i] + p[i]) ** theta_power[i] for i in range(n))
        product[tuple(z_shift[i] + p[i] for i in range(n))] = shifted_power * P
    return product
 def _indicial_roots(indicial_ideal) -> list:
    # The roots of a zero-dimensional indicial ideal with their multiplicities, as a sorted list of
    # (rho, multiplicity): each primary component belongs to one root, and its colength is that multiplicity.
    theta_gens = indicial_ideal.ring().gens()
    roots = []
    for primary, prime in indicial_ideal.complete_primary_decomposition():
        points = prime.variety()
        if len(points) != 1 or prime.vector_space_dimension() != 1:
            raise NotImplementedError(
                "The indicial ideal has a non-rational root (prime component %s), but Picard-Fuchs systems of "
                "geometric families only have rational indicials (quasi-unipotent monodromy)." % (prime,)
            )
        roots.append((tuple(points[0][theta] for theta in theta_gens), primary.vector_space_dimension()))
    return sorted(roots)
 def _leading_term(coefficients: dict, indicials) -> str:
    # The terms of a solution (a coefficients dictionary, see Period) at its z-order-smallest z-monomial,
    # e.g. 'sqrt(z1)*(log(z0) + log(z1))', to log the structure of a basis.
    n = len(indicials)
    z_gens = [SR.var('z%d' % i) for i in range(n)]
    leading = min(
        (z_index for z_dict in coefficients.values() for z_index, coeff in z_dict.items() if coeff != 0),
        key=_z_order_key,
    )
    log_polynomial = sum(
        z_dict[leading] * prod(log(z_gens[i]) ** log_index[i] for i in range(n))
        for log_index, z_dict in coefficients.items()
    )
    return str(prod(z_gens[i] ** (indicials[i] + leading[i]) for i in range(n)) * log_polynomial)
@@ -5,7 +5,7 @@
 TEST_MODELS = [
    (
        "K3_two_parameter_family",
-        lambda: ToricPolytopeProjectiveSpace([1, 1, 2, 4], model_name="K3_two_parameter_family"),
+        lambda: ToricPolytopeProjectiveSpace([4, 2, 1, 1], model_name="K3_two_parameter_family"),
        2,  # cy_dimension
        2,  # no_divs
        [("z2 - 1/4", 2), ("4096*z1^2*(4*z2 - 1) + 128*z1 - 1", 0)],
@@ -13,7 +13,7 @@ TEST_MODELS = [
    ),
    (
        "CY3_two_parameter_family",
-        lambda: ToricPolytopeProjectiveSpace([1, 1, 2, 2, 2], model_name="CY3_two_parameter_family"),
+        lambda: ToricPolytopeProjectiveSpace([2, 2, 2, 1, 1], model_name="CY3_two_parameter_family"),
        3,
        2,
        [("z2 - 1/4", 3), ("65536*z1^2*(4*z2 - 1) + 512*z1 - 1", 0)],
@@ -21,7 +21,7 @@ TEST_MODELS = [
    ),
    (
        "CY4_two_parameter_family",
-        lambda: ToricPolytopeProjectiveSpace([1, 1, 1, 1, 8, 12], model_name="CY4_two_parameter_family"),
+        lambda: ToricPolytopeProjectiveSpace([12, 8, 1, 1, 1, 1], model_name="CY4_two_parameter_family"),
        4,
        2,
        [
@@ -1,10 +1,108 @@
 import logging
 import os
 import pytest
 from sage.all import *  # noqa: F401
 load(os.path.join(os.path.dirname(__file__), "..", "sage", "periods.sage"))
 def _quintic_ideal():
    # The quintic's Picard-Fuchs operator, of rank 4 with log structure [1, 1, 1, 1] at its MUM point.
    return PFIdeal(
        [PFOperator("theta0^4 - 5*z0*(5*theta0 + 1)*(5*theta0 + 2)*(5*theta0 + 3)*(5*theta0 + 4)", no_variables=1)],
        hodge_numbers=[1, 1, 1, 1],
    )
 def _p22211_ideal():
    # The Picard-Fuchs ideal of P_{22211}[8] at its MUM point (see
    # test_ideal_finds_power_series_solution_and_recovers_first_operator), with h^{2,1} = 2.
    D1 = PFOperator(
        "theta0^2*(theta0 - 2*theta1) - 4*z0*(4*theta0 + 3)*(4*theta0 + 2)*(4*theta0 + 1)",
        no_variables=2,
    )
    D2 = PFOperator(
        "theta1^2 - z1*(2*theta1 - theta0 + 1)*(2*theta1 - theta0)",
        no_variables=2,
    )
    return PFIdeal([D1, D2], hodge_numbers=[1, 2, 2, 1])
 def _p333111_ideal():
    # The Picard-Fuchs ideal of the fourfold P_{333111}[12] at its MUM point: the GKZ operators of
    # l^(1) = (-4; 1,1,1,0,0,0,1) and l^(2) = (0; 0,0,0,1,1,1,-3), the first without its left factor theta0.
    # With the intersection numbers kappa_0000 = 36, kappa_0001 = 12, kappa_0011 = 4 (and 0 else) that
    # toric_topdata computes, the vertical h^{2,2} is 3.
    D1 = PFOperator(
        "theta0^2*(theta0 - 3*theta1) - 4*z0*(4*theta0 + 1)*(4*theta0 + 2)*(4*theta0 + 3)",
        no_variables=2,
    )
    D2 = PFOperator(
        "theta1^3 - z1*(theta0 - 3*theta1)*(theta0 - 3*theta1 - 1)*(theta0 - 3*theta1 - 2)",
        no_variables=2,
    )
    return PFIdeal([D1, D2], hodge_numbers=[1, 2, 3, 2, 1])
 def _starting_term(period):
    # A solution's lowest z-monomial (by total degree, then lexicographically) together with its coefficient
    # per log-monomial there, e.g. ((1,), {(1,): 1}) for z0*log(z0) + ...
    z_index = min(
        (z_index for z_dict in period.coefficients.values() for z_index, coeff in z_dict.items() if coeff != 0),
        key=lambda z_index: (sum(z_index), z_index),
    )
    return z_index, {
        log_index: z_dict[z_index] for log_index, z_dict in period.coefficients.items() if z_dict[z_index] != 0
    }
 def _nonzero(coefficients):
    # A coefficients dictionary without its zero entries, to compare bases that store zeros differently.
    return {
        log_index: {z_index: coeff for z_index, coeff in z_dict.items() if coeff != 0}
        for log_index, z_dict in coefficients.items()
        if any(coeff != 0 for coeff in z_dict.values())
    }
 def _truncated(f, degree):
    # A polynomial without its terms of total degree above degree.
    return sum((coeff * monomial for coeff, monomial in f if monomial.degree() <= degree), f.parent()(0))
 def _inverse_of_linear(constant, linear, degree):
    # 1/(constant + linear) for a linear form in eps, as a polynomial truncated at the given degree.
    return sum(((-linear / constant) ** j for j in range(degree + 1)), linear.parent()(0)) / constant
 def _classical_frobenius_periods(coefficient_ratio, starting_terms, order):
    # The classical Frobenius method. With coefficient_ratio(k) = c(k + eps)/c(eps) for the Gamma-function
    # coefficients c of the fundamental period (a polynomial in eps, truncated at the starting terms' degree),
    # the period starting with Q(log z) is Q(d/deps) applied to sum_k c(k + eps)/c(eps) * z^(k + eps) at eps = 0.
    n = starting_terms[0].parent().ngens()
    log_ring = PolynomialRing(QQ, n, ["L%d" % i for i in range(n)])
    logs = log_ring.gens()
    periods = [{} for _ in starting_terms]
    for k in _multi_indices(n, order):
        ratio = coefficient_ratio(k)
        for period, starting_term in zip(periods, starting_terms):
            # d^alpha/deps^alpha of ratio * exp(eps . L) at eps = 0, summed over the monomials eps^alpha of Q
            value = log_ring(0)
            for q, q_monomial in starting_term:
                alpha = q_monomial.exponents()[0]
                for r, r_monomial in ratio:
                    beta = r_monomial.exponents()[0]
                    if all(beta[i] <= alpha[i] for i in range(n)):
                        value += q * r * prod(
                            factorial(alpha[i]) // factorial(alpha[i] - beta[i]) * logs[i] ** (alpha[i] - beta[i])
                            for i in range(n)
                        )
            for coeff, log_monomial in value:
                period.setdefault(tuple(log_monomial.exponents()[0]), {})[k] = coeff
    return periods
 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
@@ -103,3 +201,247 @@ def test_single_operator_recovers_itself_from_its_holomorphic_period():
    assert ratio != 0
    assert recovered[0].operator == ratio * L.operator
 def test_indicials_of_quintic_at_mum_conifold_and_orbifold_points():
    # At the MUM point all four indicials are 0; at the conifold point they are 0, 1, 1, 2, the doubled 1
    # carrying the single logarithmic solution; at the orbifold point z0 -> 1/z0 they are 1/5, ..., 4/5.
    ideal = _quintic_ideal()
    assert ideal.indicials() == [((0,), 4)]
    assert ideal.change_coordinates("1 - 3125*z0").indicials() == [((0,), 1), ((1,), 2), ((2,), 1)]
    assert ideal.change_coordinates("1/z0").indicials() == [((QQ(k) / 5,), 1) for k in range(1, 5)]
 def test_indicials_reject_hodge_numbers_inconsistent_with_the_operators():
    # The quintic's operator has rank 4, so its indicials can't add up to the rank 5 of these Hodge numbers.
    ideal = PFIdeal(_quintic_ideal().operators, hodge_numbers=[1, 1, 1, 1, 1])
    with pytest.raises(ValueError):
        ideal.indicials()
 def test_frobenius_basis_of_quintic_at_mum_point():
    # The solutions start with 1, log(z0), log(z0)^2 and log(z0)^3, the first being the classic holomorphic
    # period. Each solution's top-log coefficient is a solution without logs, hence that period again.
    ideal = _quintic_ideal()
    basis = ideal.frobenius_basis(order=4)
    assert [_starting_term(period) for period in basis] == [((0,), {(j,): 1}) for j in range(4)]
    holomorphic = {(0,): 1, (1,): 120, (2,): 113400, (3,): 168168000, (4,): 305540235000}
    assert basis[0].coefficients == {(0,): holomorphic}
    assert all(period.coefficients[(j,)] == holomorphic for j, period in enumerate(basis))
    assert all(period.apply_operator(op).coefficients == {} for period in basis for op in ideal.operators)
 def test_frobenius_basis_of_quintic_at_conifold_point():
    # The indicials 0, 1, 1, 2 differ by integers, so all four solutions come out of one solve from exponent 0.
    # Picard-Lefschetz: exactly one of them is logarithmic, z0*log(z0) + ..., starting where the vanishing
    # period z0 + ... does.
    conifold = _quintic_ideal().change_coordinates("1 - 3125*z0")
    basis = conifold.frobenius_basis(order=4)
    assert [_starting_term(period) for period in basis] == [
        ((0,), {(0,): 1}),
        ((1,), {(0,): 1}),
        ((1,), {(1,): 1}),
        ((2,), {(0,): 1}),
    ]
    assert [max(log_index[0] for log_index in period.coefficients) for period in basis] == [0, 0, 1, 0]
    assert all(period.apply_operator(op).coefficients == {} for period in basis for op in conifold.operators)
 def test_frobenius_basis_of_quintic_at_orbifold_point():
    # Four exponent classes 1/5, ..., 4/5 mod 1, each with a single solution without logs. Order 6 reaches
    # the step at which Sage's echelon form of sparse rational matrices went wrong (see operators.sage).
    orbifold = _quintic_ideal().change_coordinates("1/z0")
    basis = orbifold.frobenius_basis(order=6)
    assert [period.indicials for period in basis] == [[QQ(k) / 5] for k in range(1, 5)]
    assert all(list(period.coefficients) == [(0,)] for period in basis)
    assert all(period.apply_operator(op).coefficients == {} for period in basis for op in orbifold.operators)
 def test_frobenius_basis_of_two_parameter_model_at_mum_point():
    # All six indicials are (0, 0), and the solutions start with 1; log(z0), log(z1); log(z0)^2,
    # log(z0)*log(z1); and a triple-log term proportional to kappa_ijk*log(z_i)*log(z_j)*log(z_k) with the
    # intersection numbers kappa_000 = 8, kappa_001 = 4: the log structure [1, 2, 2, 1] of the Hodge numbers.
    ideal = _p22211_ideal()
    assert ideal.indicials() == [((0, 0), 6)]
    basis = ideal.frobenius_basis(order=3)
    assert [_starting_term(period) for period in basis] == [
        ((0, 0), {(0, 0): 1}),
        ((0, 0), {(1, 0): 1}),
        ((0, 0), {(0, 1): 1}),
        ((0, 0), {(2, 0): 1}),
        ((0, 0), {(1, 1): 1}),
        ((0, 0), {(3, 0): 1, (2, 1): QQ(3) / 2}),
    ]
    assert all(period.apply_operator(op).coefficients == {} for period in basis for op in ideal.operators)
 def test_frobenius_basis_of_two_parameter_model_at_strong_coupling():
    # At (z0, z1) = (0, 1/4) the generators' z-free parts theta0^2*theta1 and theta1*(2*theta1 - 2*theta0 - 1)
    # vanish on the whole line theta1 = 0. Completing the indicial ideal pins the indicials down to (0, 0)
    # with multiplicity 4 and (0, 1/2) with multiplicity 2. The latter class holds sqrt(z1)*(1 + ...), the
    # power series solution of test_ideal_finds_power_series_solution_and_recovers_first_operator, and a
    # logarithmic partner sqrt(z1)*(log(z0) + log(z1) + ...).
    moved = _p22211_ideal().change_coordinates("z0", "z1 - 1/4")
    assert moved.indicials() == [((0, 0), 4), ((0, QQ(1) / 2), 2)]
    basis = moved.frobenius_basis(order=4)
    assert [period.indicials for period in basis] == [[0, 0]] * 4 + [[0, QQ(1) / 2]] * 2
    assert [_starting_term(period) for period in basis] == [
        ((0, 0), {(0, 0): 1}),
        ((0, 0), {(1, 0): 1}),
        ((0, 0), {(2, 0): 1}),
        ((0, 0), {(3, 0): 1}),
        ((0, 0), {(0, 0): 1}),
        ((0, 0), {(1, 0): 1, (0, 1): 1}),
    ]
    power_series = moved.find_power_series_solution(indicials=[0, QQ(1) / 2], order=4)
    assert basis[4].coefficients == power_series[0].coefficients
    assert all(period.apply_operator(op).coefficients == {} for period in basis for op in moved.operators)
 def test_frobenius_basis_warns_when_log_degree_is_capped_too_low(caplog):
    # The triple-log solution at the MUM point needs log-degree 3, so with max_log_degree=2 only five of the
    # six solutions survive - which is reported rather than silently returned.
    ideal = _p22211_ideal()
    with caplog.at_level(logging.WARNING):
        basis = ideal.frobenius_basis(order=2, max_log_degree=2)
    assert len(basis) == 5
    assert any("rank sum(hodge_numbers) is 6" in record.getMessage() for record in caplog.records)
 def test_frobenius_basis_of_two_parameter_fourfold(caplog):
    # At the MUM point the nine solutions have the log structure [1, 2, 3, 2, 1], and the top ones are given by
    # the intersection numbers: the quartic one by kappa_ijkl*log(z_i)*...*log(z_l), i.e. 36*log(z0)^4 +
    # 48*log(z0)^3*log(z1) + 24*log(z0)^2*log(z1)^2 up to scale, a cubic one by its derivative in log(z1).
    ideal = _p333111_ideal()
    assert ideal.indicials() == [((0, 0), 9)]
    basis = ideal.frobenius_basis(order=3)
    assert [_starting_term(period) for period in basis] == [
        ((0, 0), {(0, 0): 1}),
        ((0, 0), {(1, 0): 1}),
        ((0, 0), {(0, 1): 1}),
        ((0, 0), {(2, 0): 1}),
        ((0, 0), {(1, 1): 1}),
        ((0, 0), {(0, 2): 1}),
        ((0, 0), {(3, 0): 1, (2, 1): 1}),
        ((0, 0), {(1, 2): 1}),
        ((0, 0), {(4, 0): 1, (3, 1): QQ(4) / 3, (2, 2): QQ(2) / 3}),
    ]
    assert all(period.apply_operator(op).coefficients == {} for period in basis for op in ideal.operators)
    # At (0, -1/27), where the strong coupling divisor 1 + 27*z1 = 0 crosses the MUM divisor, all indicials lie
    # in one exponent class and there is a log(z0)^4 solution. But some indicials are (0, 1): this is no MUM
    # point, and its log structure [2, 3, 2, 1, 1] must not be compared with the Hodge numbers.
    moved = ideal.change_coordinates("z0", "z1 + 1/27")
    with caplog.at_level(logging.WARNING):
        basis = moved.frobenius_basis(order=3)
    assert moved.indicials() == [((0, 0), 5), ((0, 1), 4)]
    assert len(basis) == 9
    assert all(period.apply_operator(op).coefficients == {} for period in basis for op in moved.operators)
    assert not caplog.records
 def test_frobenius_bases_at_mum_points_agree_with_the_classical_frobenius_method():
    # An independent construction of every period at the MUM points, logarithmic ones included: the classical
    # Frobenius method, from nothing but the Gamma-function coefficients c(k) of the fundamental period. The
    # starting terms are those of the canonical bases; for P_{22211}[8] the cubic one is
    # kappa_ijk*eps_i*eps_j*eps_k/8 with the intersection numbers kappa_000 = 8, kappa_001 = 4.
    eps_ring = PolynomialRing(QQ, 1, ["e0"])
    e = eps_ring.gen(0)
    def quintic_ratio(k):
        # c(n + eps)/c(eps) for c(n) = (5n)!/(n!)^5
        ratio = eps_ring(1)
        for i in range(1, 5 * k[0] + 1):
            ratio = _truncated(ratio * (5 * e + i), 3)
        for i in range(1, k[0] + 1):
            ratio = _truncated(ratio * _inverse_of_linear(i, e, 3) ** 5, 3)
        return ratio
    classical = _classical_frobenius_periods(quintic_ratio, [e**j for j in range(4)], order=8)
    basis = _quintic_ideal().frobenius_basis(order=8)
    assert [_nonzero(period.coefficients) for period in basis] == [_nonzero(c) for c in classical]
    eps_ring = PolynomialRing(QQ, 2, ["e0", "e1"])
    e0, e1 = eps_ring.gens()
    def p22211_ratio(k):
        # c(n + eps0, m + eps1)/c(eps0, eps1) for c(n, m) = (4n)!/((n!)^3 (m!)^2 (n - 2m)!). For n < 2m the
        # last factor is Gamma(1 + delta)/Gamma(1 + delta + n - 2m) with delta = eps0 - 2*eps1, a polynomial
        # vanishing at eps = 0, which only the logarithmic periods see.
        n, m = k
        delta = e0 - 2 * e1
        ratio = eps_ring(1)
        for i in range(1, 4 * n + 1):
            ratio = _truncated(ratio * (4 * e0 + i), 3)
        for i in range(1, n + 1):
            ratio = _truncated(ratio * _inverse_of_linear(i, e0, 3) ** 3, 3)
        for i in range(1, m + 1):
            ratio = _truncated(ratio * _inverse_of_linear(i, e1, 3) ** 2, 3)
        for i in range(1, n - 2 * m + 1):
            ratio = _truncated(ratio * _inverse_of_linear(i, delta, 3), 3)
        for i in range(2 * m - n):
            ratio = _truncated(ratio * (delta - i), 3)
        return ratio
    starting_terms = [eps_ring(1), e0, e1, e0**2, e0 * e1, e0**3 + QQ(3) / 2 * e0**2 * e1]
    classical = _classical_frobenius_periods(p22211_ratio, starting_terms, order=6)
    basis = _p22211_ideal().frobenius_basis(order=6)
    assert [_nonzero(period.coefficients) for period in basis] == [_nonzero(c) for c in classical]
 def test_quintic_bases_at_mum_and_conifold_point_continue_into_each_other():
    # Both expansions converge on the lens between z0 = 0 and z0 = 1/3125, where the MUM periods must be
    # constant linear combinations T of the conifold periods. Through T the conifold monodromy
    # (log w -> log w + 2*pi*i) becomes a Picard-Lefschetz transvection in the MUM basis, and its product with
    # the MUM monodromy is the monodromy around z0 = oo, of order 5.
    CF = ComplexField(200)
    two_pi_i = 2 * CF.pi() * CF.gen()
    ideal = _quintic_ideal()
    mum = ideal.frobenius_basis(order=100)
    conifold = ideal.change_coordinates("1 - 3125*z0").frobenius_basis(order=100)
    def values(basis, points, shifted=False):
        # One row per point z0, one column per period, continued once around w = 0 if shifted.
        rows = []
        for z0 in points:
            w = z0 if basis is mum else 1 - 3125 * z0
            log_w = CF(w).log() + (two_pi_i if shifted else 0)
            rows.append([
                sum(
                    CF(coeff) * log_w ** log_index[0] * ((z_index[0] + period.indicials[0]) * log_w).exp()
                    for log_index, z_dict in period.coefficients.items()
                    for z_index, coeff in z_dict.items()
                    if coeff != 0
                )
                for period in basis
            ])
        return matrix(CF, rows)
    fit = [CF(x, y) / 3125 for x, y in [(0.4, -0.15), (0.5, 0.1), (0.6, -0.1), (0.45, 0.15)]]
    test = [CF(x, y) / 3125 for x, y in [(0.55, 0.05), (0.42, -0.05), (0.6, 0.12)]]
    T = values(conifold, fit).solve_right(values(mum, fit)).transpose()
    assert (values(mum, test) - values(conifold, test) * T.transpose()).norm() < 1e-8 * values(mum, test).norm()
    M_0 = values(mum, fit).solve_right(values(mum, fit, shifted=True)).transpose()
    M_c = T * values(conifold, fit).solve_right(values(conifold, fit, shifted=True)).transpose() * T.inverse()
    identity = identity_matrix(CF, 4)
    singular_values = (M_c - identity).change_ring(CDF).SVD()[1].diagonal()
    assert singular_values[1] < 1e-8 * singular_values[0]
    assert ((M_0 * M_c) ** 5 - identity).norm() < 1e-8