From 51a5612fe12a58e5b7263e7f2a371420ad64f99a Mon Sep 17 00:00:00 2001 From: julian Date: Mon, 17 Aug 2026 19:51:31 +0200 Subject: [PATCH] Initial setup for period computation (#13) Introduces: - class for Picard--Fuchs operators and their ideals - class for periods (their complex linear combinations in the Frobenius bases) Implements: - method to obtain operator from period - method to get power series solution (at given indicials) to PF ideal --------- Co-authored-by: Julian Piribauer Reviewed-on: https://gitea.piribauer.ch/julian/Calabi-Yau-Period-Geometry/pulls/13 --- README.md | 89 +++++- pyproject.toml | 3 + sage/period_computation.sage | 503 +++++++++++++++++++++++++++++++ sage/util.py | 15 + tests/test_period_computation.py | 194 ++++++++++++ 5 files changed, 798 insertions(+), 6 deletions(-) create mode 100644 sage/period_computation.sage create mode 100644 sage/util.py create mode 100644 tests/test_period_computation.py diff --git a/README.md b/README.md index 8ae320a..1b5b95a 100644 --- a/README.md +++ b/README.md @@ -12,12 +12,21 @@ elliptic_curve_D = ToricPolytopeProjectiveSpace([1, 2, 3], model_name="elliptic_ CY3_quintic = ToricPolytopeProjectiveSpace([1, 1, 1, 1, 1], model_name="quintic") CY3_bicubic = ToricPolytopeCICY([[3, 3]], model_name="bi-cubic") -CICY3_two_parameter_manual_nef = ToricPolytope([[1,0,0,0,0,0],[0,1,0,0,0,0],[0,0,1,0,0,0], - [0,0,0,1,0,0],[0,0,0,0,1,0],[0,0,0,0,0,1], - [-1,-1,0,0,0,0],[0,0,-1,-1,-1,-1]], - nef_partition=[[0,1,2,3], [4,5,6,7]]) +CICY3_two_parameter_manual_nef = ToricPolytope( + [ + [1, 0, 0, 0, 0, 0], + [0, 1, 0, 0, 0, 0], + [0, 0, 1, 0, 0, 0], + [0, 0, 0, 1, 0, 0], + [0, 0, 0, 0, 1, 0], + [0, 0, 0, 0, 0, 1], + [-1, -1, 0, 0, 0, 0], + [0, 0, -1, -1, -1, -1], + ], + nef_partition=[[0, 1, 2, 3], [4, 5, 6, 7]], +) -CICY5_two_parameter = ToricPolytopeCICY([[6, 1],[0, 2]]) +CICY5_two_parameter = ToricPolytopeCICY([[6, 1], [0, 2]]) ``` The discriminant factors and topological data, e.g. for the quintic, can then be computed with the methods below. @@ -71,4 +80,72 @@ DEBUG:__main__: The result is saved in the folder `data/topdata` as a JSON file — giving a model name helps keeping track of these outputs. -Note that for Calabi–Yau dimensions larger than four, the additional \ No newline at end of file +Note that for Calabi–Yau dimensions larger than four, the additional + +## period_computation + +`period_computation.sage` provides classes for working with Picard–Fuchs operators and their +period solutions: `PFOperator`, `PFIdeal` and `Period`, together with Ansatz variants of the first and +last (`PFOperatorAnsatz`, `PeriodAnsatz`) used to search for unknown operators or periods of a given +z- and theta-degree. + +A `PFOperator` is parsed from a string in the variables `z0, ..., z` and `theta0, ..., theta`, +the logarithmic derivatives theta_i = z_i d/dz_i. For example, the quintic's Picard–Fuchs operator: + +```python +L = PFOperator( + "theta0^4 - 3125*z0*theta0^4 - 6250*z0*theta0^3 - 4375*z0*theta0^2 - 1250*z0*theta0 - 120*z0", + no_variables=1, +) + +L.simplify().operator_string +``` + +```term +'-5*(5*theta0 + 4)*(5*theta0 + 3)*(5*theta0 + 2)*(5*theta0 + 1)*z0 + theta0^4' +``` + +An operator (or a `PFIdeal` of several) can be solved for its power series solution at given indicial +exponents and order. + +```python +ideal = PFIdeal([L]) +period = ideal.find_power_series_solution(indicials=[0], order=3)[0] + +period.period_string +``` + +```term +'168168000*z0^3 + 113400*z0^2 + 120*z0 + 1' +``` + +The reverse direction is supported too: given a `Period`, `find_annihilating_operators` searches for +`PFOperator`s of a given z- and theta-degree that annihilate it, by solving an Ansatz of unknown +coefficients via linear algebra. Both directions extend to several moduli, e.g. for the two-parameter +model P_{2,2,2,1,1}[8]: + +```python +M1 = PFOperator( + "theta1*(-2*theta0 + 2*theta1 - 1) + 2*(theta0 - 2*theta1 - 1)*(theta0 - 2*theta1)*z1", + no_variables=2, +) +M2 = PFOperator( + "theta0^2*(2*(theta0 - 2*theta1)*z1 - theta1) - 16*(2*theta0 + 1)*(4*theta0 + 1)*(4*theta0 + 3)*z0*z1", + no_variables=2, +) +ideal = PFIdeal([M1, M2]) + +period = ideal.find_power_series_solution(indicials=[0, 1 / 2], order=6)[0] +recovered = period.find_annihilating_operators(z_degree=1, theta_degree=2) + +period.period_string +recovered[0].simplify().operator_string +``` + +```term +'-1/45045*(60886425600*z0^3*z1^3 - 2767564800*z0^2*z1^4 + 100638720*z0*z1^5 - 14192640*z1^6 + 830269440*z0^2*z1^3 - 30750720*z0*z1^4 + 4193280*z1^5 - 242161920*z0^2*z1^2 + 9884160*z0*z1^3 - 1281280*z1^4 - 3459456*z0*z1^2 + 411840*z1^3 + 1441440*z0*z1 - 144144*z1^2 + 60060*z1 - 45045)*sqrt(z1)' +'-2*(theta0 - 2*theta1)*(theta0 - 2*theta1 - 1)*z1 + (2*theta0 - 2*theta1 + 1)*theta1' +``` + +`recovered[0]` is, up to scale, `M1` — recovered purely from `M1`, `M2`'s shared power series +solution. \ No newline at end of file diff --git a/pyproject.toml b/pyproject.toml index 4b852c2..61e6732 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -15,3 +15,6 @@ ignore = [ # test_smoke.py star-imports sage.all and loads a .sage file, so ruff can't see where its names come from. "tests/test_topdata_and_disc.py" = ["F403", "F405"] + +# Same star-import + load() pattern as above. +"tests/test_period_computation.py" = ["F403", "F405"] diff --git a/sage/period_computation.sage b/sage/period_computation.sage new file mode 100644 index 0000000..70d5a65 --- /dev/null +++ b/sage/period_computation.sage @@ -0,0 +1,503 @@ +import copy +import logging +from sage.all import sage_eval, PolynomialRing, QQ, SR, log, matrix, prod, 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 + +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. + """ + def __init__(self, operators: list): + 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.") + self.no_variables = operators[0].no_variables + self.operators = operators + 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 find_power_series_solution(self, indicials: list, order: int) -> list: + """ + Finds power series solutions (no logs) to this PFIdeal at given indicial exponents/order. + + 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. + """ + 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] + + 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) + ] + + logger.info( + "Found %d power series solution(s) for indicials %s at order %d.", len(solutions), indicials, order + ) + return solutions + + +class Period: + """ + A class representing a period as a formal power series in z-variables and their logs. + + The coefficients are stored in a dictionary of dictionaries, where the first key is the multi-index + of the logarithmic part and the second key is the multi-index of the z-variables. + So, for example, the coefficient of + + log(z0)^2 * log(z1) * z0^3 * z1^2 would be stored as coefficients[(2, 1)][(3, 2)] + + """ + + def _initialise_ring(self): + z_names = ['z%d' % i for i in range(self.no_variables)] + log_names = ['L%d' % i for i in range(self.no_variables)] + ring = PolynomialRing(QQ, z_names + log_names) + z_gens = ring.gens()[:self.no_variables] + log_gens = ring.gens()[self.no_variables:] + + return ring, z_gens, log_gens # type: (PolynomialRing, list, list) + + def _period_from_string(self, period_string: str) -> dict: + ring, z_gens, log_gens = self._initialise_ring() + log_of_z = dict(zip(z_gens, log_gens)) + + def log(zi): + try: + return log_of_z[zi] + except (KeyError, TypeError): + raise ValueError( + "log(...) may only be applied to one of the z-variables z0, ..., z%d" + % (self.no_variables - 1) + ) + + parse_locals = dict(ring.gens_dict()) + parse_locals['log'] = log + + try: + expression = ring(sage_eval(period_string, locals=parse_locals)) + except Exception as e: + raise ValueError("Invalid period string: %s" % e) + + coefficients = {} + for coeff, monomial in expression: + exponents = monomial.exponents()[0] + z_index = tuple(exponents[:self.no_variables]) + log_index = tuple(exponents[self.no_variables:]) + coefficients.setdefault(log_index, {})[z_index] = coeff + + return coefficients + + def _period_to_string(self) -> str: + ring, z_gens, log_gens = self._initialise_ring() + expression = ring(0) + for log_index, z_dict in self.coefficients.items(): + for z_index, coeff in z_dict.items(): + monomial = coeff + for i in range(self.no_variables): + monomial *= (log_gens[i] ** log_index[i]) * (z_gens[i] ** z_index[i]) + expression += monomial + + log_substitutions = {log_gens[i]: log(SR(z_gens[i])) for i in range(self.no_variables)} + result = SR(expression).subs(log_substitutions) + + # A nonzero indicial ρ_i means the coefficients above are for z_i^k, but the + # actual solution is z_i^(ρ_i + k); make that explicit in the printed form. + indicial_prefactor = prod( + (SR(z_gens[i]) ** self.indicials[i] for i in range(self.no_variables) if self.indicials[i] != 0), + SR(1), + ) + if indicial_prefactor != 1: + result *= indicial_prefactor + + return str(result) + + def _max_z_degree(self) -> int: + # Highest total z-degree (sum of the z-multi-index) among all coefficients. + if not self.coefficients: + return 0 + return max( + sum(z_index) + for z_dict in self.coefficients.values() + for z_index in z_dict + ) + + def _truncate_coefficients(self, order: int) -> dict: + # Drop every (log_index, z_index) entry whose total z-degree exceeds order, + # removing it from the dictionary rather than merely zeroing it out. + truncated = {} + for log_index, z_dict in self.coefficients.items(): + kept = {z_index: coeff for z_index, coeff in z_dict.items() if sum(z_index) <= order} + if kept: + truncated[log_index] = kept + return truncated + + def _apply_theta(self, coefficients: dict, index: int) -> dict: + """ + Applies the logarithmic derivative theta_i = z_i * d/dz_i once to a coefficients dict of the same + shape as self.coefficients. It uses the product rule + + theta_i(z^a log(z)^k) = a_i * z^a log(z)^k + k_i * z^a log(z)^(k - e_i). + + z_index entries are offsets from self.indicials: a stored z_index of a really means + z^(a + self.indicials[index]), so theta_i's eigenvalue is a_i + self.indicials[index] + rather than the bare a_i (self.indicials is all-zero unless explicitly given, in which + case this reduces to the ordinary power-series rule). + """ + result = {} + for log_index, z_dict in coefficients.items(): + for z_index, coeff in z_dict.items(): + a_i = z_index[index] + self.indicials[index] + if a_i != 0: + inner = result.setdefault(log_index, {}) + inner[z_index] = inner.get(z_index, 0) + a_i * coeff + + k_i = log_index[index] + if k_i != 0: + lowered_log_index = log_index[:index] + (k_i - 1,) + log_index[index + 1:] + inner = result.setdefault(lowered_log_index, {}) + inner[z_index] = inner.get(z_index, 0) + k_i * coeff + return result + + def apply_operator(self, pf_operator: PFOperator) -> "Period": + """ + Applies a PFOperator to this period and returns the result as a new Period, + truncated to self.order. Each monomial of the operator is normalised as + coeff * z^p * theta^q (see PFOperator's docstring), so theta^q is applied + to the period first and the result is then multiplied by coeff * z^p. + """ + if pf_operator.no_variables != self.no_variables: + raise ValueError( + "Variable count mismatch: period has %d variable(s), operator has %d." + % (self.no_variables, pf_operator.no_variables) + ) + + result_coefficients = {} + for monomial_coeff, monomial in pf_operator.operator: + exponents = monomial.exponents()[0] + z_exponents = exponents[:self.no_variables] + theta_exponents = exponents[self.no_variables:] + + term = self.coefficients + for i, power in enumerate(theta_exponents): + for _ in range(power): + term = self._apply_theta(term, i) + + for log_index, z_dict in term.items(): + inner = result_coefficients.setdefault(log_index, {}) + for z_index, coeff in z_dict.items(): + shifted_z_index = tuple(z_index[i] + z_exponents[i] for i in range(self.no_variables)) + inner[shifted_z_index] = inner.get(shifted_z_index, 0) + monomial_coeff * coeff + + result_coefficients = { + log_index: {z_index: c for z_index, c in z_dict.items() if c != 0} + for log_index, z_dict in result_coefficients.items() + } + result_coefficients = {log_index: z_dict for log_index, z_dict in result_coefficients.items() if z_dict} + + return Period( + no_variables=self.no_variables, + coefficients=result_coefficients, + order=self.order, + indicials=self.indicials, + ) + + def find_annihilating_operators(self, z_degree: int, theta_degree: int) -> list: + """ + Finds operators annihilating this period among PFOperator Ansätze of the given + z- and theta-degree. + """ + ansatz = PFOperatorAnsatz(self.no_variables, theta_degree, z_degree) + unknowns = ansatz.unknowns + + equations = [ + coeff + for z_dict in self.apply_operator(ansatz).coefficients.values() + for coeff in z_dict.values() + ] + + coefficient_rows = [[SR(equation).coefficient(b) for b in unknowns] for equation in equations] + kernel_basis = matrix(QQ, coefficient_rows, ncols=len(unknowns)).right_kernel().basis() + + operators = [] + for basis_vector in kernel_basis: + solution = {unknowns[j]: basis_vector[j] for j in range(len(unknowns))} + solved_operator = ansatz.operator.map_coefficients(lambda c: SR(c).subs(solution)) + operators.append(PFOperator(str(solved_operator), no_variables=self.no_variables)) + + return operators + + def __init__( + self, + no_variables: int = 1, + coefficients: dict = None, + period_string: str = None, + order: int = None, + indicials: list = None, + ): + self.no_variables = no_variables + # indicials[i] is the (rational) Frobenius exponent ρ_i of z_i: a stored + # z-index of a really represents z^(a + indicials[i]). + self.indicials = list(indicials) if indicials is not None else [0] * no_variables + + # Use coefficients or period_string to initialize the period + if period_string is not None: + if coefficients is not None: + raise ValueError("Provide either coefficients or period_string, not both.") + self.coefficients = self._period_from_string(period_string) + self.period_string = period_string + logger.debug("Using string for initialisation.") + elif coefficients is None: + self.coefficients = {} + else: + self.coefficients = coefficients + self.period_string = self._period_to_string() + logger.debug("Using coefficients for initialisation.") + + if order is None: + self.order = self._max_z_degree() + logger.debug("Order not provided, using maximum z-degree: %d", self.order) + else: + self.coefficients = self._truncate_coefficients(order) + self.order = order + self.period_string = self._period_to_string() + logger.debug("Truncated coefficients to order %d.", self.order) + + logger.info("Initialised Period in %d variables at order %d.", self.no_variables, self.order) + + +class PeriodAnsatz(Period): + """ + A class for period Ansätze, characterised by number of variables and their z- and log-multi-degrees. + """ + + def __init__(self, no_variables: int, z_degree: int, log_degree: int, indicials: list = None): + self.no_variables = no_variables + self.z_degree = z_degree + self.log_degree = log_degree + + log_indices = _multi_indices(self.no_variables, log_degree) + z_indices = _multi_indices(self.no_variables, z_degree) + + # Every (log_index, z_index) pair allowed by the degree bounds gets its own + # fresh unknown, to be solved for once a PFOperator is applied to the Ansatz. + coefficients = { + log_index: { + z_index: var("a_" + "_".join(str(x) for x in log_index + z_index)) + for z_index in z_indices + } + for log_index in log_indices + } + + super().__init__(no_variables=no_variables, coefficients=coefficients, order=z_degree, indicials=indicials) + self.expansion_coefficients = self.coefficients + logger.info( + "Initialised PeriodAnsatz with %d unknown coefficient(s).", + sum(len(z_dict) for z_dict in self.expansion_coefficients.values()), + ) \ No newline at end of file diff --git a/sage/util.py b/sage/util.py new file mode 100644 index 0000000..2301b1c --- /dev/null +++ b/sage/util.py @@ -0,0 +1,15 @@ +""" +Utility functions for the Calabi-Yau period geometry project. +""" + +def _multi_indices(no_variables, total_degree: int) -> list: + # All no_variables-tuples of non-negative integers whose entries sum to at most total_degree. + def helper(remaining_vars, remaining_degree): + if remaining_vars == 0: + yield () + return + for first in range(remaining_degree + 1): + for rest in helper(remaining_vars - 1, remaining_degree - first): + yield (first,) + rest + + return list(helper(no_variables, total_degree)) diff --git a/tests/test_period_computation.py b/tests/test_period_computation.py new file mode 100644 index 0000000..81a2270 --- /dev/null +++ b/tests/test_period_computation.py @@ -0,0 +1,194 @@ +import os + +from sage.all import * # noqa: F401 + +load(os.path.join(os.path.dirname(__file__), "..", "sage", "period_computation.sage")) + + +def test_apply_operator_univariate(): + # theta^2 - z applied to log(z) should give -z*log(z), since theta(log z) = 1 + # and theta^2(log z) = theta(1) = 0. + period = Period(no_variables=1, coefficients={(1,): {(0,): 1}}, order=5) + op = PFOperator("theta0^2 - z0", no_variables=1) + + result = period.apply_operator(op) + + assert result.coefficients == {(1,): {(1,): -1}} + + +def test_apply_operator_mixes_variables(): + # z0*theta1 applied to z0*log(z1): theta1 strips log(z1) down to a bare 1 + # (leaving z0 untouched), then multiplying by z0 gives z0^2. + period = Period(no_variables=2, coefficients={(0, 1): {(1, 0): 1}}, order=5) + op = PFOperator("z0*theta1", no_variables=2) + + result = period.apply_operator(op) + + assert result.coefficients == {(0, 0): {(2, 0): 1}} + + +def test_apply_operator_truncates_to_order(): + # Multiplying by z0^2 pushes some terms above the period's order, so they + # should be dropped rather than kept with a nonzero coefficient. + period = Period(no_variables=1, coefficients={(0,): {(0,): 1, (1,): 1, (2,): 1}}, order=2) + op = PFOperator("z0^2", no_variables=1) + + result = period.apply_operator(op) + + assert result.coefficients == {(0,): {(2,): 1}} + assert result.order == 2 + + +def test_find_annihilating_operators_recovers_theta_squared(): + # theta0^2 annihilates log(z0): theta0(log z0) = 1, theta0^2(log z0) = theta0(1) = 0. + # Among degree-(0, 2) Ansatze this should be the only solution, up to scaling. + period = Period(no_variables=1, coefficients={(1,): {(0,): 1}}, order=5) + + ops = period.find_annihilating_operators(z_degree=0, theta_degree=2) + + # A one-dimensional solution space means exactly one basis operator. + assert len(ops) == 1 + assert period.apply_operator(ops[0]).coefficients == {} + assert str(ops[0].operator) == "theta0^2" + + +def test_find_annihilating_operators_recovers_geometric_series_operator(): + # sum_{k=0}^{4} z0^k is annihilated (up to truncation order) by + # (1 - z0)*theta0 - z0, i.e. theta0 - z0*theta0 - z0. + period = Period( + no_variables=1, + coefficients={(0,): {(0,): 1, (1,): 1, (2,): 1, (3,): 1, (4,): 1}}, + order=4, + ) + + ops = period.find_annihilating_operators(z_degree=1, theta_degree=1) + + assert len(ops) == 1 + assert period.apply_operator(ops[0]).coefficients == {} + + +def test_find_annihilating_operators_returns_empty_list_when_no_solution_exists(): + # No degree-0 (constant) operator other than the zero operator can annihilate a + # nonzero constant period, so the linear system's only solution is trivial. + period = Period(no_variables=1, coefficients={(0,): {(0,): 1}}, order=0) + + ops = period.find_annihilating_operators(z_degree=0, theta_degree=0) + + assert ops == [] + + +def test_simplify_factorises_theta_polynomial_per_z_monomial(): + # theta0^4 - 5*z0*(5*theta0+1)*(5*theta0+2)*(5*theta0+3)*(5*theta0+4), expanded, is the + # quintic's Picard-Fuchs operator. simplify() should recover the factorised form: the + # z0^0 part (theta0^4) has no theta-factor to pull out, while the z0^1 part factorises + # into the four linear pieces. + expanded = "theta0^4 - 3125*z0*theta0^4 - 6250*z0*theta0^3 - 4375*z0*theta0^2 - 1250*z0*theta0 - 120*z0" + op = PFOperator(expanded, no_variables=1) + + simplified = op.simplify() + + # The underlying (expanded) operator is unchanged - only the display string differs. + assert simplified.operator == op.operator + assert simplified.operator_string == "-5*(5*theta0 + 4)*(5*theta0 + 3)*(5*theta0 + 2)*(5*theta0 + 1)*z0 + theta0^4" + + +def test_find_power_series_solution_recovers_quintic_period(): + # The quintic's Picard-Fuchs operator has indicial equation theta0^4 = 0 at z0 = 0 + # (a quadruple root at 0), so its holomorphic (non-logarithmic) power series solution + # is found at indicial 0. Up to normalisation this is the classic quintic period + # 1 + 120*z0 + 113400*z0^2 + ... + op = PFOperator( + "theta0^4 - 5*z0*(5*theta0 + 1)*(5*theta0 + 2)*(5*theta0 + 3)*(5*theta0 + 4)", + no_variables=1, + ) + ideal = PFIdeal([op]) + + solutions = ideal.find_power_series_solution(indicials=[0], order=3) + + assert len(solutions) == 1 + assert solutions[0].coefficients == {(0,): {(0,): 1, (1,): 120, (2,): 113400, (3,): 168168000}} + + +def test_find_power_series_solution_handles_rational_indicial(): + # theta0 - 1/2 kills z0^(1/2 + k) only when k = 0, since its eigenvalue is 1/2 + k; + # so at indicial 1/2 there is exactly one solution (a constant multiple of sqrt(z0)), + # while at indicial 0 no power series solution exists at all. + op = PFOperator("theta0 - 1/2", no_variables=1) + ideal = PFIdeal([op]) + + solutions = ideal.find_power_series_solution(indicials=[1/2], order=3) + assert len(solutions) == 1 + assert solutions[0].coefficients == {(0,): {(0,): 1, (1,): 0, (2,): 0, (3,): 0}} + assert solutions[0].period_string == "sqrt(z0)" + + assert ideal.find_power_series_solution(indicials=[0], order=3) == [] + + +def test_ideal_finds_power_series_solution_and_recovers_first_operator(): + # Picard--Fuchs ideal for P_{22211}[8]: + # M1 = theta2*(-2*theta1+2*theta2-1) + 2*(theta1-2*theta2-1)*(theta1-2*theta2)*z2 + # M2 = theta1^2*(2*(theta1-2*theta2)*z2-theta2) - 16*(2*theta1+1)*(4*theta1+1)*(4*theta1+3)*z1*z2 + M1 = PFOperator( + "theta1*(-2*theta0 + 2*theta1 - 1) + 2*(theta0 - 2*theta1 - 1)*(theta0 - 2*theta1)*z1", + no_variables=2, + ) + M2 = PFOperator( + "theta0^2*(2*(theta0 - 2*theta1)*z1 - theta1) - 16*(2*theta0 + 1)*(4*theta0 + 1)*(4*theta0 + 3)*z0*z1", + no_variables=2, + ) + ideal = PFIdeal([M1, M2]) + + solutions = ideal.find_power_series_solution(indicials=[0, QQ(1) / 2], order=6) + assert len(solutions) == 1 + + period = solutions[0] + assert period.apply_operator(M1).coefficients == {} + assert period.apply_operator(M2).coefficients == {} + # Normalisation convention: the free parameter at the leading (0, 0) coefficient is 1. + assert period.coefficients[(0, 0)][(0, 0)] == 1 + assert period.coefficients[(0, 0)][(1, 1)] == -32 + + # M1 lives entirely within z_degree <= 1, theta_degree <= 2 (a single bare factor of + # z2, quadratic in theta), so this is the natural Ansatz level to look for it at. + recovered = period.find_annihilating_operators(z_degree=1, theta_degree=2) + assert len(recovered) == 1 + + # recovered[0] should be a scalar multiple of M1 - compare via the coefficient of the + # bare theta2 (theta1 in code) monomial, which is nonzero in M1. + ratio = ( + recovered[0].operator.monomial_coefficient(M1.theta_gens[1]) + / M1.operator.monomial_coefficient(M1.theta_gens[1]) + ) + assert ratio != 0 + assert recovered[0].operator == ratio * M1.operator + + +def test_single_operator_recovers_itself_from_its_holomorphic_period(): + # Ensuring methods work for operator 4.2.1: + L = PFOperator( + "theta0^4 - 4*z0*(2*theta0 + 1)^2*(7*theta0^2 + 7*theta0 + 2) " + "- 128*z0^2*(2*theta0 + 1)^2*(2*theta0 + 3)^2", + no_variables=1, + ) + ideal = PFIdeal([L]) + + solutions = ideal.find_power_series_solution(indicials=[0], order=15) + assert len(solutions) == 1 + + period = solutions[0] + assert period.apply_operator(L).coefficients == {} + assert period.coefficients[(0,)][(0,)] == 1 + assert period.coefficients[(0,)][(1,)] == 8 + assert period.coefficients[(0,)][(2,)] == 360 + + # L lives entirely within z_degree <= 2, theta_degree <= 4, its stated level. + recovered = period.find_annihilating_operators(z_degree=2, theta_degree=4) + assert len(recovered) == 1 + + # recovered[0] should be a scalar multiple of L - compare via the coefficient of the + # bare theta0^4 monomial, which is 1 in L. + ratio = recovered[0].operator.monomial_coefficient(L.theta_gens[0] ** 4) / L.operator.monomial_coefficient( + L.theta_gens[0] ** 4 + ) + assert ratio != 0 + assert recovered[0].operator == ratio * L.operator