""" Picard-Fuchs operators and ideals in the logarithmic derivatives theta_i = z_i * d/dz_i. Contents -------- PFOperator __init__ Construct from a string in z- and theta-variables simplify Group by z-monomials and factorise theta-polynomials change_coordinates Rewrite the operator in new coordinates _operator_from_string Parse a string into the operator polynomial _theta_parts Theta-polynomial multiplying each z-monomial PFOperatorAnsatz(PFOperator) __init__ Ansatz of given theta- and z-multi-degrees PFIdeal __init__ Construct from a list of PFOperators and optional Hodge numbers add_operator Append a PFOperator remove_operator Remove a PFOperator change_coordinates Rewrite all operators in new coordinates indicial_ideal Indicial ideal at z = 0, completed by linear algebra indicials Indicials (exponents) at z = 0 with multiplicities frobenius_basis Basis of all solutions at z = 0, logarithmic ones included find_power_series_solution Power series solutions at given indicials/order _operator_theta_parts Theta-parts of all operators in a common theta-ring _complete_indicial_ideal Indicial ideal and the extra operators needed to see it _solve_exponent_class Solutions with given indicials, up to a degree in log(z) Helpers _z_order_key Graded order on z-monomials (initial forms, solving order) _left_multiply Normal-ordered product z^a * theta^c * operator _indicial_roots Roots of an indicial ideal with their multiplicities _leading_term Lowest term of a solution, for logging """ import copy import logging from sage.all import ( sage_eval, PolynomialRing, QQ, SR, ZZ, binomial, block_matrix, floor, function, log, matrix, prod, solve, var ) load("sage/util.py") # Logger logger = logging.getLogger(__name__) logging.basicConfig(level=logging.DEBUG) class PFOperator: """ A class representing a Picard-Fuchs operator in a given number of variables (moduli). Independent of the given order, the variables are assumed to be left of the derivatives. The variables extra_locals are needed for symbolic coefficients used for operator Ansätze. """ def _operator_from_string(self, operator_string: str, extra_locals: dict = None): z_names = ['z%d' % i for i in range(self.no_variables)] theta_names = ['theta%d' % i for i in range(self.no_variables)] base_ring = SR if extra_locals else QQ self.ring = PolynomialRing(base_ring, z_names + theta_names) self.z_gens = self.ring.gens()[:self.no_variables] self.theta_gens = self.ring.gens()[self.no_variables:] parse_locals = dict(self.ring.gens_dict()) if extra_locals: parse_locals.update(extra_locals) try: operator = self.ring(sage_eval(operator_string, locals=parse_locals)) return operator except Exception as e: raise ValueError("Invalid operator string: %s" % e) def __init__(self, operator_string: str, no_variables: int = 1, extra_locals: dict = None): self.no_variables = no_variables self.operator_string = operator_string self.operator = self._operator_from_string(operator_string, extra_locals=extra_locals) logger.info("Initialised PFOperator: %s", self.operator) def simplify(self) -> "PFOperator": """ Groups the operator's terms by z-monomial and factorises the theta-polynomial multiplying each z-monomial, e.g. turning a computed quintic operator into the well-known theta0^4 - 5*z0*(5*theta0 + 1)*(5*theta0 + 2)*(5*theta0 + 3)*(5*theta0 + 4). """ z_monomial_theta_parts = {} for coeff, monomial in self.operator: exponents = monomial.exponents()[0] z_exponents = exponents[:self.no_variables] theta_exponents = exponents[self.no_variables:] theta_monomial = SR(1) for gen, exp in zip(self.theta_gens, theta_exponents): theta_monomial *= SR(gen) ** exp z_monomial_theta_parts[z_exponents] = z_monomial_theta_parts.get(z_exponents, SR(0)) + SR(coeff) * theta_monomial simplified_expr = SR(0) for z_exponents, theta_part in z_monomial_theta_parts.items(): z_monomial = SR(1) for gen, exp in zip(self.z_gens, z_exponents): z_monomial *= SR(gen) ** exp simplified_expr += theta_part.factor() * z_monomial simplified = copy.copy(self) simplified.operator_string = str(simplified_expr) logger.debug("Simplified PFOperator to: %s", simplified.operator_string) return simplified def change_coordinates(self, *new_coordinates: str) -> "PFOperator": """ Rewrites this operator in new coordinates w0, ..., w{n-1} := new_coordinates(z0, ..., z{n-1}), where new_coordinates is given as n expressions in the *old* coordinates z0, ..., z{n-1}; the inverse change of coordinates needed to do this is solved for automatically (the two aren't independent - one determines the other). The new origin w = 0 is the point z0, ..., z{n-1} at which new_coordinates vanishes. E.g. pf_operator.change_coordinates("1 - 5*z0", "1/z1") moves the origin to z0 = 1/5, z1 = infinity. """ n = self.no_variables if len(new_coordinates) != n: raise ValueError( "Expected %d new coordinate expression(s), got %d." % (n, len(new_coordinates)) ) z_gens = [var('z%d' % i) for i in range(n)] parse_locals = {str(z): z for z in z_gens} try: new_coords = [SR(sage_eval(expr, locals=parse_locals)) for expr in new_coordinates] except Exception as e: raise ValueError("Invalid coordinate expression: %s" % e) w_gens = [var('w%d' % i) for i in range(n)] equations = [w_gens[i] == new_coords[i] for i in range(n)] solutions = solve(equations, z_gens, solution_dict=True) if not solutions: raise ValueError("Could not invert the given change of coordinates %s." % (new_coordinates,)) inverse_coords = [solutions[0][z] for z in z_gens] Y = function('Y')(*w_gens) y = Y.subs({w_gens[i]: new_coords[i] for i in range(n)}) def apply_theta(expr, i): return z_gens[i] * expr.diff(z_gens[i]) total = SR(0) for coeff, monomial in self.operator: exponents = monomial.exponents()[0] z_exponents = exponents[:n] theta_exponents = exponents[n:] term = y for i, power in enumerate(theta_exponents): for _ in range(power): term = apply_theta(term, i) z_monomial = prod(z_gens[i] ** z_exponents[i] for i in range(n)) total += SR(coeff) * z_monomial * term total = total.subs({z_gens[i]: inverse_coords[i] for i in range(n)}) total = total.simplify_full().expand() total_theta_degree = max(sum(monomial.exponents()[0][n:]) for _, monomial in self.operator) theta_ring = PolynomialRing(QQ, n, ['theta%d' % i for i in range(n)]) theta_gens = theta_ring.gens() # Peel off the coefficient of each derivative order k, highest total order first (so # a coefficient can never still contain an as-yet-unextracted higher derivative as a # factor), and rewrite w^k*D^k Y as a falling-factorial polynomial in theta applied # to Y - the Euler-operator identity theta*(theta-1)*...*(theta-k+1) = w^k*d^k/dw^k, # taken variable by variable. Each term's coefficient is left as a general rational # function of w for now; only once every term has been collected is the whole # operator scaled by their common denominator (see below). remainder = total terms = [] for k in sorted(_multi_indices(n, total_theta_degree), key=lambda k: -sum(k)): if sum(k) == 0: dterm = Y else: args = [arg for i in range(n) if k[i] for arg in (w_gens[i], k[i])] dterm = Y.diff(*args) coeff = remainder.coefficient(dterm) if coeff == 0: continue remainder -= coeff * dterm falling_factorial = theta_ring(1) for i in range(n): for j in range(k[i]): falling_factorial *= (theta_gens[i] - j) w_monomial = prod(w_gens[i] ** k[i] for i in range(n)) terms.append(((coeff / w_monomial).simplify_rational(), falling_factorial)) remainder = remainder.simplify_full() if remainder != 0: raise ValueError("Residual nonzero after change of coordinates: %s" % remainder) if not terms: raise ValueError("Change of coordinates produced the zero operator.") w_ring = PolynomialRing(QQ, n, ['w%d' % i for i in range(n)]) common_denominator = w_ring(1) for rational_coeff, _ in terms: common_denominator = common_denominator.lcm(w_ring(rational_coeff.denominator())) result_terms = {} for rational_coeff, falling_factorial in terms: scaled = w_ring((rational_coeff * SR(common_denominator)).simplify_rational()) for w_coeff, w_monomial in scaled: p = w_monomial.exponents()[0] for theta_coeff, theta_monomial in falling_factorial: q = theta_monomial.exponents()[0] key = (p, q) result_terms[key] = result_terms.get(key, QQ(0)) + w_coeff * theta_coeff content = gcd([QQ(coeff) for coeff in result_terms.values() if coeff != 0]) if content not in (0, 1): result_terms = {key: coeff / content for key, coeff in result_terms.items()} def _format_monomial(coeff, p, q): factors = [str(QQ(coeff))] for i in range(n): if p[i]: factors.append("z%d^%d" % (i, p[i]) if p[i] != 1 else "z%d" % i) for i in range(n): if q[i]: factors.append("theta%d^%d" % (i, q[i]) if q[i] != 1 else "theta%d" % i) return "*".join(factors) terms_string = [ _format_monomial(coeff, p, q) for (p, q), coeff in result_terms.items() if coeff != 0 ] if not terms_string: raise ValueError("Change of coordinates produced the zero operator.") operator_string = " + ".join(terms_string) new_operator = PFOperator(operator_string, no_variables=n) logger.debug("Changed coordinates of PFOperator to: %s", new_operator.operator_string) return new_operator def _theta_parts(self, theta_ring) -> dict: """ Splits this operator as sum_p z^p * P_p(theta) and returns {p: P_p}, each P_p an element of theta_ring = QQ[theta0, ..., theta{n-1}]; z-monomials whose theta-polynomial vanishes are dropped. theta_ring has to be built as PolynomialRing(QQ, n, names): without the explicit n, a single variable gives a univariate ring, whose exponents() and monomial() work differently. """ if self.ring.base_ring() != QQ: raise ValueError("Only operators with rational coefficients can be split into theta-parts.") theta_gens = theta_ring.gens() parts = {} for coeff, monomial in self.operator: exponents = monomial.exponents()[0] z_exponents = tuple(exponents[:self.no_variables]) theta_monomial = prod(theta_gens[i] ** exponents[self.no_variables + i] for i in range(self.no_variables)) parts[z_exponents] = parts.get(z_exponents, theta_ring(0)) + coeff * theta_monomial return {p: P for p, P in parts.items() if P != 0} class PFOperatorAnsatz(PFOperator): """ A class for Picard-Fuchs operator Ansätze, characterised by number of variables and their z- and theta-multi-degrees. """ def __init__(self, no_variables: int, theta_degree: int, z_degree: int): self.no_variables = no_variables self.theta_degree = theta_degree self.z_degree = z_degree z_indices = _multi_indices(self.no_variables, z_degree) theta_indices = _multi_indices(self.no_variables, theta_degree) # Every (z_index, theta_index) pair allowed by the degree bounds gets its own # fresh unknown, to be solved for once the Ansatz is applied to a period. operator_terms = [ (var("b_" + "_".join(str(x) for x in z_index + theta_index)), (z_index, theta_index)) for z_index in z_indices for theta_index in theta_indices ] self.unknowns = [coeff for coeff, _ in operator_terms] def _monomial_factors(index, name): return [f"{name}{i}^{exp}" for i, exp in enumerate(index) if exp > 0] operator_string = " + ".join( "*".join([str(coeff), *_monomial_factors(z_index, "z"), *_monomial_factors(theta_index, "theta")]) for coeff, (z_index, theta_index) in operator_terms ) extra_locals = {str(coeff): coeff for coeff in self.unknowns} super().__init__(operator_string=operator_string, no_variables=no_variables, extra_locals=extra_locals) class PFIdeal: """ A class representing a Picard-Fuchs ideal, which is a collection of PFOperators. hodge_numbers optionally lists the numbers of solutions of top log-degree 0, 1, ..., n at a point of maximal unipotent monodromy, i.e. h^{n,0}, ..., h^{0,n} of the (horizontal part of the) middle cohomology - e.g. [1, 2, 2, 1] for a two-parameter Calabi-Yau threefold. Their sum is the holonomic rank, the number of solutions at every point, and n bounds the degree in log(z) of every solution at every point (monodromy theorem: N^(n+1) = 0 for a weight-n variation of Hodge structure). Without them, indicials and frobenius_basis still work, just without the certificates and checks these allow. """ def __init__(self, operators: list, hodge_numbers: list = None): if not all(op.no_variables == operators[0].no_variables for op in operators): raise ValueError("All operators must have the same number of variables.") if hodge_numbers is not None and (not hodge_numbers or any(h not in ZZ or h < 1 for h in hodge_numbers)): raise ValueError("hodge_numbers must be a non-empty list of positive integers, got %s." % (hodge_numbers,)) self.no_variables = operators[0].no_variables self.operators = operators self.hodge_numbers = list(hodge_numbers) if hodge_numbers is not None else None logger.info("Initialised PFIdeal with %d operator(s).", len(self.operators)) def add_operator(self, pf_operator: PFOperator): self.operators.append(pf_operator) logger.info("Added PFOperator to PFIdeal: %s", pf_operator.operator_string) def remove_operator(self, pf_operator: PFOperator): self.operators.remove(pf_operator) logger.info("Removed PFOperator from PFIdeal: %s", pf_operator.operator_string) def change_coordinates(self, *new_coordinates: str) -> "PFIdeal": """ Moves the origin of the whole system of differential equations to a new point, by returning a new PFIdeal whose operators are each rewritten in new_coordinates; see PFOperator.change_coordinates for what new_coordinates means. """ new_operators = [pf_operator.change_coordinates(*new_coordinates) for pf_operator in self.operators] logger.info("Changed coordinates of PFIdeal with %d operator(s) to %s.", len(new_operators), new_coordinates) return PFIdeal(new_operators, hodge_numbers=self.hodge_numbers) def indicial_ideal(self, max_multiplier_degree: int = 4): """ Returns the indicial ideal of this PFIdeal at z = 0, an ideal of QQ[theta0, ..., theta{n-1}] whose roots are the indicials (exponents) of the solutions there; see _complete_indicial_ideal. """ indicial_ideal, _ = self._complete_indicial_ideal(max_multiplier_degree) return indicial_ideal def indicials(self, max_multiplier_degree: int = 4) -> list: """ Returns the indicials (exponents) of this PFIdeal at z = 0, i.e. the roots of its indicial ideal, as a sorted list of (rho, multiplicity) pairs, the multiplicities adding up to the holonomic rank. E.g. the quintic has [((0,), 4)] at its MUM point and [((0,), 1), ((1,), 2), ((2,), 1)] at its conifold point. """ return _indicial_roots(self.indicial_ideal(max_multiplier_degree)) def frobenius_basis(self, order: int, max_log_degree: int = None, max_multiplier_degree: int = 4) -> list: """ Computes a basis of all solutions of this PFIdeal at z = 0, logarithmic ones included, as Periods to the given order. For the basis at another point, move that point to the origin first with change_coordinates. The solutions are found one exponent class at a time: the indicials are grouped by their residue mod Z^n, and each class is solved once, from its componentwise smallest exponent rho_min (which becomes the Periods' indicials), so that solutions starting at rho_min + k with k > 0 - e.g. at the indicials 0, 1, 1, 2 of a conifold point - come out of the same solve. Their degree in log(z) is at most max_log_degree if given, and otherwise at most the class' total multiplicity minus 1 (the logarithms of the monodromies' unipotent parts act on the class' solutions as commuting nilpotent maps, so any product of that many of them vanishes) and at most n = len(hodge_numbers) - 1. Within each class the basis is canonical (see _solve_exponent_class): every solution has its own leading term with coefficient 1, at which the class' other solutions vanish. Coefficients are those of plain log(z_i)^b, without factors of 2*pi*i or factorials. This assumes the origin is a normal-crossing point of the discriminant, where every solution is a sum of terms z^rho * (power series) * (polynomial in log z). Elsewhere fewer solutions than the rank turn up, which is logged as a warning - as is a log structure at a MUM point that differs from hodge_numbers. """ n = self.no_variables indicial_ideal, operator_parts = self._complete_indicial_ideal(max_multiplier_degree) theta_ring = indicial_ideal.ring() roots = _indicial_roots(indicial_ideal) classes = {} for rho, multiplicity in roots: classes.setdefault(tuple(r - floor(r) for r in rho), []).append((rho, multiplicity)) exponent_classes = sorted( ((tuple(min(rho[i] for rho, _ in members) for i in range(n)), members) for members in classes.values()), key=lambda exponent_class: exponent_class[0], ) cy_dimension = len(self.hodge_numbers) - 1 if self.hodge_numbers is not None else None basis = [] for rho_min, members in exponent_classes: class_size = sum(multiplicity for _, multiplicity in members) # A solution starting at the indicial rho only turns up once the solve reaches z^(rho - rho_min). required_order = max(sum(rho[i] - rho_min[i] for i in range(n)) for rho, _ in members) if order < required_order: raise ValueError( "Order %d is too low to reach the indicials %s; it must be at least %d." % (order, [rho for rho, _ in members], required_order) ) if max_log_degree is not None: log_degree = max_log_degree elif cy_dimension is not None: log_degree = min(class_size - 1, cy_dimension) else: log_degree = class_size - 1 solutions = self._solve_exponent_class(operator_parts, theta_ring, rho_min, order, log_degree) logger.info( "Exponent class of %s: %d solution(s) with leading terms %s.", rho_min, len(solutions), [_leading_term(coefficients, rho_min) for coefficients in solutions], ) if len(solutions) != class_size: logger.warning( "Found %d solution(s) in the exponent class of %s, whose indicials have total multiplicity " "%d: z = 0 may not be a normal-crossing point of the discriminant, or the degree in log(z), " "capped at %d, is too low.", len(solutions), rho_min, class_size, log_degree, ) basis.extend( Period(no_variables=n, coefficients=coefficients, order=order, indicials=list(rho_min)) for coefficients in solutions ) if self.hodge_numbers is not None: if len(basis) != sum(self.hodge_numbers): logger.warning( "Found %d solution(s), but the rank sum(hodge_numbers) is %d.", len(basis), sum(self.hodge_numbers) ) # At a point of maximal unipotent monodromy - all indicials 0, with a solution of top log-degree # n - the numbers of solutions of each top log-degree are the Hodge numbers. A single exponent # class at 0 isn't enough: at (0, -1/27) of P_{333111}[12], on the MUM divisor z0 = 0, the # indicials are (0, 0) and (0, 1) and there is a log(z0)^4, but no MUM point. top_log_degrees = [max(sum(log_index) for log_index in period.coefficients) for period in basis] if ( len(roots) == 1 and all(r == 0 for r in roots[0][0]) and max(top_log_degrees, default=-1) == cy_dimension ): log_structure = [top_log_degrees.count(degree) for degree in range(cy_dimension + 1)] if log_structure != self.hodge_numbers: logger.warning( "The numbers of solutions per top log-degree %s at this MUM point differ from " "hodge_numbers %s.", log_structure, self.hodge_numbers, ) logger.info("Found a Frobenius basis of %d solution(s) at order %d.", len(basis), order) return basis def find_power_series_solution(self, indicials: list, order: int) -> list: """ Finds the power series solutions (no logs) to this PFIdeal at given indicial exponents/order, i.e. all solutions of the form z^indicials * (power series in z); frobenius_basis finds all solutions, logarithmic ones included. This is _solve_exponent_class with no room for logarithms, on the operators as given (without completing their indicial ideal). A combination of solutions that could only be continued with a log(z)-term, at a resonance, is not a power series solution and is dropped. """ if len(indicials) != self.no_variables: raise ValueError( "Indicials must have length no_variables=%d, got %d." % (self.no_variables, len(indicials)) ) indicials = [QQ(rho) for rho in indicials] theta_ring, operator_parts = self._operator_theta_parts() solutions = [ Period(no_variables=self.no_variables, coefficients=coefficients, order=order, indicials=indicials) for coefficients in self._solve_exponent_class(operator_parts, theta_ring, indicials, order, 0) ] logger.info( "Found %d power series solution(s) for indicials %s at order %d.", len(solutions), indicials, order ) return solutions def _operator_theta_parts(self) -> tuple: # The theta-ring QQ[theta0, ..., theta{n-1}] and, in it, the theta-parts (see # PFOperator._theta_parts) of every nonzero operator of this PFIdeal. theta_ring = PolynomialRing(QQ, self.no_variables, ['theta%d' % i for i in range(self.no_variables)]) operator_parts = [parts for parts in (op._theta_parts(theta_ring) for op in self.operators) if parts] if not operator_parts: raise ValueError("The PFIdeal contains no nonzero operator.") return theta_ring, operator_parts def _complete_indicial_ideal(self, max_multiplier_degree: int = 4) -> tuple: """ Computes the indicial ideal of this PFIdeal at z = 0, together with the operators needed to see it. The indicial ideal consists of the theta-polynomials Q for which some operator in the (left) ideal has initial part z^e * Q(theta), i.e. this as its part at the lowest z-monomial w.r.t. _z_order_key. Its roots are the indicials (exponents) of the solutions at z = 0, and for a regular holonomic system, as every Picard-Fuchs system is, its colength - the number of roots counted with multiplicity - is the holonomic rank. See M. Saito, B. Sturmfels, N. Takayama, "Gröbner Deformations of Hypergeometric Differential Equations", Algorithms and Computation in Mathematics vol. 6, Springer, 2000, ch. 2. The initial parts of the generators alone can generate a strictly smaller ideal. E.g. for P_{22211}[8] moved to (z0, z1) = (0, 1/4) they are theta0^2*theta1 and theta1*(2*theta1 - 2*theta0 - 1), which vanish on the whole line theta1 = 0. The missing elements are initial parts of combinations of the generators in which the generators' own initial parts cancel (S-pairs), and are found here by linear algebra: the left multiples z^a * theta^c * L of the generators L with |a| + |c| <= d span a finite-dimensional space (a Macaulay matrix), whose echelon form w.r.t. the order on z-monomials exposes every initial part it contains. d grows until the result is zero-dimensional and complete: of colength rank = sum(hodge_numbers) if those are known, and otherwise unchanged from d - 1. Returns (indicial ideal, operator parts): the theta-parts (see PFOperator._theta_parts) of the generators, followed by those of the few extra operators whose initial parts were needed - solving with all of them sees exactly the right starting terms. """ n = self.no_variables theta_ring, generator_parts = self._operator_theta_parts() rank = sum(self.hodge_numbers) if self.hodge_numbers is not None else None generator_initials = [parts[min(parts, key=_z_order_key)] for parts in generator_parts] # A single ordinary differential operator needs no completion: initial parts are multiplicative, so # its indicial ideal is generated by its own initial part (Fuchs' indicial polynomial). single_ode = n == 1 and len(generator_parts) == 1 previous = None for degree in range(max_multiplier_degree + 1): # One row per left multiple, one column per (z-monomial, theta-monomial), the z-monomials in # z-order, so that each row of the echelon form starts at its initial z-monomial. rows = [ _left_multiply(parts, multiplier[:n], multiplier[n:], theta_ring) for parts in generator_parts for multiplier in _multi_indices(2 * n, degree) ] columns = sorted( {(p, tuple(q)) for row in rows for p, P in row.items() for q in P.dict()}, key=lambda column: (_z_order_key(column[0]), column[1]), ) column_position = {column: j for j, column in enumerate(columns)} # Dense on purpose, here and in _solve_exponent_class: Sage's echelon form of *sparse* rational # matrices can be wrong (seen with Sage 9.7 for a 1 x 2 matrix met when solving the quintic at # its orbifold point). macaulay = matrix(QQ, len(rows), len(columns)) for i, row in enumerate(rows): for p, P in row.items(): for q, coeff in P.dict().items(): macaulay[i, column_position[(p, tuple(q))]] = coeff echelon = macaulay.echelon_form() # Keep, greedily, every echelon row whose initial part isn't in the ideal yet, then drop again # those made redundant by rows kept later: every extra operator costs a full pass of the # solver, and very few are needed (one for P_{22211}[8] at (0, 1/4)). extra_parts, extra_initials = [], [] current = theta_ring.ideal(generator_initials) for i, pivot in enumerate(echelon.pivots()): row = echelon.row(i).dict() initial_z_index = columns[pivot][0] initial = sum( coeff * theta_ring.monomial(*columns[j][1]) for j, coeff in row.items() if columns[j][0] == initial_z_index ) if current.reduce(initial) == 0: continue parts = {} for j, coeff in row.items(): p, q = columns[j] parts[p] = parts.get(p, theta_ring(0)) + coeff * theta_ring.monomial(*q) extra_parts.append(parts) extra_initials.append(initial) current = theta_ring.ideal(generator_initials + extra_initials) for i in reversed(range(len(extra_parts))): others = generator_initials + extra_initials[:i] + extra_initials[i + 1:] if theta_ring.ideal(others).reduce(extra_initials[i]) == 0: del extra_parts[i], extra_initials[i] indicial_ideal = theta_ring.ideal(generator_initials + extra_initials) finite = indicial_ideal.dimension() <= 0 colength = indicial_ideal.vector_space_dimension() if finite else None logger.debug( "Indicial ideal at multiplier degree %d with %d extra operator(s): %s.", degree, len(extra_parts), "colength %d" % colength if finite else "not zero-dimensional", ) if finite and rank is not None: if colength == rank: break if colength < rank: raise ValueError( "The indicial ideal at z = 0 has colength %d, less than the rank sum(hodge_numbers) = %d: " "an operator is not in the Picard-Fuchs ideal, or hodge_numbers is wrong." % (colength, rank) ) elif finite and (single_ode or indicial_ideal == previous): if not single_ode: logger.warning( "Taking the indicial ideal at z = 0 as complete since it didn't change from multiplier " "degree %d to %d; give the PFIdeal hodge_numbers to certify this by its rank.", degree - 1, degree, ) break previous = indicial_ideal else: raise ValueError( "Could not complete the indicial ideal at z = 0 up to multiplier degree %d: the operators may not " "generate the whole Picard-Fuchs ideal (whose rank then exceeds sum(hodge_numbers)), or a higher " "max_multiplier_degree is needed." % max_multiplier_degree ) logger.info( "Completed the indicial ideal at z = 0 at multiplier degree %d with %d extra operator(s): colength %d.", degree, len(extra_parts), colength, ) return indicial_ideal, generator_parts + extra_parts def _solve_exponent_class( self, operator_parts: list, theta_ring, indicials, order: int, max_log_degree: int ) -> list: """ Finds all solutions y = sum_k z^(indicials + k) * g_k(log z) of the operators given by their theta-parts (see PFOperator._theta_parts), to the given order: k runs over the multi-indices >= 0, and every g_k is a polynomial in L_i = log(z_i) of total degree at most max_log_degree. Returns one coefficients dictionary {log_index: {z_index: coeff}} (see Period) per solution of a canonical basis. Since theta_i acts on z^(indicials + k) * g(L) as indicials_i + k_i + d/dL_i, a term c*z^p*theta^q sends it to c*z^(indicials + k + p) * (indicials + k + d/dL)^q g(L): it shifts k up by p (p >= 0), never down or sideways. So, as in Frobenius' method, the solutions can be found z-monomial by z-monomial - only with every unknown a polynomial g_k rather than a number, on whose coefficients a theta-polynomial P acts by the matrix of P(indicials + k + d/dL). Logarithms, indicials differing by integers (resonance) and several operators constraining each other are then all handled by the kernels and consistency conditions of these small matrices, with no derivatives with respect to the exponent and no case analysis. The basis is canonical: in the reduced echelon form of all solutions, each has its own leading coefficient 1 at which all others vanish, so in particular no solution's tail contains another one's leading term - as for the canonical series solutions of M. Saito, B. Sturmfels, N. Takayama, "Gröbner Deformations of Hypergeometric Differential Equations", Algorithms and Computation in Mathematics vol. 6, Springer, 2000, ch. 2. """ n = self.no_variables theta_gens = theta_ring.gens() indicials = [QQ(rho) for rho in indicials] # Monomials L^b spanning the polynomials of total degree <= max_log_degree in L, lowest degree first # (1, L0, L1, L0^2, L0*L1, ...); this order decides which coefficient of a solution is normalised. log_indices = sorted(_multi_indices(n, max_log_degree), key=lambda b: (sum(b), tuple(-x for x in b))) log_position = {log_index: i for i, log_index in enumerate(log_indices)} size = len(log_indices) # By Taylor's formula, P(s + d/dL) sends L^b to the sum over b' <= b of # binom(b, b') * (d^(b - b')P/dtheta^(b - b'))(s) * L^b'. With this pattern of entries and every # theta-polynomial's derivatives computed once, each matrix below costs only a few evaluations. action_pattern = [ ( log_position[lower], log_position[upper], tuple(upper[i] - lower[i] for i in range(n)), prod(binomial(upper[i], lower[i]) for i in range(n)), ) for upper in log_indices for lower in log_indices if all(lower[i] <= upper[i] for i in range(n)) ] def taylor_derivatives(P): derivatives = {} for alpha in _multi_indices(n, max_log_degree): derivative = P for i in range(n): if alpha[i]: derivative = derivative.derivative(theta_gens[i], alpha[i]) if derivative != 0: derivatives[alpha] = derivative return derivatives def action_matrix(derivatives, shift): values = {alpha: derivative(*shift) for alpha, derivative in derivatives.items()} entries = {} for row, column, alpha, binomial_factor in action_pattern: if values.get(alpha, 0) != 0: entries[(row, column)] = binomial_factor * values[alpha] # matrix() would make a sparse matrix from a dictionary; see the comment on dense matrices in # _complete_indicial_ideal. return matrix(QQ, size, size, entries, sparse=False) operator_derivatives = [{p: taylor_derivatives(P) for p, P in parts.items()} for parts in operator_parts] supports = [sorted(parts, key=_z_order_key) for parts in operator_parts] z_indices = sorted(_multi_indices(n, order), key=_z_order_key) # solved[k] holds g_k as a size x dimension matrix, whose column j holds g_k's coefficients # (w.r.t. log_indices) in the j-th of the `dimension` independent solutions found so far. dimension = 0 solved = {} for k in z_indices: # Collect the equations whose latest unknown is g_k. The coefficient of z^(indicials + m) in L(y) # involves g_(m - p) for all p <= m (componentwise) in L's support, the latest of them in z-order # being g_(m - p') for the smallest such p'. Taking m = k + p' for those p' that are the smallest # element p <= k + p' of the support thus visits every equation exactly once - also for operators # whose lowest z-monomial isn't 1, like the extra ones from _complete_indicial_ideal. diagonal_blocks, known_blocks = [], [] for derivatives, support in zip(operator_derivatives, supports): for position, pivot_exponent in enumerate(support): m = tuple(k[i] + pivot_exponent[i] for i in range(n)) if any(all(p[i] <= m[i] for i in range(n)) for p in support[:position]): continue diagonal_blocks.append( action_matrix(derivatives[pivot_exponent], [indicials[i] + k[i] for i in range(n)]) ) known = matrix(QQ, size, dimension) if dimension: for p in support[position + 1:]: if all(p[i] <= m[i] for i in range(n)): j = tuple(m[i] - p[i] for i in range(n)) # earlier than k in z-order: solved known += action_matrix(derivatives[p], [indicials[i] + j[i] for i in range(n)]) * solved[j] known_blocks.append(known) diagonal = block_matrix(QQ, [[block] for block in diagonal_blocks], subdivide=False) known = block_matrix(QQ, [[block] for block in known_blocks], subdivide=False) # Row-reduce [diagonal | known]: rows with a pivot in the diagonal part give g_k in terms of the # solutions so far, rows without are linear conditions on those solutions, and the diagonal # part's free columns start new solutions at z^(indicials + k), where the indicial ideal vanishes. reduced = diagonal.augment(known).rref() conditions = [row[size:] for row in reduced.rows() if row[:size] == 0 and row[size:] != 0] if conditions: # Some combinations of the solutions so far can't be continued past this point - e.g. as they # would need a higher power of log(z) than max_log_degree. Keep those that can, by rewriting # every g_j in a basis of the combinations satisfying all conditions. restriction = matrix(QQ, conditions).right_kernel_matrix().transpose() for j in solved: solved[j] = solved[j] * restriction dimension = restriction.ncols() reduced = diagonal.augment(known * restriction).rref() pivots = [column for column in reduced.pivots() if column < size] free = [column for column in range(size) if column not in pivots] # The free coefficients are 0 in the continuations of the solutions so far, so that no tail # contains the leading term of a solution starting here; each new solution has one of them 1. particular = matrix(QQ, size, dimension) new_solutions = matrix(QQ, size, len(free)) for row_index, column in enumerate(pivots): particular.set_row(column, -reduced.row(row_index)[size:]) for j, free_column in enumerate(free): new_solutions[column, j] = -reduced[row_index, free_column] for j, free_column in enumerate(free): new_solutions[free_column, j] = 1 if free: for j in solved: solved[j] = solved[j].augment(matrix(QQ, size, len(free))) solved[k] = particular.augment(new_solutions) dimension += len(free) if dimension == 0: return [] # Canonical basis: the reduced echelon form of all solutions' coefficients, ordered by z-monomial (in # z-order) and then log-monomial. Each solution then has its own leading coefficient 1 at which all # others vanish, independently of the order in which the solutions turned up above. coefficient_labels = [(k, log_index) for k in z_indices for log_index in log_indices] echelon = block_matrix(QQ, [[solved[k]] for k in z_indices], subdivide=False).transpose().rref() solutions = [] for solution_vector in echelon.rows()[:dimension]: by_log_index = {} for (k, log_index), coeff in zip(coefficient_labels, solution_vector): by_log_index.setdefault(log_index, {})[k] = coeff # Every z-index up to order is kept, zeros included, for each log-monomial that occurs. solutions.append({ log_index: z_dict for log_index, z_dict in by_log_index.items() if any(coeff != 0 for coeff in z_dict.values()) }) return solutions def _z_order_key(z_index) -> tuple: # The order on z-monomials used throughout: by total degree, ties broken lexicographically. Being graded, # it is compatible with truncating at a total degree; being a monomial order (p < q implies p + r < q + r), # it acts like a generic weight vector, so an operator's initial (lowest) part is a single z-monomial # times a theta-polynomial, and solving in this order only ever needs earlier unknowns. return (sum(z_index),) + tuple(z_index) def _left_multiply(parts: dict, z_shift, theta_power, theta_ring) -> dict: # Normal-ordered product (z^a * theta^c) * (sum_p z^p * P_p(theta)) = sum_p z^(a + p) * (theta + p)^c * P_p(theta) # of theta-parts (see PFOperator._theta_parts), by theta_i * z^p = z^p * (theta_i + p_i). theta_gens = theta_ring.gens() n = len(theta_gens) product = {} for p, P in parts.items(): shifted_power = prod((theta_gens[i] + p[i]) ** theta_power[i] for i in range(n)) product[tuple(z_shift[i] + p[i] for i in range(n))] = shifted_power * P return product def _indicial_roots(indicial_ideal) -> list: # The roots of a zero-dimensional indicial ideal with their multiplicities, as a sorted list of # (rho, multiplicity): each primary component belongs to one root, and its colength is that multiplicity. theta_gens = indicial_ideal.ring().gens() roots = [] for primary, prime in indicial_ideal.complete_primary_decomposition(): points = prime.variety() if len(points) != 1 or prime.vector_space_dimension() != 1: raise NotImplementedError( "The indicial ideal has a non-rational root (prime component %s), but Picard-Fuchs systems of " "geometric families only have rational indicials (quasi-unipotent monodromy)." % (prime,) ) roots.append((tuple(points[0][theta] for theta in theta_gens), primary.vector_space_dimension())) return sorted(roots) def _leading_term(coefficients: dict, indicials) -> str: # The terms of a solution (a coefficients dictionary, see Period) at its z-order-smallest z-monomial, # e.g. 'sqrt(z1)*(log(z0) + log(z1))', to log the structure of a basis. n = len(indicials) z_gens = [SR.var('z%d' % i) for i in range(n)] leading = min( (z_index for z_dict in coefficients.values() for z_index, coeff in z_dict.items() if coeff != 0), key=_z_order_key, ) log_polynomial = sum( z_dict[leading] * prod(log(z_gens[i]) ** log_index[i] for i in range(n)) for log_index, z_dict in coefficients.items() ) return str(prod(z_gens[i] ** (indicials[i] + leading[i]) for i in range(n)) * log_polynomial)