Implementing Frobenius basis for PFIdeals

Refactor and document operators and periods in sage module (#21 )
--------- Co-authored-by: Julian Piribauer <julian.piribauer@gmail.com> Reviewed-on: #21
2026-11-01 16:46:13 +01:00 · 2026-10-31 22:22:38 +01:00 · 2026-10-27 14:43:37 +01:00 · 2026-09-09 22:03:29 +02:00 · 2026-08-17 19:51:31 +02:00 · 2026-07-29 20:30:12 +02:00
13 changed files with 2141 additions and 6 deletions
@@ -0,0 +1,17 @@
+# Custom CI image for the arm64 (Raspberry Pi) Gitea Actions runner,
+# since the official sagemath/sagemath image is amd64-only.
+#
+# Build & push commands:
+#   docker build -t gitea.piribauer.ch/julian/sage-ci:latest -f .gitea/ci-image/Dockerfile .
+#   docker login gitea.piribauer.ch -u julian
+#   docker push gitea.piribauer.ch/julian/sage-ci:latest
+
+FROM condaforge/miniforge3:latest
+
+RUN mamba install -y -c conda-forge sage ruff pytest \
+    && mamba clean -afy
+
+RUN mamba install -y -c conda-forge nodejs \
+    && mamba clean -afy
+
+WORKDIR /workspace
@@ -0,0 +1,26 @@
+#!/usr/bin/env bash
+
+# This script is run in the CI container to preparse all .sage files into .sage.py files
+# and can also be run locally for local linting/testing.
+
+set -e
+
+if command -v python >/dev/null 2>&1 && python -c "import sage.repl" >/dev/null 2>&1; then
+  run_python() { python "$@"; }
+else
+  run_python() { sage --python "$@"; }
+fi
+
+for f in $(find . -name "*.sage" -not -path "./playground/*"); do
+  run_python -c "
+import sys
+from sage.repl.preparse import preparse_file
+
+path = sys.argv[1]
+with open(path) as fh:
+    src = fh.read()
+with open(path + '.py', 'w') as fh:
+    fh.write('from sage.all import *  # noqa: F401,F403\n')
+    fh.write(preparse_file(src))
+" "$f"
+done
@@ -0,0 +1,37 @@
+name: CI
+
+on:
+  push:
+    branches: ["**"]
+  pull_request:
+    branches: ["**"]
+
+jobs:
+  lint:
+    runs-on: ubuntu-latest
+    # Custom image (see .gitea/ci-image/Dockerfile)
+    container:
+      image: gitea.piribauer.ch/julian/sage-ci:latest
+
+    steps:
+      - name: Checkout
+        uses: actions/checkout@v4
+
+      - name: Preparse .sage files
+        run: .gitea/preparse.sh
+
+      - name: ruff check
+        run: python -m ruff check --no-respect-gitignore .
+
+  test:
+    runs-on: ubuntu-latest
+    needs: lint
+    container:
+      image: gitea.piribauer.ch/julian/sage-ci:latest
+
+    steps:
+      - name: Checkout
+        uses: actions/checkout@v4
+
+      - name: Run tests
+        run: python -m pytest tests/
@@ -0,0 +1,10 @@
+# Sage preparser output (regenerated from .sage sources, not hand-maintained)
+*.sage.py
+
+# Python / tooling caches
+__pycache__/
+.mypy_cache/
+.ruff_cache/
+.pytest_cache/
+
+.venv/
@@ -1,2 +1,191 @@
 # Calabi-Yau-Period-Geometry

+This project collects code used for analysing Calabi&ndash;Yau families.
+It allows for computation of discriminant loci and topological data of manifolds defined as hypersurfaces or complete intersections in toric ambient spaces.
+
+## toric_topdata
+
+Supported initialisations are represented by the following examples.
+
+```python
+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")
+CICY3_two_parameter_manual_nef = ToricPolytope(
+    [
+        [1, 0, 0, 0, 0, 0],
+        [0, 1, 0, 0, 0, 0],
+        [0, 0, 1, 0, 0, 0],
+        [0, 0, 0, 1, 0, 0],
+        [0, 0, 0, 0, 1, 0],
+        [0, 0, 0, 0, 0, 1],
+        [-1, -1, 0, 0, 0, 0],
+        [0, 0, -1, -1, -1, -1],
+    ],
+    nef_partition=[[0, 1, 2, 3], [4, 5, 6, 7]],
+)
+
+CICY5_two_parameter = ToricPolytopeCICY([[6, 1], [0, 2]])
+```
+
+The discriminant factors and topological data, e.g. for the quintic, can then be computed with the methods below.
+```python
+CY3_quintic.disc()
+
+CY3_quintic.topdata()
+```
+
+We list the output of above two lines.
+
+```term
+INFO:__main__:Discriminant factors: [[z1 + 1/3125, 0]]
+INFO:__main__:
+
+--- Polytope and GLSM table ------------------
+ 1|1|1|1|1|| 1
+
+-1|1|0|0|0|| 0
+-1|0|1|0|0|| 0
+-1|0|0|1|0|| 0
+-1|0|0|0|1|| 0
+--------------
+ 1,1,1,1,1; -5
+
+--- L-vectors (GLSM charges) -----------------
+[[1, 1, 1, 1, 1, -5]]
+
+--- Intersection numbers CY ------------------
+{(0, 0, 0): 5}
+
+--- Intersection ring ------------------------
+5*t0^3
+
+--- Intersection ring no multiplicities ------
+5*t0^3
+
+--- Chern polynomials ------------------------
+[[0], [10*t0^3], [-40*t0^3]]
+
+--- Integrated Chern classes -----------------
+[[0], [50], [-200]]
+DEBUG:__main__:
+
+--- Kähler cone generators (ambient space) ---
+['[z4]']
+
+--- Intersection numbers ambient space -------
+{(0, 0, 0, 0): 1}
+```
+
+The result is saved in the folder `data/topdata` as a JSON file &mdash; giving a model name helps keeping
+track of these outputs.
+Note that for Calabi&ndash;Yau dimensions larger than four, the additional
+
+## Operators and periods
+
+`operators.sage` and `periods.sage` provide classes for working with Picard&ndash;Fuchs operators and their
+period solutions: `PFOperator`, `PFIdeal` and `Period`, together with Ansatz variants of the first and
+last (`PFOperatorAnsatz`, `PeriodAnsatz`) used to search for unknown operators or periods of a given
+z- and theta-degree. The operator classes live in `operators.sage`, the period classes in `periods.sage`;
+since the two depend on each other, load `periods.sage`, which also loads `operators.sage`.
+
+A `PFOperator` is parsed from a string in the variables `z0, ..., z<n-1>` and `theta0, ..., theta<n-1>`,
+the logarithmic derivatives theta_i = z_i d/dz_i. For example, the quintic's Picard&ndash;Fuchs operator:
+
+```python
+L = PFOperator(
+    "theta0^4 - 3125*z0*theta0^4 - 6250*z0*theta0^3 - 4375*z0*theta0^2 - 1250*z0*theta0 - 120*z0",
+    no_variables=1,
+)
+
+L.simplify().operator_string
+```
+
+```term
+'-5*(5*theta0 + 4)*(5*theta0 + 3)*(5*theta0 + 2)*(5*theta0 + 1)*z0 + theta0^4'
+```
+
+An operator (or a `PFIdeal` of several) can be solved for its power series solution at given indicial
+exponents and order.
+
+```python
+ideal = PFIdeal([L])
+period = ideal.find_power_series_solution(indicials=[0], order=3)[0]
+
+period.period_string
+```
+
+```term
+'168168000*z0^3 + 113400*z0^2 + 120*z0 + 1'
+```
+
+The reverse direction is supported too: given a `Period`, `find_annihilating_operators` searches for
+`PFOperator`s of a given z- and theta-degree that annihilate it, by solving an Ansatz of unknown
+coefficients via linear algebra. Both directions extend to several moduli, e.g. for the two-parameter
+model P_{2,2,2,1,1}[8]:
+
+```python
+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])
+
+period = ideal.find_power_series_solution(indicials=[0, 1 / 2], order=6)[0]
+recovered = period.find_annihilating_operators(z_degree=1, theta_degree=2)
+
+period.period_string
+recovered[0].simplify().operator_string
+```
+
+```term
+'-1/45045*(60886425600*z0^3*z1^3 - 2767564800*z0^2*z1^4 + 100638720*z0*z1^5 - 14192640*z1^6 + 830269440*z0^2*z1^3 - 30750720*z0*z1^4 + 4193280*z1^5 - 242161920*z0^2*z1^2 + 9884160*z0*z1^3 - 1281280*z1^4 - 3459456*z0*z1^2 + 411840*z1^3 + 1441440*z0*z1 - 144144*z1^2 + 60060*z1 - 45045)*sqrt(z1)'
+'-2*(theta0 - 2*theta1)*(theta0 - 2*theta1 - 1)*z1 + (2*theta0 - 2*theta1 + 1)*theta1'
+```
+
+`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.
@@ -0,0 +1,20 @@
+[tool.ruff]
+line-length = 120
+target-version = "py311"
+
+[tool.ruff.lint]
+# E: pycodestyle errors, F: pyflakes (undefined/unused names), I: isort (import order), W: pycodestyle warnings
+select = ["E", "F", "I", "W"]
+ignore = [
+    "E501", # symbolic-math expressions routinely exceed a "normal" line length
+]
+
+[tool.ruff.lint.per-file-ignores]
+# *.sage.py is preparser output: star-import globals, semicolon preamble, and literal-substitution artifacts trip static analysis.
+"*.sage.py" = ["F403", "F405", "F821", "E741", "E402", "E702", "I001", "W291", "W293", "W292"]
+
+# test_smoke.py star-imports sage.all and loads a .sage file, so ruff can't see where its names come from.
+"tests/test_topdata_and_disc.py" = ["F403", "F405"]
+
+# Same star-import + load() pattern as above.
+"tests/test_period_computation.py" = ["F403", "F405"]
@@ -0,0 +1,811 @@
+"""
+Picard-Fuchs operators and ideals in the logarithmic derivatives theta_i = z_i * d/dz_i.
+
+Contents
+--------
+PFOperator
+    __init__                      Construct from a string in z- and theta-variables
+    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 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, ZZ, binomial, block_matrix, floor, function, log, matrix, prod, solve, var
+)
+
+load("sage/util.py")
+
+# Logger
+logger = logging.getLogger(__name__)
+logging.basicConfig(level=logging.DEBUG)
+
+class PFOperator:
+    """
+    A class representing a Picard-Fuchs operator in a given number of variables (moduli).
+    Independent of the given order, the variables are assumed to be left of
+    the derivatives.
+
+    The variables extra_locals are needed for symbolic coefficients used for
+    operator Ansätze.
+    """
+    def _operator_from_string(self, operator_string: str, extra_locals: dict = None):
+        z_names = ['z%d' % i for i in range(self.no_variables)]
+        theta_names = ['theta%d' % i for i in range(self.no_variables)]
+        base_ring = SR if extra_locals else QQ
+        self.ring = PolynomialRing(base_ring, z_names + theta_names)
+        self.z_gens = self.ring.gens()[:self.no_variables]
+        self.theta_gens = self.ring.gens()[self.no_variables:]
+
+        parse_locals = dict(self.ring.gens_dict())
+        if extra_locals:
+            parse_locals.update(extra_locals)
+
+        try:
+            operator = self.ring(sage_eval(operator_string, locals=parse_locals))
+            return operator
+        except Exception as e:
+            raise ValueError("Invalid operator string: %s" % e)
+
+    def __init__(self, operator_string: str, no_variables: int = 1, extra_locals: dict = None):
+        self.no_variables = no_variables
+        self.operator_string = operator_string
+        self.operator = self._operator_from_string(operator_string, extra_locals=extra_locals)
+        logger.info("Initialised PFOperator: %s", self.operator)
+
+    def simplify(self) -> "PFOperator":
+        """
+        Groups the operator's terms by z-monomial and factorises the theta-polynomial
+        multiplying each z-monomial, e.g. turning a computed quintic operator into the
+        well-known theta0^4 - 5*z0*(5*theta0 + 1)*(5*theta0 + 2)*(5*theta0 + 3)*(5*theta0 + 4).
+        """
+        z_monomial_theta_parts = {}
+        for coeff, monomial in self.operator:
+            exponents = monomial.exponents()[0]
+            z_exponents = exponents[:self.no_variables]
+            theta_exponents = exponents[self.no_variables:]
+
+            theta_monomial = SR(1)
+            for gen, exp in zip(self.theta_gens, theta_exponents):
+                theta_monomial *= SR(gen) ** exp
+
+            z_monomial_theta_parts[z_exponents] = z_monomial_theta_parts.get(z_exponents, SR(0)) + SR(coeff) * theta_monomial
+
+        simplified_expr = SR(0)
+        for z_exponents, theta_part in z_monomial_theta_parts.items():
+            z_monomial = SR(1)
+            for gen, exp in zip(self.z_gens, z_exponents):
+                z_monomial *= SR(gen) ** exp
+            simplified_expr += theta_part.factor() * z_monomial
+
+        simplified = copy.copy(self)
+        simplified.operator_string = str(simplified_expr)
+        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
+
+    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.
+    """
+    def __init__(self, no_variables: int, theta_degree: int, z_degree: int):
+        self.no_variables = no_variables
+        self.theta_degree = theta_degree
+        self.z_degree = z_degree
+
+        z_indices = _multi_indices(self.no_variables, z_degree)
+        theta_indices = _multi_indices(self.no_variables, theta_degree)
+
+        # Every (z_index, theta_index) pair allowed by the degree bounds gets its own
+        # fresh unknown, to be solved for once the Ansatz is applied to a period.
+        operator_terms = [
+            (var("b_" + "_".join(str(x) for x in z_index + theta_index)), (z_index, theta_index))
+            for z_index in z_indices
+            for theta_index in theta_indices
+        ]
+        self.unknowns = [coeff for coeff, _ in operator_terms]
+
+        def _monomial_factors(index, name):
+            return [f"{name}{i}^{exp}" for i, exp in enumerate(index) if exp > 0]
+
+        operator_string = " + ".join(
+            "*".join([str(coeff), *_monomial_factors(z_index, "z"), *_monomial_factors(theta_index, "theta")])
+            for coeff, (z_index, theta_index) in operator_terms
+        )
+
+        extra_locals = {str(coeff): coeff for coeff in self.unknowns}
+        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.
+
+    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, 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):
+        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 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, 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 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.
+
+        This is _solve_exponent_class with no room for logarithms, on the operators as given (without
+        completing their indicial ideal). A combination of solutions that could only be continued with a
+        log(z)-term, at a resonance, is not a power series solution and is dropped.
+        """
+        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()
+
+        solutions = [
+            Period(no_variables=self.no_variables, coefficients=coefficients, order=order, indicials=indicials)
+            for coefficients in self._solve_exponent_class(operator_parts, theta_ring, indicials, order, 0)
+        ]
+
+        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)
@@ -0,0 +1,349 @@
+"""
+Periods of Calabi-Yau manifolds as formal power series in z-variables and their logs.
+
+Contents
+--------
+Period
+    __init__                      Construct from a string or a coefficients dictionary
+    from_l_vectors (classmethod)  Holomorphic period from GLSM/Mori charge vectors
+    apply_operator                Apply a PFOperator, truncated to self.order
+    find_annihilating_operators   Find PFOperator Ansätze annihilating this period
+    _initialise_ring              Polynomial ring in z- and log-variables
+    _period_from_string           Parse a string into a coefficients dictionary
+    _period_to_string             Render the coefficients as a string
+    _max_z_degree                 Highest total z-degree among the coefficients
+    _truncate_coefficients        Drop terms above a given total z-degree
+    _apply_theta                  Apply theta_i = z_i * d/dz_i to a coefficients dictionary
+
+PeriodAnsatz(Period)
+    __init__                      Ansatz of given z- and log-multi-degrees
+"""
+import logging
+from sage.all import sage_eval, PolynomialRing, QQ, SR, ceil, factorial, floor, log, matrix, prod, var
+
+load("sage/operators.sage")
+
+# Logger
+logger = logging.getLogger(__name__)
+logging.basicConfig(level=logging.DEBUG)
+
+
+class Period:
+    """
+    A class representing a period as a formal power series in z-variables and their logs.
+
+    The coefficients are stored in a dictionary of dictionaries, where the first key is the multi-index
+    of the logarithmic part and the second key is the multi-index of the z-variables. 
+    So, for example, the coefficient of
+
+        log(z0)^2 * log(z1) * z0^3 * z1^2 would be stored as coefficients[(2, 1)][(3, 2)]
+
+    """
+
+    def _initialise_ring(self):
+        z_names = ['z%d' % i for i in range(self.no_variables)]
+        log_names = ['L%d' % i for i in range(self.no_variables)]
+        ring = PolynomialRing(QQ, z_names + log_names)
+        z_gens = ring.gens()[:self.no_variables]
+        log_gens = ring.gens()[self.no_variables:]
+
+        return ring, z_gens, log_gens  # type: (PolynomialRing, list, list)
+
+    def _period_from_string(self, period_string: str) -> dict:
+        ring, z_gens, log_gens = self._initialise_ring()
+        log_of_z = dict(zip(z_gens, log_gens))
+
+        def log(zi):
+            try:
+                return log_of_z[zi]
+            except (KeyError, TypeError):
+                raise ValueError(
+                    "log(...) may only be applied to one of the z-variables z0, ..., z%d"
+                    % (self.no_variables - 1)
+                )
+
+        parse_locals = dict(ring.gens_dict())
+        parse_locals['log'] = log
+
+        try:
+            expression = ring(sage_eval(period_string, locals=parse_locals))
+        except Exception as e:
+            raise ValueError("Invalid period string: %s" % e)
+
+        coefficients = {}
+        for coeff, monomial in expression:
+            exponents = monomial.exponents()[0]
+            z_index = tuple(exponents[:self.no_variables])
+            log_index = tuple(exponents[self.no_variables:])
+            coefficients.setdefault(log_index, {})[z_index] = coeff
+
+        return coefficients
+    
+    def _period_to_string(self) -> str:
+        ring, z_gens, log_gens = self._initialise_ring()
+        expression = ring(0)
+        for log_index, z_dict in self.coefficients.items():
+            for z_index, coeff in z_dict.items():
+                monomial = coeff
+                for i in range(self.no_variables):
+                    monomial *= (log_gens[i] ** log_index[i]) * (z_gens[i] ** z_index[i])
+                expression += monomial
+
+        log_substitutions = {log_gens[i]: log(SR(z_gens[i])) for i in range(self.no_variables)}
+        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.
+        if not self.coefficients:
+            return 0
+        return max(
+            sum(z_index)
+            for z_dict in self.coefficients.values()
+            for z_index in z_dict
+        )
+
+    def _truncate_coefficients(self, order: int) -> dict:
+        # Drop every (log_index, z_index) entry whose total z-degree exceeds order,
+        # removing it from the dictionary rather than merely zeroing it out.
+        truncated = {}
+        for log_index, z_dict in self.coefficients.items():
+            kept = {z_index: coeff for z_index, coeff in z_dict.items() if sum(z_index) <= order}
+            if kept:
+                truncated[log_index] = kept
+        return truncated
+
+    def _apply_theta(self, coefficients: dict, index: int) -> dict:
+        """
+        Applies the logarithmic derivative theta_i = z_i * d/dz_i once to a coefficients dict of the same
+        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] + self.indicials[index]
+                if a_i != 0:
+                    inner = result.setdefault(log_index, {})
+                    inner[z_index] = inner.get(z_index, 0) + a_i * coeff
+
+                k_i = log_index[index]
+                if k_i != 0:
+                    lowered_log_index = log_index[:index] + (k_i - 1,) + log_index[index + 1:]
+                    inner = result.setdefault(lowered_log_index, {})
+                    inner[z_index] = inner.get(z_index, 0) + k_i * coeff
+        return result
+
+    def apply_operator(self, pf_operator: PFOperator) -> "Period":
+        """
+        Applies a PFOperator to this period and returns the result as a new Period,
+        truncated to self.order. Each monomial of the operator is normalised as
+        coeff * z^p * theta^q (see PFOperator's docstring), so theta^q is applied
+        to the period first and the result is then multiplied by coeff * z^p.
+        """
+        if pf_operator.no_variables != self.no_variables:
+            raise ValueError(
+                "Variable count mismatch: period has %d variable(s), operator has %d."
+                % (self.no_variables, pf_operator.no_variables)
+            )
+
+        result_coefficients = {}
+        for monomial_coeff, monomial in pf_operator.operator:
+            exponents = monomial.exponents()[0]
+            z_exponents = exponents[:self.no_variables]
+            theta_exponents = exponents[self.no_variables:]
+
+            term = self.coefficients
+            for i, power in enumerate(theta_exponents):
+                for _ in range(power):
+                    term = self._apply_theta(term, i)
+
+            for log_index, z_dict in term.items():
+                inner = result_coefficients.setdefault(log_index, {})
+                for z_index, coeff in z_dict.items():
+                    shifted_z_index = tuple(z_index[i] + z_exponents[i] for i in range(self.no_variables))
+                    inner[shifted_z_index] = inner.get(shifted_z_index, 0) + monomial_coeff * coeff
+
+        result_coefficients = {
+            log_index: {z_index: c for z_index, c in z_dict.items() if c != 0}
+            for log_index, z_dict in result_coefficients.items()
+        }
+        result_coefficients = {log_index: z_dict for log_index, z_dict in result_coefficients.items() if z_dict}
+
+        return 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:
+        """
+        Finds operators annihilating this period among PFOperator Ansätze of the given
+        z- and theta-degree.
+        """
+        ansatz = PFOperatorAnsatz(self.no_variables, theta_degree, z_degree)
+        unknowns = ansatz.unknowns
+
+        equations = [
+            coeff
+            for z_dict in self.apply_operator(ansatz).coefficients.values()
+            for coeff in z_dict.values()
+        ]
+
+        coefficient_rows = [[SR(equation).coefficient(b) for b in unknowns] for equation in equations]
+        kernel_basis = matrix(QQ, coefficient_rows, ncols=len(unknowns)).right_kernel().basis()
+
+        operators = []
+        for basis_vector in kernel_basis:
+            solution = {unknowns[j]: basis_vector[j] for j in range(len(unknowns))}
+            solved_operator = ansatz.operator.map_coefficients(lambda c: SR(c).subs(solution))
+            operators.append(PFOperator(str(solved_operator), no_variables=self.no_variables))
+
+        return operators
+
+    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:
+            if coefficients is not None:
+                raise ValueError("Provide either coefficients or period_string, not both.")
+            self.coefficients = self._period_from_string(period_string)
+            self.period_string = period_string
+            logger.debug("Using string for initialisation.")
+        elif coefficients is None:
+            self.coefficients = {}
+        else:
+            self.coefficients = coefficients
+            self.period_string = self._period_to_string()
+            logger.debug("Using coefficients for initialisation.")
+
+        if order is None:
+            self.order = self._max_z_degree()
+            logger.debug("Order not provided, using maximum z-degree: %d", self.order)
+        else:
+            self.coefficients = self._truncate_coefficients(order)
+            self.order = order
+            self.period_string = self._period_to_string()
+            logger.debug("Truncated coefficients to order %d.", self.order)
+
+        logger.info("Initialised Period in %d variables at order %d.", self.no_variables, self.order)
+
+    @classmethod
+    def from_l_vectors(cls, l_vectors: list, order: int) -> "Period":
+        """
+        Computes the holomorphic (fundamental) period directly from the GLSM/Mori charge
+        vectors l^(1), ..., l^(m) of a toric Calabi-Yau (m = number of moduli), via the
+        hypergeometric-type sum
+
+            varpi_0(z) = sum_{a in Z_{>=0}^m}
+                (prod_j Gamma(1 - l_num^(j).a)) / (prod_i Gamma(1 + l_den^(i).a)) * z^a.
+       """
+        if any(len(l_vector) != 2 for l_vector in l_vectors):
+            raise ValueError("Each charge vector must be given as [numerator part, denominator part].")
+        m = len(l_vectors)
+        numerator_charges = [l_vector[0] for l_vector in l_vectors]
+        denominator_charges = [l_vector[1] for l_vector in l_vectors]
+        len0 = len(numerator_charges[0]) if m else 0
+        len1 = len(denominator_charges[0]) if m else 0
+
+        # Each denominator Gamma-factor is linear in a; pivot it on the smallest index j with
+        # a nonzero coefficient there. Since a is built up from a_{m-1} down to a_0, every
+        # other index still appearing in that linear form is then already fixed, turning the
+        # factor's pole at non-positive integers into a concrete bound on a_j.
+        pivots = [[] for _ in range(m)]
+        for i in range(len1):
+            charges = [denominator_charges[k][i] for k in range(m)]
+            j = next(k for k in range(m) if charges[k] != 0)
+            pivots[j].append(charges)
+
+        coefficients = {}
+        a = [0] * m
+
+        def recurse(j, num_partial, den_partial, degree):
+            if j < 0:
+                numerator = prod(factorial(-x) for x in num_partial)
+                denominator = prod(factorial(x) for x in den_partial)
+                coefficients[tuple(a)] = QQ(numerator) / QQ(denominator)
+                return
+
+            lower, upper = 0, order - degree
+            for charges in pivots[j]:
+                c = charges[j]
+                bound = -QQ(sum(charges[k] * a[k] for k in range(j + 1, m))) / QQ(c)
+                if c < 0:
+                    upper = min(upper, floor(bound))
+                else:
+                    lower = max(lower, ceil(bound))
+
+            for a_j in range(lower, upper + 1):
+                a[j] = a_j
+                recurse(
+                    j - 1,
+                    [num_partial[jj] + numerator_charges[j][jj] * a_j for jj in range(len0)],
+                    [den_partial[ii] + denominator_charges[j][ii] * a_j for ii in range(len1)],
+                    degree + a_j,
+                )
+
+        recurse(m - 1, [0] * len0, [0] * len1, 0)
+
+        logger.info("Computed holomorphic period from l-vectors for %d moduli at order %d.", m, order)
+        return cls(no_variables=m, coefficients={tuple([0] * m): coefficients}, order=order)
+
+
+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, indicials: list = None):
+        self.no_variables = no_variables
+        self.z_degree = z_degree
+        self.log_degree = log_degree
+
+        log_indices = _multi_indices(self.no_variables, log_degree)
+        z_indices = _multi_indices(self.no_variables, z_degree)
+
+        # Every (log_index, z_index) pair allowed by the degree bounds gets its own
+        # fresh unknown, to be solved for once a PFOperator is applied to the Ansatz.
+        coefficients = {
+            log_index: {
+                z_index: var("a_" + "_".join(str(x) for x in log_index + z_index))
+                for z_index in z_indices
+            }
+            for log_index in log_indices
+        }
+
+        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).",
+            sum(len(z_dict) for z_dict in self.expansion_coefficients.values()),
+        )
@@ -1,8 +1,52 @@
+"""
+Topological data of Calabi-Yau hypersurfaces and CICYs in toric ambient spaces.
+
+Contents
+--------
+Polytope
+    __init__                             Construct from a list of defining points
+    _validate_points_and_get_dimension   Check the points and return their dimension
+    _points_to_polytope                  Lattice points of the convex hull, minus codim-1 interiors
+
+ToricPolytope(Polytope)
+    __init__                             Construct with triangulation, nef partition and Mori cone
+    set_triangulation                    Select a fine star triangulation by index
+    set_nef_partition                    Set and validate the nef partition
+    set_Mori_cone                        Set the Mori cone (default or from given l-vectors)
+    disc                                 A-discriminant (singular loci in Batyrev coordinates)
+    topdata                              Compute topological data and write it to JSON
+    _set_fine_star_triangulations        Fine star triangulations of the polytope
+    _collect_linear_dependencies         Linear relations among points in all faces
+    _append_origin_weight                Add the origin weight to each relation
+    _setup_discriminant_symbols          Symbolic variables for the discriminant
+    _build_equation_system               Equation system for one block of relations
+    _eliminate_and_normalize             Eliminate lambda-variables and normalise
+    _append_unique_discriminants         Collect discriminant factors without duplicates
+    _compute_Kahler_generators           Kähler cone generators dual to the Mori cone
+    _compute_intersection_ring           Intersection ring polynomial of the CY
+    _intersection_substitution_rules     Kähler monomial -> intersection number rules
+    _replace_intersection_monomials      Substitute intersection numbers into an expression
+    _compute_Chern_character             Total Chern class via adjunction
+    _compute_Chern_polynomials           Chern polynomials with intersection numbers
+    _check_for_fibrations                Detect possible elliptic fibrations
+    _get_polytope_and_glsm_table         Table of polytope points and GLSM charges
+    _get_lvec_rows                       GLSM charge l-vectors as lists of ints
+    _get_output                          Dictionary of all topological data
+    _to_json_safe (staticmethod)         Convert values into JSON-safe primitives
+    _write_output_json                   Write the output dictionary to JSON
+
+ToricPolytopeCICY(ToricPolytope)
+    __init__                             Construct from a CICY configuration matrix
+    _CICY_to_points                      Configuration matrix -> points and nef partition
+
+ToricPolytopeProjectiveSpace(ToricPolytope)
+    __init__                             Construct from weighted projective space weights
+    _compute_points_from_weights         Weights -> polytope points
+"""
 import numpy as np
 import logging
 import os
 import json
-from datetime import datetime, timezone
 import re

 # Logger
@@ -102,6 +146,7 @@ class ToricPolytope(Polytope):
        if no_triangulation < 0 or no_triangulation >= len(self.triangulations):
            raise IndexError("Invalid triangulation index {} (available: 0..{}).".format(no_triangulation, len(self.triangulations) - 1))
        
+        if len(self.triangulations) > 1:
            logger.info("Using triangulation index %d of [0..%d].", no_triangulation, len(self.triangulations) - 1)
        self.triangulation = self.triangulations[no_triangulation]
    
@@ -233,7 +278,7 @@ class ToricPolytope(Polytope):
    def _setup_discriminant_symbols(self):
        mori_rays = self.Mori_cone.rays()

-        a_vars = [var("a_{}".format(u), latex_name="a_{{}}".format(u)) for u in (1..len(mori_rays))]
+        a_vars = [var("a_{}".format(u), latex_name="a_{{{}}}".format(u)) for u in (1..len(mori_rays))]
        z_vars = var('z', n=len(mori_rays)+1, latex_name='z') # z[0] is superfluous
        lambda_vars = var('l', n=len(mori_rays)+1, latex_name='l') # l[0] is superfluous
        a_row = matrix(a_vars)
@@ -275,7 +320,7 @@ class ToricPolytope(Polytope):
            if polynomials[0] == 0:
                polynomials = maxima.eliminate(equation_system, lambda_symbols[:-1]).sage()
                reverse_solve = True
-        except:
+        except Exception:
            # In one-parameter cases there may be nothing to eliminate.
            polynomials = equation_system
            logger.debug("No elimination needed for the equation system: %s", equation_system)
@@ -289,7 +334,7 @@ class ToricPolytope(Polytope):
            else:
                last_lambda = lambda_symbols[-1]
                polynomials = [poly / last_lambda ** (poly.degree(last_lambda)) for poly in polynomials]
-        except:
+        except Exception:
            pass

        return polynomials
@@ -740,3 +785,34 @@ class ToricPolytopeCICY(ToricPolytope):
            lvec=lvec
        )
    
+class ToricPolytopeProjectiveSpace(ToricPolytope):
+    """
+    Computes topological data for hypersurfaces in toric ambient spaces.
+    """
+    def _compute_points_from_weights(self, weights):
+        # Convert the weights of the projective space to the corresponding points
+        dimension = len(weights) - 1
+        points = identity_matrix(dimension)
+        try:
+            weights.remove(1)
+        except ValueError:
+            raise ValueError("The weights must include a 1 for the projective space.")
+
+        points = points.insert_row(0, [-w for w in weights])
+        return list(points)
+    
+    def __init__(
+            self,
+            weights,
+            no_triangulation = 0,
+            model_name = None,
+            nef_partition = None,
+            lvec = None, # Note that the inner point is the last entry
+        ):
+        super().__init__(
+            self._compute_points_from_weights(weights),
+            no_triangulation=no_triangulation,
+            model_name=model_name,
+            nef_partition=nef_partition,
+            lvec=lvec
+        )
@@ -0,0 +1,15 @@
+"""
+Utility functions for the Calabi-Yau period geometry project.
+"""
+
+def _multi_indices(no_variables, total_degree: int) -> list:
+    # All no_variables-tuples of non-negative integers whose entries sum to at most total_degree.
+    def helper(remaining_vars, remaining_degree):
+        if remaining_vars == 0:
+            yield ()
+            return
+        for first in range(remaining_degree + 1):
+            for rest in helper(remaining_vars - 1, remaining_degree - first):
+                yield (first,) + rest
+
+    return list(helper(no_variables, total_degree))
@@ -0,0 +1,104 @@
+# Entries are (id, factory, expected_cy_dimension, expected_no_divs, expected_disc,
+# expected_intersection_numbers_CY), where expected_disc is a list of
+# (str(discriminant_factor), codimension) pairs as returned by ToricPolytope.disc().
+
+TEST_MODELS = [
+    (
+        "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)],
+        {(0, 0): 4, (0, 1): 2, (1, 1): 0},
+    ),
+    (
+        "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)],
+        {(0, 0, 0): 8, (0, 0, 1): 4, (0, 1, 1): 0, (1, 1, 1): 0},
+    ),
+    (
+        "CY4_two_parameter_family",
+        lambda: ToricPolytopeProjectiveSpace([12, 8, 1, 1, 1, 1], model_name="CY4_two_parameter_family"),
+        4,
+        2,
+        [
+            ("z2 - 1/256", 2),
+            ("34828517376*z1^4*(256*z2 - 1) + 322486272*z1^3 - 1119744*z1^2 + 1728*z1 - 1", 0),
+        ],
+        {(0, 0, 0, 0): 64, (0, 0, 0, 1): 16, (0, 0, 1, 1): 4, (0, 1, 1, 1): 1, (1, 1, 1, 1): 0},
+    ),
+    (
+        "CICY3_one_parameter",
+        lambda: ToricPolytopeCICY([[3, 3]], model_name="CICY3_one_parameter"),
+        3,
+        1,
+        [("z1 - 1/46656", 0)],
+        {(0, 0, 0): 9},
+    ),
+    (
+        "CICY3_two_parameter",
+        lambda: ToricPolytopeCICY([[3], [3]], model_name="CICY3_two_parameter"),
+        3,
+        2,
+        [
+            (
+                "-19683*z1^3 - 2187*(27*z1 + 1)*z2^2 - 19683*z2^3 - 2187*z1^2 "
+                "- 81*(729*z1^2 - 189*z1 + 1)*z2 - 81*z1 - 1",
+                0,
+            )
+        ],
+        {(0, 0, 0): 0, (0, 0, 1): 3, (0, 1, 1): 3, (1, 1, 1): 0},
+    ),
+    (
+        "CICY5_two_parameter",
+        lambda: ToricPolytopeCICY([[6, 1], [0, 2]], model_name="CICY5_two_parameter"),
+        5,
+        2,
+        [
+            (
+                "16384*z2^7 - 28672*z2^6 + 21504*z2^5 - 448*(1647086*z1 - 5)*z2^3 - 8960*z2^4 "
+                "- 112*(8235430*z1 + 3)*z2^2 - 678223072849*z1^2 - 28*(4941258*z1 - 1)*z2 - 1647086*z1 - 1",
+                0,
+            )
+        ],
+        {
+            (0, 0, 0, 0, 0): 6,
+            (0, 0, 0, 0, 1): 12,
+            (0, 0, 0, 1, 1): 0,
+            (0, 0, 1, 1, 1): 0,
+            (0, 1, 1, 1, 1): 0,
+            (1, 1, 1, 1, 1): 0,
+        },
+    ),
+    (
+        "CICY3_two_parameter_manual",
+        lambda: ToricPolytope(
+            [
+                [1, 0, 0, 0, 0, 0],
+                [0, 1, 0, 0, 0, 0],
+                [0, 0, 1, 0, 0, 0],
+                [0, 0, 0, 1, 0, 0],
+                [0, 0, 0, 0, 1, 0],
+                [0, 0, 0, 0, 0, 1],
+                [-1, -1, 0, 0, 0, 0],
+                [0, 0, -1, -1, -1, -1],
+            ],
+            nef_partition=[[0, 1, 2, 3], [4, 5, 6, 7]],
+            model_name="CICY3_two_parameter_manual",
+        ),
+        4,
+        2,
+        [
+            (
+                "-14348907*z1^5 - 2657205*z1^4 - 196830*z1^3 - 29296875*(270*z1 + 1)*z2^2 "
+                "- 30517578125*z2^3 - 7290*z1^2 + 9375*(98415*z1^3 - 32805*z1^2 + 810*z1 - 1)*z2 "
+                "- 135*z1 - 1",
+                0,
+            )
+        ],
+        {(0, 0, 0, 0): 0, (0, 0, 0, 1): 0, (0, 0, 1, 1): 6, (0, 1, 1, 1): 8, (1, 1, 1, 1): 2},
+    ),
+]
@@ -0,0 +1,447 @@
+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
+    # 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_power_series_solution_and_recovers_first_operator():
+    # 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,
+    )
+    D2 = PFOperator(
+        "theta1^2 - z1*(2*theta1 - theta0 + 1)*(2*theta1 - theta0)",
+        no_variables=2,
+    )
+    moved = PFIdeal([D1, D2]).change_coordinates("z0", "z1 - 1/4")
+    op1, op2 = moved.operators
+
+    solutions = moved.find_power_series_solution(indicials=[0, QQ(1) / 2], order=14)
+    assert len(solutions) == 1
+
+    period = solutions[0]
+    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, 0)] == 0
+
+    recovered = period.find_annihilating_operators(z_degree=1, theta_degree=2)
+    assert len(recovered) == 1
+    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():
+    # Ensuring methods work for operator 4.2.1:
+    L = PFOperator(
+        "theta0^4 - 4*z0*(2*theta0 + 1)^2*(7*theta0^2 + 7*theta0 + 2) "
+        "- 128*z0^2*(2*theta0 + 1)^2*(2*theta0 + 3)^2",
+        no_variables=1,
+    )
+    ideal = PFIdeal([L])
+
+    solutions = ideal.find_power_series_solution(indicials=[0], order=15)
+    assert len(solutions) == 1
+
+    period = solutions[0]
+    assert period.apply_operator(L).coefficients == {}
+    assert period.coefficients[(0,)][(0,)] == 1
+    assert period.coefficients[(0,)][(1,)] == 8
+    assert period.coefficients[(0,)][(2,)] == 360
+
+    # L lives entirely within z_degree <= 2, theta_degree <= 4, its stated level.
+    recovered = period.find_annihilating_operators(z_degree=2, theta_degree=4)
+    assert len(recovered) == 1
+
+    # recovered[0] should be a scalar multiple of L - compare via the coefficient of the
+    # bare theta0^4 monomial, which is 1 in L.
+    ratio = recovered[0].operator.monomial_coefficient(L.theta_gens[0] ** 4) / L.operator.monomial_coefficient(
+        L.theta_gens[0] ** 4
+    )
+    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
@@ -0,0 +1,34 @@
+import os
+
+import pytest
+from sage.all import *  # noqa: F401
+from sage.geometry.lattice_polytope import set_palp_dimension
+
+load(os.path.join(os.path.dirname(__file__), "..", "sage", "toric_topdata.sage"))
+load(os.path.join(os.path.dirname(__file__), "test_models.sage"))
+
+set_palp_dimension(11)
+
+TEST_MODELS_PARAMETRISED = [pytest.param(*row[1:], id=row[0]) for row in TEST_MODELS]
+
+
+@pytest.mark.parametrize(
+    "make_model, expected_cy_dimension, expected_no_divs, expected_disc, expected_intersection_numbers_CY",
+    TEST_MODELS_PARAMETRISED,
+)
+def test_topdata_and_disc(
+    make_model,
+    expected_cy_dimension,
+    expected_no_divs,
+    expected_disc,
+    expected_intersection_numbers_CY,
+):
+    polytope = make_model()
+
+    discriminants = polytope.disc()
+    assert [(str(factor), codim) for factor, codim in discriminants] == expected_disc
+
+    polytope.topdata()
+    assert polytope.cy_dimension == expected_cy_dimension
+    assert polytope.no_divs == expected_no_divs
+    assert polytope.intersection_numbers_CY == expected_intersection_numbers_CY
Author	SHA1	Message	Date
Julian Piribauer	7217598d73	Implementing Frobenius basis for PFIdeals CI / lint (push) Successful in 1m20s Details CI / test (push) Successful in 3m49s Details	2026-11-01 16:46:13 +01:00
julianandJulian Piribauer	d3412f1ff6	Refactor and document operators and periods in sage module (#21 ) CI / lint (push) Successful in 52s Details CI / test (push) Successful in 7m2s Details --------- Co-authored-by: Julian Piribauer <julian.piribauer@gmail.com> Reviewed-on: #21	2026-10-31 22:22:38 +01:00
julianandJulian Piribauer	6600379681	Adding holo from l-vectors (#20 ) CI / lint (push) Successful in 24s Details CI / test (push) Successful in 2m32s Details --------- Co-authored-by: Julian Piribauer <julian.piribauer@gmail.com> Reviewed-on: #20	2026-10-27 14:43:37 +01:00
julianandJulian Piribauer	f5f960b42e	Implementing moving origin of ideal (#19 ) CI / lint (push) Successful in 20s Details CI / test (push) Successful in 1m59s Details --------- Co-authored-by: Julian Piribauer <julian.piribauer@gmail.com> Reviewed-on: #19	2026-09-09 22:03:29 +02:00
julianandJulian Piribauer	51a5612fe1	Initial setup for period computation (#13 ) CI / lint (push) Successful in 14s Details CI / test (push) Successful in 1m32s Details Introduces: - class for Picard--Fuchs operators and their ideals - class for periods (their complex linear combinations in the Frobenius bases) Implements: - method to obtain operator from period - method to get power series solution (at given indicials) to PF ideal --------- Co-authored-by: Julian Piribauer <julian.piribauer@gmail.com> Reviewed-on: #13	2026-08-17 19:51:31 +02:00
julianandJulian Piribauer	59dbb8f8bf	Updating README (#11 ) CI / lint (push) Canceled after 0s Details CI / test (push) Canceled after 0s Details --------- Co-authored-by: Julian Piribauer <julian.piribauer@gmail.com> Reviewed-on: #11	2026-07-29 20:30:12 +02:00
julianandJulian Piribauer	c6a2926ad6	5 pipeline with unit tests (#10 ) CI / lint (push) Successful in 12s Details CI / test (push) Successful in 1m30s Details --------- Co-authored-by: Julian Piribauer <julian.piribauer@gmail.com> Co-authored-by: julian <julian.piribauer@gmail.com> Reviewed-on: #10	2026-07-29 20:01:57 +02:00
julianandJulian Piribauer	bb830640fc	Adding class for projective ambient spaces (#9 ) --------- Co-authored-by: Julian Piribauer <julian.piribauer@gmail.com> Reviewed-on: #9	2026-07-26 13:22:08 +02:00