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

Compare commits

merge into: julian/Calabi-Yau-Period-Geometry:master
Branches Tags
julian/Calabi-Yau-Period-Geometry:master julian/Calabi-Yau-Period-Geometry:frobenius-basis
...
pull from: julian/Calabi-Yau-Period-Geometry:frobenius-basis
Branches Tags
julian/Calabi-Yau-Period-Geometry:frobenius-basis julian/Calabi-Yau-Period-Geometry:master
1 Commits
master ... frobenius-basis
Author SHA1 Message Date
Julian Piribauer 7217598d73 Implementing Frobenius basis for PFIdeals
CI / test (push) Successful in 3m49s
Details
CI / lint (push) Successful in 1m20s
Details
2026-11-01 16:46:13 +01:00
4 changed files with 895 additions and 110 deletions
Expand all files Collapse all files
README.md
+40 -1
View File
@@ -8,7 +8,7 @@ It allows for computation of discriminant loci and topological data of manifolds
Supported initialisations are represented by the following examples.
```python
elliptic_curve_D = ToricPolytopeProjectiveSpace([1, 2, 3], model_name="elliptic_curve_D")
elliptic_curve_D = ToricPolytopeProjectiveSpace([3, 2, 1], model_name="elliptic_curve_D")
CY3_quintic = ToricPolytopeProjectiveSpace([1, 1, 1, 1, 1], model_name="quintic")
CY3_bicubic = ToricPolytopeCICY([[3, 3]], model_name="bi-cubic")
@@ -150,3 +150,42 @@ recovered[0].simplify().operator_string
`recovered[0]` is, up to scale, `M1` — recovered purely from `M1`, `M2`'s shared power series
solution.
### Frobenius bases
A `PFIdeal` also computes its complete Frobenius basis — all solutions, logarithmic ones included —
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` — 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 — 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πi.
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)
tests/test_models.sage
+3 -3
View File
@@ -5,7 +5,7 @@
TEST_MODELS = [
(
"K3_two_parameter_family",
lambda: ToricPolytopeProjectiveSpace([1, 1, 2, 4], model_name="K3_two_parameter_family"),
lambda: ToricPolytopeProjectiveSpace([4, 2, 1, 1], model_name="K3_two_parameter_family"),
2, # cy_dimension
2, # no_divs
[("z2 - 1/4", 2), ("4096*z1^2*(4*z2 - 1) + 128*z1 - 1", 0)],
@@ -13,7 +13,7 @@ TEST_MODELS = [
),
(
"CY3_two_parameter_family",
lambda: ToricPolytopeProjectiveSpace([1, 1, 2, 2, 2], model_name="CY3_two_parameter_family"),
lambda: ToricPolytopeProjectiveSpace([2, 2, 2, 1, 1], model_name="CY3_two_parameter_family"),
3,
2,
[("z2 - 1/4", 3), ("65536*z1^2*(4*z2 - 1) + 512*z1 - 1", 0)],
@@ -21,7 +21,7 @@ TEST_MODELS = [
),
(
"CY4_two_parameter_family",
lambda: ToricPolytopeProjectiveSpace([1, 1, 1, 1, 8, 12], model_name="CY4_two_parameter_family"),
lambda: ToricPolytopeProjectiveSpace([12, 8, 1, 1, 1, 1], model_name="CY4_two_parameter_family"),
4,
2,
[
tests/test_period_computation.py
+342
View File
@@ -1,10 +1,108 @@
import logging
import os
import pytest
from sage.all import * # noqa: F401
load(os.path.join(os.path.dirname(__file__), "..", "sage", "periods.sage"))
def _quintic_ideal():
# The quintic's Picard-Fuchs operator, of rank 4 with log structure [1, 1, 1, 1] at its MUM point.
return PFIdeal(
[PFOperator("theta0^4 - 5*z0*(5*theta0 + 1)*(5*theta0 + 2)*(5*theta0 + 3)*(5*theta0 + 4)", no_variables=1)],
hodge_numbers=[1, 1, 1, 1],
)
def _p22211_ideal():
# The Picard-Fuchs ideal of P_{22211}[8] at its MUM point (see
# test_ideal_finds_power_series_solution_and_recovers_first_operator), with h^{2,1} = 2.
D1 = PFOperator(
"theta0^2*(theta0 - 2*theta1) - 4*z0*(4*theta0 + 3)*(4*theta0 + 2)*(4*theta0 + 1)",
no_variables=2,
)
D2 = PFOperator(
"theta1^2 - z1*(2*theta1 - theta0 + 1)*(2*theta1 - theta0)",
no_variables=2,
)
return PFIdeal([D1, D2], hodge_numbers=[1, 2, 2, 1])
def _p333111_ideal():
# The Picard-Fuchs ideal of the fourfold P_{333111}[12] at its MUM point: the GKZ operators of
# l^(1) = (-4; 1,1,1,0,0,0,1) and l^(2) = (0; 0,0,0,1,1,1,-3), the first without its left factor theta0.
# With the intersection numbers kappa_0000 = 36, kappa_0001 = 12, kappa_0011 = 4 (and 0 else) that
# toric_topdata computes, the vertical h^{2,2} is 3.
D1 = PFOperator(
"theta0^2*(theta0 - 3*theta1) - 4*z0*(4*theta0 + 1)*(4*theta0 + 2)*(4*theta0 + 3)",
no_variables=2,
)
D2 = PFOperator(
"theta1^3 - z1*(theta0 - 3*theta1)*(theta0 - 3*theta1 - 1)*(theta0 - 3*theta1 - 2)",
no_variables=2,
)
return PFIdeal([D1, D2], hodge_numbers=[1, 2, 3, 2, 1])
def _starting_term(period):
# A solution's lowest z-monomial (by total degree, then lexicographically) together with its coefficient
# per log-monomial there, e.g. ((1,), {(1,): 1}) for z0*log(z0) + ...
z_index = min(
(z_index for z_dict in period.coefficients.values() for z_index, coeff in z_dict.items() if coeff != 0),
key=lambda z_index: (sum(z_index), z_index),
)
return z_index, {
log_index: z_dict[z_index] for log_index, z_dict in period.coefficients.items() if z_dict[z_index] != 0
}
def _nonzero(coefficients):
# A coefficients dictionary without its zero entries, to compare bases that store zeros differently.
return {
log_index: {z_index: coeff for z_index, coeff in z_dict.items() if coeff != 0}
for log_index, z_dict in coefficients.items()
if any(coeff != 0 for coeff in z_dict.values())
}
def _truncated(f, degree):
# A polynomial without its terms of total degree above degree.
return sum((coeff * monomial for coeff, monomial in f if monomial.degree() <= degree), f.parent()(0))
def _inverse_of_linear(constant, linear, degree):
# 1/(constant + linear) for a linear form in eps, as a polynomial truncated at the given degree.
return sum(((-linear / constant) ** j for j in range(degree + 1)), linear.parent()(0)) / constant
def _classical_frobenius_periods(coefficient_ratio, starting_terms, order):
# The classical Frobenius method. With coefficient_ratio(k) = c(k + eps)/c(eps) for the Gamma-function
# coefficients c of the fundamental period (a polynomial in eps, truncated at the starting terms' degree),
# the period starting with Q(log z) is Q(d/deps) applied to sum_k c(k + eps)/c(eps) * z^(k + eps) at eps = 0.
n = starting_terms[0].parent().ngens()
log_ring = PolynomialRing(QQ, n, ["L%d" % i for i in range(n)])
logs = log_ring.gens()
periods = [{} for _ in starting_terms]
for k in _multi_indices(n, order):
ratio = coefficient_ratio(k)
for period, starting_term in zip(periods, starting_terms):
# d^alpha/deps^alpha of ratio * exp(eps . L) at eps = 0, summed over the monomials eps^alpha of Q
value = log_ring(0)
for q, q_monomial in starting_term:
alpha = q_monomial.exponents()[0]
for r, r_monomial in ratio:
beta = r_monomial.exponents()[0]
if all(beta[i] <= alpha[i] for i in range(n)):
value += q * r * prod(
factorial(alpha[i]) // factorial(alpha[i] - beta[i]) * logs[i] ** (alpha[i] - beta[i])
for i in range(n)
)
for coeff, log_monomial in value:
period.setdefault(tuple(log_monomial.exponents()[0]), {})[k] = coeff
return periods
def test_find_power_series_solution_recovers_quintic_period():
# The quintic's Picard-Fuchs operator has indicial equation theta0^4 = 0 at z0 = 0
# (a quadruple root at 0), so its holomorphic (non-logarithmic) power series solution
@@ -103,3 +201,247 @@ def test_single_operator_recovers_itself_from_its_holomorphic_period():
assert ratio != 0
assert recovered[0].operator == ratio * L.operator
def test_indicials_of_quintic_at_mum_conifold_and_orbifold_points():
# At the MUM point all four indicials are 0; at the conifold point they are 0, 1, 1, 2, the doubled 1
# carrying the single logarithmic solution; at the orbifold point z0 -> 1/z0 they are 1/5, ..., 4/5.
ideal = _quintic_ideal()
assert ideal.indicials() == [((0,), 4)]
assert ideal.change_coordinates("1 - 3125*z0").indicials() == [((0,), 1), ((1,), 2), ((2,), 1)]
assert ideal.change_coordinates("1/z0").indicials() == [((QQ(k) / 5,), 1) for k in range(1, 5)]
def test_indicials_reject_hodge_numbers_inconsistent_with_the_operators():
# The quintic's operator has rank 4, so its indicials can't add up to the rank 5 of these Hodge numbers.
ideal = PFIdeal(_quintic_ideal().operators, hodge_numbers=[1, 1, 1, 1, 1])
with pytest.raises(ValueError):
ideal.indicials()
def test_frobenius_basis_of_quintic_at_mum_point():
# The solutions start with 1, log(z0), log(z0)^2 and log(z0)^3, the first being the classic holomorphic
# period. Each solution's top-log coefficient is a solution without logs, hence that period again.
ideal = _quintic_ideal()
basis = ideal.frobenius_basis(order=4)
assert [_starting_term(period) for period in basis] == [((0,), {(j,): 1}) for j in range(4)]
holomorphic = {(0,): 1, (1,): 120, (2,): 113400, (3,): 168168000, (4,): 305540235000}
assert basis[0].coefficients == {(0,): holomorphic}
assert all(period.coefficients[(j,)] == holomorphic for j, period in enumerate(basis))
assert all(period.apply_operator(op).coefficients == {} for period in basis for op in ideal.operators)
def test_frobenius_basis_of_quintic_at_conifold_point():
# The indicials 0, 1, 1, 2 differ by integers, so all four solutions come out of one solve from exponent 0.
# Picard-Lefschetz: exactly one of them is logarithmic, z0*log(z0) + ..., starting where the vanishing
# period z0 + ... does.
conifold = _quintic_ideal().change_coordinates("1 - 3125*z0")
basis = conifold.frobenius_basis(order=4)
assert [_starting_term(period) for period in basis] == [
((0,), {(0,): 1}),
((1,), {(0,): 1}),
((1,), {(1,): 1}),
((2,), {(0,): 1}),
]
assert [max(log_index[0] for log_index in period.coefficients) for period in basis] == [0, 0, 1, 0]
assert all(period.apply_operator(op).coefficients == {} for period in basis for op in conifold.operators)
def test_frobenius_basis_of_quintic_at_orbifold_point():
# Four exponent classes 1/5, ..., 4/5 mod 1, each with a single solution without logs. Order 6 reaches
# the step at which Sage's echelon form of sparse rational matrices went wrong (see operators.sage).
orbifold = _quintic_ideal().change_coordinates("1/z0")
basis = orbifold.frobenius_basis(order=6)
assert [period.indicials for period in basis] == [[QQ(k) / 5] for k in range(1, 5)]
assert all(list(period.coefficients) == [(0,)] for period in basis)
assert all(period.apply_operator(op).coefficients == {} for period in basis for op in orbifold.operators)
def test_frobenius_basis_of_two_parameter_model_at_mum_point():
# All six indicials are (0, 0), and the solutions start with 1; log(z0), log(z1); log(z0)^2,
# log(z0)*log(z1); and a triple-log term proportional to kappa_ijk*log(z_i)*log(z_j)*log(z_k) with the
# intersection numbers kappa_000 = 8, kappa_001 = 4: the log structure [1, 2, 2, 1] of the Hodge numbers.
ideal = _p22211_ideal()
assert ideal.indicials() == [((0, 0), 6)]
basis = ideal.frobenius_basis(order=3)
assert [_starting_term(period) for period in basis] == [
((0, 0), {(0, 0): 1}),
((0, 0), {(1, 0): 1}),
((0, 0), {(0, 1): 1}),
((0, 0), {(2, 0): 1}),
((0, 0), {(1, 1): 1}),
((0, 0), {(3, 0): 1, (2, 1): QQ(3) / 2}),
]
assert all(period.apply_operator(op).coefficients == {} for period in basis for op in ideal.operators)
def test_frobenius_basis_of_two_parameter_model_at_strong_coupling():
# At (z0, z1) = (0, 1/4) the generators' z-free parts theta0^2*theta1 and theta1*(2*theta1 - 2*theta0 - 1)
# vanish on the whole line theta1 = 0. Completing the indicial ideal pins the indicials down to (0, 0)
# with multiplicity 4 and (0, 1/2) with multiplicity 2. The latter class holds sqrt(z1)*(1 + ...), the
# power series solution of test_ideal_finds_power_series_solution_and_recovers_first_operator, and a
# logarithmic partner sqrt(z1)*(log(z0) + log(z1) + ...).
moved = _p22211_ideal().change_coordinates("z0", "z1 - 1/4")
assert moved.indicials() == [((0, 0), 4), ((0, QQ(1) / 2), 2)]
basis = moved.frobenius_basis(order=4)
assert [period.indicials for period in basis] == [[0, 0]] * 4 + [[0, QQ(1) / 2]] * 2
assert [_starting_term(period) for period in basis] == [
((0, 0), {(0, 0): 1}),
((0, 0), {(1, 0): 1}),
((0, 0), {(2, 0): 1}),
((0, 0), {(3, 0): 1}),
((0, 0), {(0, 0): 1}),
((0, 0), {(1, 0): 1, (0, 1): 1}),
]
power_series = moved.find_power_series_solution(indicials=[0, QQ(1) / 2], order=4)
assert basis[4].coefficients == power_series[0].coefficients
assert all(period.apply_operator(op).coefficients == {} for period in basis for op in moved.operators)
def test_frobenius_basis_warns_when_log_degree_is_capped_too_low(caplog):
# The triple-log solution at the MUM point needs log-degree 3, so with max_log_degree=2 only five of the
# six solutions survive - which is reported rather than silently returned.
ideal = _p22211_ideal()
with caplog.at_level(logging.WARNING):
basis = ideal.frobenius_basis(order=2, max_log_degree=2)
assert len(basis) == 5
assert any("rank sum(hodge_numbers) is 6" in record.getMessage() for record in caplog.records)
def test_frobenius_basis_of_two_parameter_fourfold(caplog):
# At the MUM point the nine solutions have the log structure [1, 2, 3, 2, 1], and the top ones are given by
# the intersection numbers: the quartic one by kappa_ijkl*log(z_i)*...*log(z_l), i.e. 36*log(z0)^4 +
# 48*log(z0)^3*log(z1) + 24*log(z0)^2*log(z1)^2 up to scale, a cubic one by its derivative in log(z1).
ideal = _p333111_ideal()
assert ideal.indicials() == [((0, 0), 9)]
basis = ideal.frobenius_basis(order=3)
assert [_starting_term(period) for period in basis] == [
((0, 0), {(0, 0): 1}),
((0, 0), {(1, 0): 1}),
((0, 0), {(0, 1): 1}),
((0, 0), {(2, 0): 1}),
((0, 0), {(1, 1): 1}),
((0, 0), {(0, 2): 1}),
((0, 0), {(3, 0): 1, (2, 1): 1}),
((0, 0), {(1, 2): 1}),
((0, 0), {(4, 0): 1, (3, 1): QQ(4) / 3, (2, 2): QQ(2) / 3}),
]
assert all(period.apply_operator(op).coefficients == {} for period in basis for op in ideal.operators)
# At (0, -1/27), where the strong coupling divisor 1 + 27*z1 = 0 crosses the MUM divisor, all indicials lie
# in one exponent class and there is a log(z0)^4 solution. But some indicials are (0, 1): this is no MUM
# point, and its log structure [2, 3, 2, 1, 1] must not be compared with the Hodge numbers.
moved = ideal.change_coordinates("z0", "z1 + 1/27")
with caplog.at_level(logging.WARNING):
basis = moved.frobenius_basis(order=3)
assert moved.indicials() == [((0, 0), 5), ((0, 1), 4)]
assert len(basis) == 9
assert all(period.apply_operator(op).coefficients == {} for period in basis for op in moved.operators)
assert not caplog.records
def test_frobenius_bases_at_mum_points_agree_with_the_classical_frobenius_method():
# An independent construction of every period at the MUM points, logarithmic ones included: the classical
# Frobenius method, from nothing but the Gamma-function coefficients c(k) of the fundamental period. The
# starting terms are those of the canonical bases; for P_{22211}[8] the cubic one is
# kappa_ijk*eps_i*eps_j*eps_k/8 with the intersection numbers kappa_000 = 8, kappa_001 = 4.
eps_ring = PolynomialRing(QQ, 1, ["e0"])
e = eps_ring.gen(0)
def quintic_ratio(k):
# c(n + eps)/c(eps) for c(n) = (5n)!/(n!)^5
ratio = eps_ring(1)
for i in range(1, 5 * k[0] + 1):
ratio = _truncated(ratio * (5 * e + i), 3)
for i in range(1, k[0] + 1):
ratio = _truncated(ratio * _inverse_of_linear(i, e, 3) ** 5, 3)
return ratio
classical = _classical_frobenius_periods(quintic_ratio, [e**j for j in range(4)], order=8)
basis = _quintic_ideal().frobenius_basis(order=8)
assert [_nonzero(period.coefficients) for period in basis] == [_nonzero(c) for c in classical]
eps_ring = PolynomialRing(QQ, 2, ["e0", "e1"])
e0, e1 = eps_ring.gens()
def p22211_ratio(k):
# c(n + eps0, m + eps1)/c(eps0, eps1) for c(n, m) = (4n)!/((n!)^3 (m!)^2 (n - 2m)!). For n < 2m the
# last factor is Gamma(1 + delta)/Gamma(1 + delta + n - 2m) with delta = eps0 - 2*eps1, a polynomial
# vanishing at eps = 0, which only the logarithmic periods see.
n, m = k
delta = e0 - 2 * e1
ratio = eps_ring(1)
for i in range(1, 4 * n + 1):
ratio = _truncated(ratio * (4 * e0 + i), 3)
for i in range(1, n + 1):
ratio = _truncated(ratio * _inverse_of_linear(i, e0, 3) ** 3, 3)
for i in range(1, m + 1):
ratio = _truncated(ratio * _inverse_of_linear(i, e1, 3) ** 2, 3)
for i in range(1, n - 2 * m + 1):
ratio = _truncated(ratio * _inverse_of_linear(i, delta, 3), 3)
for i in range(2 * m - n):
ratio = _truncated(ratio * (delta - i), 3)
return ratio
starting_terms = [eps_ring(1), e0, e1, e0**2, e0 * e1, e0**3 + QQ(3) / 2 * e0**2 * e1]
classical = _classical_frobenius_periods(p22211_ratio, starting_terms, order=6)
basis = _p22211_ideal().frobenius_basis(order=6)
assert [_nonzero(period.coefficients) for period in basis] == [_nonzero(c) for c in classical]
def test_quintic_bases_at_mum_and_conifold_point_continue_into_each_other():
# Both expansions converge on the lens between z0 = 0 and z0 = 1/3125, where the MUM periods must be
# constant linear combinations T of the conifold periods. Through T the conifold monodromy
# (log w -> log w + 2*pi*i) becomes a Picard-Lefschetz transvection in the MUM basis, and its product with
# the MUM monodromy is the monodromy around z0 = oo, of order 5.
CF = ComplexField(200)
two_pi_i = 2 * CF.pi() * CF.gen()
ideal = _quintic_ideal()
mum = ideal.frobenius_basis(order=100)
conifold = ideal.change_coordinates("1 - 3125*z0").frobenius_basis(order=100)
def values(basis, points, shifted=False):
# One row per point z0, one column per period, continued once around w = 0 if shifted.
rows = []
for z0 in points:
w = z0 if basis is mum else 1 - 3125 * z0
log_w = CF(w).log() + (two_pi_i if shifted else 0)
rows.append([
sum(
CF(coeff) * log_w ** log_index[0] * ((z_index[0] + period.indicials[0]) * log_w).exp()
for log_index, z_dict in period.coefficients.items()
for z_index, coeff in z_dict.items()
if coeff != 0
)
for period in basis
])
return matrix(CF, rows)
fit = [CF(x, y) / 3125 for x, y in [(0.4, -0.15), (0.5, 0.1), (0.6, -0.1), (0.45, 0.15)]]
test = [CF(x, y) / 3125 for x, y in [(0.55, 0.05), (0.42, -0.05), (0.6, 0.12)]]
T = values(conifold, fit).solve_right(values(mum, fit)).transpose()
assert (values(mum, test) - values(conifold, test) * T.transpose()).norm() < 1e-8 * values(mum, test).norm()
M_0 = values(mum, fit).solve_right(values(mum, fit, shifted=True)).transpose()
M_c = T * values(conifold, fit).solve_right(values(conifold, fit, shifted=True)).transpose() * T.inverse()
identity = identity_matrix(CF, 4)
singular_values = (M_c - identity).change_ring(CDF).SVD()[1].diagonal()
assert singular_values[1] < 1e-8 * singular_values[0]
assert ((M_0 * M_c) ** 5 - identity).norm() < 1e-8