diff --git a/README.md b/README.md index 101220e..4f26c52 100644 --- a/README.md +++ b/README.md @@ -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") @@ -149,4 +149,43 @@ recovered[0].simplify().operator_string ``` `recovered[0]` is, up to scale, `M1` — recovered purely from `M1`, `M2`'s shared power series -solution. \ No newline at end of file +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. \ No newline at end of file diff --git a/sage/operators.sage b/sage/operators.sage index 34b6e2e..c56fe6f 100644 --- a/sage/operators.sage +++ b/sage/operators.sage @@ -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) diff --git a/tests/test_models.sage b/tests/test_models.sage index d69e8c1..15b04a6 100644 --- a/tests/test_models.sage +++ b/tests/test_models.sage @@ -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, [ diff --git a/tests/test_period_computation.py b/tests/test_period_computation.py index 8df924c..2016dce 100644 --- a/tests/test_period_computation.py +++ b/tests/test_period_computation.py @@ -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