julian/Calabi-Yau-Period-Geometry
1
0
Fork 0
Code Issues 5 Pull Requests Actions Packages Projects Releases Wiki Activity

Implementing Frobenius basis for PFIdeals
CI / lint (push) Successful in 1m20s
Details
CI / test (push) Successful in 3m49s
Details

Browse Source
This commit is contained in:
Julian Piribauer
2026-11-01 16:46:13 +01:00
parent d3412f1ff6
commit 7217598d73
4 changed files with 895 additions and 110 deletions
Download Patch File Download Diff File Expand all files Collapse all files
sage/operators.sage
+509 -105
View File
@@ -8,20 +8,35 @@ PFOperator
simplify Group by z-monomials and factorise theta-polynomials
change_coordinates Rewrite the operator in new coordinates
_operator_from_string Parse a string into the operator polynomial
_theta_parts Theta-polynomial multiplying each z-monomial
PFOperatorAnsatz(PFOperator)
__init__ Ansatz of given theta- and z-multi-degrees
PFIdeal
__init__ Construct from a list of PFOperators
__init__ Construct from a list of PFOperators and optional Hodge numbers
add_operator Append a PFOperator
remove_operator Remove a PFOperator
change_coordinates Rewrite all operators in new coordinates
indicial_ideal Indicial ideal at z = 0, completed by linear algebra
indicials Indicials (exponents) at z = 0 with multiplicities
frobenius_basis Basis of all solutions at z = 0, logarithmic ones included
find_power_series_solution Power series solutions at given indicials/order
_operator_theta_parts Theta-parts of all operators in a common theta-ring
_complete_indicial_ideal Indicial ideal and the extra operators needed to see it
_solve_exponent_class Solutions with given indicials, up to a degree in log(z)
Helpers
_z_order_key Graded order on z-monomials (initial forms, solving order)
_left_multiply Normal-ordered product z^a * theta^c * operator
_indicial_roots Roots of an indicial ideal with their multiplicities
_leading_term Lowest term of a solution, for logging
"""
import copy
import logging
from sage.all import sage_eval, PolynomialRing, QQ, SR, function, prod, solve, var
from sage.all import (
sage_eval, PolynomialRing, QQ, SR, ZZ, binomial, block_matrix, floor, function, log, matrix, prod, solve, var
)
load("sage/util.py")
@@ -225,6 +240,24 @@ class PFOperator:
logger.debug("Changed coordinates of PFOperator to: %s", new_operator.operator_string)
return new_operator
def _theta_parts(self, theta_ring) -> dict:
"""
Splits this operator as sum_p z^p * P_p(theta) and returns {p: P_p}, each P_p an element of
theta_ring = QQ[theta0, ..., theta{n-1}]; z-monomials whose theta-polynomial vanishes are dropped.
theta_ring has to be built as PolynomialRing(QQ, n, names): without the explicit n, a single
variable gives a univariate ring, whose exponents() and monomial() work differently.
"""
if self.ring.base_ring() != QQ:
raise ValueError("Only operators with rational coefficients can be split into theta-parts.")
theta_gens = theta_ring.gens()
parts = {}
for coeff, monomial in self.operator:
exponents = monomial.exponents()[0]
z_exponents = tuple(exponents[:self.no_variables])
theta_monomial = prod(theta_gens[i] ** exponents[self.no_variables + i] for i in range(self.no_variables))
parts[z_exponents] = parts.get(z_exponents, theta_ring(0)) + coeff * theta_monomial
return {p: P for p, P in parts.items() if P != 0}
class PFOperatorAnsatz(PFOperator):
"""
A class for Picard-Fuchs operator Ansätze, characterised by number of variables and their z- and theta-multi-degrees.
@@ -261,12 +294,22 @@ class PFOperatorAnsatz(PFOperator):
class PFIdeal:
"""
A class representing a Picard-Fuchs ideal, which is a collection of PFOperators.
hodge_numbers optionally lists the numbers of solutions of top log-degree 0, 1, ..., n at a point of
maximal unipotent monodromy, i.e. h^{n,0}, ..., h^{0,n} of the (horizontal part of the) middle
cohomology - e.g. [1, 2, 2, 1] for a two-parameter Calabi-Yau threefold. Their sum is the holonomic
rank, the number of solutions at every point, and n bounds the degree in log(z) of every solution at
every point (monodromy theorem: N^(n+1) = 0 for a weight-n variation of Hodge structure). Without them,
indicials and frobenius_basis still work, just without the certificates and checks these allow.
"""
def __init__(self, operators: list):
def __init__(self, operators: list, hodge_numbers: list = None):
if not all(op.no_variables == operators[0].no_variables for op in operators):
raise ValueError("All operators must have the same number of variables.")
if hodge_numbers is not None and (not hodge_numbers or any(h not in ZZ or h < 1 for h in hodge_numbers)):
raise ValueError("hodge_numbers must be a non-empty list of positive integers, got %s." % (hodge_numbers,))
self.no_variables = operators[0].no_variables
self.operators = operators
self.hodge_numbers = list(hodge_numbers) if hodge_numbers is not None else None
logger.info("Initialised PFIdeal with %d operator(s).", len(self.operators))
def add_operator(self, pf_operator: PFOperator):
@@ -285,123 +328,484 @@ class PFIdeal:
"""
new_operators = [pf_operator.change_coordinates(*new_coordinates) for pf_operator in self.operators]
logger.info("Changed coordinates of PFIdeal with %d operator(s) to %s.", len(new_operators), new_coordinates)
return PFIdeal(new_operators)
return PFIdeal(new_operators, hodge_numbers=self.hodge_numbers)
def indicial_ideal(self, max_multiplier_degree: int = 4):
"""
Returns the indicial ideal of this PFIdeal at z = 0, an ideal of QQ[theta0, ..., theta{n-1}] whose
roots are the indicials (exponents) of the solutions there; see _complete_indicial_ideal.
"""
indicial_ideal, _ = self._complete_indicial_ideal(max_multiplier_degree)
return indicial_ideal
def indicials(self, max_multiplier_degree: int = 4) -> list:
"""
Returns the indicials (exponents) of this PFIdeal at z = 0, i.e. the roots of its indicial ideal, as
a sorted list of (rho, multiplicity) pairs, the multiplicities adding up to the holonomic rank. E.g.
the quintic has [((0,), 4)] at its MUM point and [((0,), 1), ((1,), 2), ((2,), 1)] at its conifold point.
"""
return _indicial_roots(self.indicial_ideal(max_multiplier_degree))
def frobenius_basis(self, order: int, max_log_degree: int = None, max_multiplier_degree: int = 4) -> list:
"""
Computes a basis of all solutions of this PFIdeal at z = 0, logarithmic ones included, as Periods
to the given order. For the basis at another point, move that point to the origin first with
change_coordinates.
The solutions are found one exponent class at a time: the indicials are grouped by their residue
mod Z^n, and each class is solved once, from its componentwise smallest exponent rho_min (which
becomes the Periods' indicials), so that solutions starting at rho_min + k with k > 0 - e.g. at the
indicials 0, 1, 1, 2 of a conifold point - come out of the same solve. Their degree in log(z) is at
most max_log_degree if given, and otherwise at most the class' total multiplicity minus 1 (the
logarithms of the monodromies' unipotent parts act on the class' solutions as commuting nilpotent
maps, so any product of that many of them vanishes) and at most n = len(hodge_numbers) - 1.
Within each class the basis is canonical (see _solve_exponent_class): every solution has its own
leading term with coefficient 1, at which the class' other solutions vanish. Coefficients are those
of plain log(z_i)^b, without factors of 2*pi*i or factorials.
This assumes the origin is a normal-crossing point of the discriminant, where every solution is a
sum of terms z^rho * (power series) * (polynomial in log z). Elsewhere fewer solutions than the
rank turn up, which is logged as a warning - as is a log structure at a MUM point that differs from
hodge_numbers.
"""
n = self.no_variables
indicial_ideal, operator_parts = self._complete_indicial_ideal(max_multiplier_degree)
theta_ring = indicial_ideal.ring()
roots = _indicial_roots(indicial_ideal)
classes = {}
for rho, multiplicity in roots:
classes.setdefault(tuple(r - floor(r) for r in rho), []).append((rho, multiplicity))
exponent_classes = sorted(
((tuple(min(rho[i] for rho, _ in members) for i in range(n)), members) for members in classes.values()),
key=lambda exponent_class: exponent_class[0],
)
cy_dimension = len(self.hodge_numbers) - 1 if self.hodge_numbers is not None else None
basis = []
for rho_min, members in exponent_classes:
class_size = sum(multiplicity for _, multiplicity in members)
# A solution starting at the indicial rho only turns up once the solve reaches z^(rho - rho_min).
required_order = max(sum(rho[i] - rho_min[i] for i in range(n)) for rho, _ in members)
if order < required_order:
raise ValueError(
"Order %d is too low to reach the indicials %s; it must be at least %d."
% (order, [rho for rho, _ in members], required_order)
)
if max_log_degree is not None:
log_degree = max_log_degree
elif cy_dimension is not None:
log_degree = min(class_size - 1, cy_dimension)
else:
log_degree = class_size - 1
solutions = self._solve_exponent_class(operator_parts, theta_ring, rho_min, order, log_degree)
logger.info(
"Exponent class of %s: %d solution(s) with leading terms %s.",
rho_min, len(solutions), [_leading_term(coefficients, rho_min) for coefficients in solutions],
)
if len(solutions) != class_size:
logger.warning(
"Found %d solution(s) in the exponent class of %s, whose indicials have total multiplicity "
"%d: z = 0 may not be a normal-crossing point of the discriminant, or the degree in log(z), "
"capped at %d, is too low.",
len(solutions), rho_min, class_size, log_degree,
)
basis.extend(
Period(no_variables=n, coefficients=coefficients, order=order, indicials=list(rho_min))
for coefficients in solutions
)
if self.hodge_numbers is not None:
if len(basis) != sum(self.hodge_numbers):
logger.warning(
"Found %d solution(s), but the rank sum(hodge_numbers) is %d.", len(basis), sum(self.hodge_numbers)
)
# At a point of maximal unipotent monodromy - all indicials 0, with a solution of top log-degree
# n - the numbers of solutions of each top log-degree are the Hodge numbers. A single exponent
# class at 0 isn't enough: at (0, -1/27) of P_{333111}[12], on the MUM divisor z0 = 0, the
# indicials are (0, 0) and (0, 1) and there is a log(z0)^4, but no MUM point.
top_log_degrees = [max(sum(log_index) for log_index in period.coefficients) for period in basis]
if (
len(roots) == 1
and all(r == 0 for r in roots[0][0])
and max(top_log_degrees, default=-1) == cy_dimension
):
log_structure = [top_log_degrees.count(degree) for degree in range(cy_dimension + 1)]
if log_structure != self.hodge_numbers:
logger.warning(
"The numbers of solutions per top log-degree %s at this MUM point differ from "
"hodge_numbers %s.", log_structure, self.hodge_numbers,
)
logger.info("Found a Frobenius basis of %d solution(s) at order %d.", len(basis), order)
return basis
def find_power_series_solution(self, indicials: list, order: int) -> list:
"""
Finds power series solutions (no logs) to this PFIdeal at given indicial exponents/order.
Finds the power series solutions (no logs) to this PFIdeal at given indicial exponents/order, i.e.
all solutions of the form z^indicials * (power series in z); frobenius_basis finds all solutions,
logarithmic ones included.
A term c*z^p*theta^q sends a_k*z^k (true exponent k+indicials) to
c*(k+indicials)^q*a_k*z^(k+p): it shifts index k up by p (p>=0), never down or sideways.
E.g. z0*theta0 sends a_k*z0^k to (k+rho0)*a_k*z0^(k+1).
This is the multivariate Frobenius method: canonical series solutions of a regular
holonomic D-ideal via its indicial ideal. See M. Saito, B. Sturmfels, N. Takayama,
"Gröbner Deformations of Hypergeometric Differential Equations", Algorithms and
Computation in Mathematics vol. 6, Springer, 2000, chs. 2-3.
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()
operator_terms = [
[
(monomial.exponents()[0][:self.no_variables], monomial.exponents()[0][self.no_variables:], coeff)
for coeff, monomial in pf_operator.operator
]
for pf_operator in self.operators
]
def eigenvalue(k, q):
# theta_i^q_i acts on z_i^(k_i + indicials[i]) as multiplication by
# (k_i + indicials[i])^q_i; theta^q's combined eigenvalue is the product over i.
value = QQ(1)
for i in range(self.no_variables):
if q[i]:
value *= (k[i] + indicials[i]) ** q[i]
return value
# Process multi-indices in order of increasing total degree: since every operator
# monomial has p >= 0 (componentwise), the z^m coefficient of L(y) only ever
# depends on a_k for k <= m, so by this point every k < m has already been solved.
z_indices = sorted(_multi_indices(self.no_variables, order), key=sum)
# solved[k] holds a_k written as a vector of coefficients over the `dimension`
# independent solutions found so far (a basis of the solution space up to k).
dimension = 0
solved = {}
for m in z_indices:
# For each operator, "coefficient of z^m in L(y) = 0" splits into a diagonal
# part (the p=0, theta-only monomials, whose unknown is a_m itself) plus a
# known part contributed by already-solved a_k with k = m - p, p > 0.
diagonals = []
known_contributions = []
for terms in operator_terms:
diagonal = QQ(0)
contribution = [QQ(0)] * dimension
for p, q, coeff in terms:
k = tuple(m[i] - p[i] for i in range(self.no_variables))
if any(ki < 0 for ki in k):
continue # this monomial would need a_k for a negative multi-index k: no such term
ev = eigenvalue(k, q)
if ev == 0:
continue # theta^q kills z_i^(k_i + indicials[i]) here, so this monomial contributes nothing
if k == m:
diagonal += coeff * ev # p = 0: coefficient multiplying the still-unknown a_m
else:
k_vector = solved[k] # p > 0: a_k is already known, add its contribution
for i in range(dimension):
contribution[i] += coeff * ev * k_vector[i]
diagonals.append(diagonal)
known_contributions.append(contribution)
# An operator with a nonzero diagonal lets us solve a_m = -(known part)/diagonal
# directly; this is exactly the indicial equation being nonzero at m + indicials.
active = next((r for r, d in enumerate(diagonals) if d != 0), None)
if active is not None:
value = [-known_contributions[active][i] / diagonals[active] for i in range(dimension)]
# Every operator's equation at m must independently be satisfied by this
# same a_m; disagreement means the ideal is inconsistent with these indicials.
for r, d in enumerate(diagonals):
if any(d * value[i] + known_contributions[r][i] != 0 for i in range(dimension)):
raise ValueError(
"Inconsistent Picard-Fuchs ideal or indicial exponents %s at multidegree %s."
% (indicials, m)
)
solved[m] = value
else:
# Resonance: every operator's indicial part vanishes at m + indicials, so a_m
# cannot be pinned down by this equation. If the already-known lower-degree
# data still forces a nonzero constraint here, satisfying it would require a
# log(z)-term solution, which this method (deliberately) does not compute.
if any(x != 0 for contribution in known_contributions for x in contribution):
raise ValueError(
"Resonance at multidegree %s for indicials %s would require a logarithmic "
"solution, which find_power_series_solution does not compute." % (m, indicials)
)
# Otherwise a_m is genuinely free: it starts a new independent solution, so
# extend every previously solved coefficient with a 0 in this new direction.
dimension += 1
for v in solved.values():
v.append(QQ(0))
solved[m] = [QQ(0)] * (dimension - 1) + [QQ(1)]
log_index = tuple([0] * self.no_variables)
solutions = [
Period(
no_variables=self.no_variables,
coefficients={log_index: {m: solved[m][i] for m in z_indices}},
order=order,
indicials=indicials,
)
for i in range(dimension)
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)