diff --git a/src/doc/en/reference/functions/index.rst b/src/doc/en/reference/functions/index.rst index 3916eb0c260..2a0991472e3 100644 --- a/src/doc/en/reference/functions/index.rst +++ b/src/doc/en/reference/functions/index.rst @@ -21,6 +21,7 @@ Built-in Functions sage/functions/other sage/functions/special sage/functions/hypergeometric + sage/functions/hypergeometric_algebraic sage/functions/jacobi sage/functions/airy sage/functions/bessel diff --git a/src/doc/en/reference/references/index.rst b/src/doc/en/reference/references/index.rst index 1f1cd007290..f0a9a602081 100644 --- a/src/doc/en/reference/references/index.rst +++ b/src/doc/en/reference/references/index.rst @@ -1432,6 +1432,9 @@ REFERENCES: An Algorithmic Approach, Algorithms and Computation in Mathematics, Volume 20, Springer (2007) +.. [But2010] Peter Butkovič, *Max-linear systems. Theory and algorithms.* + Springer Monographs in Mathematics. London: Springer. xvii, 272 p. (2010). + .. [Buell89] Duncan A. Buell. *Binary Quadratic Forms: Classical Theory and Modern Computations.* Springer, 1989. @@ -1650,6 +1653,10 @@ REFERENCES: IV. The quotient groups of the lower central series, Ann. of Math. 68 (1958) 81--95. +.. [CFV2025] Xavier Caruso, Florian Fürnsinn, Daniel Vargas-Montoya, + *Galois groups of reductions modulo p of D-finite series*, + :arxiv:`2504.09429`, 2025. + .. [CFZ2000] \J. Cassaigne, S. Ferenczi, L.Q. Zamboni, *Imbalances in Arnoux-Rauzy sequences*, Ann. Inst. Fourier (Grenoble) 50 (2000) 1265--1276. @@ -1724,6 +1731,9 @@ REFERENCES: .. [ChenDB] Eric Chen, Online database of two-weight codes, http://moodle.tec.hkr.se/~chen/research/2-weight-codes/search.php +.. [Chr1986] Gilles Christol, *Fonctions hypergéométriques bornées*. + Groupe de travail d’analyse ultramétrique 14 (1986-87), exp. no 8, p. 1-16 + .. [CHK2001] Keith D. Cooper, Timothy J. Harvey and Ken Kennedy. *A Simple, Fast Dominance Algorithm*, Software practice and Experience, 4:1-10 (2001). @@ -2881,6 +2891,11 @@ REFERENCES: toric varieties defined by atomic lattices*. Inventiones Mathematicae. **155** (2004), no. 3, pp. 515-536. +.. [FY2024] Florian Fürnsinn, Sergey Yurkevich. + *Algebraicity of hypergeometric functions with arbitrary parameters*, + Bulletin of the London Mathematical Society. **56** (2024), + pp. 2824-2846. :doi:`10.1112/blms.13103`, :arxiv:`2308.12855` (2023). + .. [FZ2001] \S. Fomin and A. Zelevinsky. *Cluster algebras I. Foundations*, \J. Amer. Math. Soc. **15** (2002), no. 2, pp. 497-529. :arxiv:`math/0104151` (2001). diff --git a/src/doc/en/reference/semirings/index.rst b/src/doc/en/reference/semirings/index.rst index b40e71c54e1..3bbb10edd0c 100644 --- a/src/doc/en/reference/semirings/index.rst +++ b/src/doc/en/reference/semirings/index.rst @@ -6,5 +6,9 @@ Standard Semirings sage/rings/semirings/non_negative_integer_semiring sage/rings/semirings/tropical_semiring + sage/rings/semirings/tropical_polynomial + sage/rings/semirings/tropical_mpolynomial + sage/rings/semirings/tropical_matrix + sage/rings/semirings/tropical_variety .. include:: ../footer.txt diff --git a/src/sage/functions/hypergeometric.py b/src/sage/functions/hypergeometric.py index 34cb77058a3..6bf1fb51f42 100644 --- a/src/sage/functions/hypergeometric.py +++ b/src/sage/functions/hypergeometric.py @@ -4,6 +4,10 @@ This module implements manipulation of infinite hypergeometric series represented in standard parametric form (as `\,_pF_q` functions). +For a more algebraic treatment of hypergeometric functions +(including reduction modulo primes and `p`-adic properties), +we refer to :mod:`sage.functions.hypergeometric_algebraic`. + AUTHORS: - Fredrik Johansson (2010): initial version @@ -301,10 +305,19 @@ def __call__(self, a, b, z, **kwargs): sage: hypergeometric([2, 3, 4], [4, 1], 1) hypergeometric((2, 3, 4), (4, 1), 1) """ - return BuiltinFunction.__call__(self, - SR._force_pyobject(a), - SR._force_pyobject(b), - z, **kwargs) + from sage.rings.polynomial.polynomial_element import Polynomial + from sage.rings.power_series_ring_element import PowerSeries + if isinstance(z, (Polynomial, PowerSeries)): + if not z.is_gen(): + raise NotImplementedError("the argument must be the generator of the polynomial ring") + S = z.parent() + from sage.functions.hypergeometric_algebraic import HypergeometricFunctions + return HypergeometricFunctions(S.base_ring(), S.variable_name())(a, b) + else: + return BuiltinFunction.__call__(self, + SR._force_pyobject(a), + SR._force_pyobject(b), + z, **kwargs) def _print_latex_(self, a, b, z): r""" diff --git a/src/sage/functions/hypergeometric_algebraic.py b/src/sage/functions/hypergeometric_algebraic.py new file mode 100644 index 00000000000..ecbd015eb4c --- /dev/null +++ b/src/sage/functions/hypergeometric_algebraic.py @@ -0,0 +1,1510 @@ +r""" +Hypergeometric functions over arbitrary rings + +[Tutorial] + +AUTHORS: + +- Xavier Caruso, Florian Fürnsinn (2025-10): initial version +""" + +# *************************************************************************** +# Copyright (C) 2025 Xavier Caruso +# Florian Fürnsinn +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 2 of the License, or +# (at your option) any later version. +# https://www.gnu.org/licenses/ +# *************************************************************************** + +import operator + +from sage.misc.latex import latex +from sage.misc.latex import latex_variable_name + +from sage.misc.misc_c import prod +from sage.misc.functional import log +from sage.functions.other import floor, ceil +from sage.functions.hypergeometric import hypergeometric +from sage.arith.misc import gcd +from sage.matrix.constructor import matrix + +from sage.structure.unique_representation import UniqueRepresentation +from sage.structure.parent import Parent +from sage.structure.element import Element +from sage.structure.element import coerce_binop +from sage.structure.category_object import normalize_names + +from sage.categories.action import Action +from sage.categories.pushout import pushout +from sage.categories.map import Map +from sage.categories.finite_fields import FiniteFields + +from sage.matrix.special import companion_matrix +from sage.matrix.special import identity_matrix +from sage.combinat.subset import Subsets + +from sage.rings.infinity import infinity +from sage.symbolic.ring import SR +from sage.sets.primes import Primes +from sage.rings.integer_ring import ZZ +from sage.rings.rational_field import QQ +from sage.rings.finite_rings.finite_field_constructor import FiniteField +from sage.rings.padics.padic_generic import pAdicGeneric +from sage.rings.padics.factory import Qp +from sage.rings.number_field.number_field import CyclotomicField + +from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing +from sage.rings.power_series_ring import PowerSeriesRing +from sage.rings.tate_algebra import TateAlgebra +from sage.rings.polynomial.ore_polynomial_ring import OrePolynomialRing + +from sage.functions.hypergeometric_parameters import HypergeometricParameters + + +# Do we want to implement polynomial linear combinaison +# of hypergeometric functions? +# Advantages: +# . reductions mod p of hypergeometric functions have this form in general +# . many methods can be extended to this context +# Difficulty: +# . not sure we can handle easily simplifications! + +class HypergeometricAlgebraic(Element): + r""" + Class for hypergeometric functions over arbitrary base rings. + """ + def __init__(self, parent, arg1, arg2=None, scalar=None): + r""" + Initialize this hypergeometric function. + + INPUT: + + - ``parent`` -- the parent of this function + + - ``arg1``, ``arg2`` -- arguments defining this hypergeometric + function, they can be: + - the top and bottom paramters + - a hypergeometric function and ``None`` + - an instance of the class :class:`HypergeometricParameters` and ``None`` + + - ``scalar`` -- an element in the base ring, the scalar by + which the hypergeometric function is multiplied + + TESTS:: + + sage: S. = QQ[] + sage: h = hypergeometric((1/2, 1/3), (1,), x) + sage: type(h) + + sage: TestSuite(h).run() + """ + Element.__init__(self, parent) + base = parent.base_ring() + if scalar is None: + scalar = base.one() + else: + scalar = base(scalar) + if scalar == 0: + parameters = None + elif isinstance(arg1, HypergeometricAlgebraic): + parameters = arg1._parameters + scalar *= base(arg1._scalar) + elif isinstance(arg1, HypergeometricParameters): + parameters = arg1 + else: + parameters = HypergeometricParameters(arg1, arg2) + char = self.parent()._char + if scalar: + if any(b in ZZ and b < 0 for b in parameters.bottom): + raise ValueError("the parameters %s do not define a hypergeometric function" % parameters) + if char > 0: + val, _ = parameters.valuation_position(char) + if val < 0: + raise ValueError("the parameters %s do not define a hypergeometric function in characteristic %s" % (parameters, char)) + self._scalar = scalar + self._parameters = parameters + self._coeffs = [scalar] + self._char = char + + def __hash__(self): + return hash((self.base_ring(), self._parameters, self._scalar)) + + def __eq__(self, other): + return (isinstance(other, HypergeometricAlgebraic) + and self.base_ring() is other.base_ring() + and self._parameters == other._parameters + and self._scalar == other._scalar) + + def _repr_(self): + if self._parameters is None: + return "0" + scalar = self._scalar + if scalar == 1: + s = "" + elif scalar._is_atomic(): + scalar = str(scalar) + if scalar == "-1": + s = "-" + else: + s = scalar + "*" + else: + s = "(%s)*" % scalar + s += "hypergeometric(%s, %s, %s)" % (self.top(), self.bottom(), self.parent().variable_name()) + return s + + def _latex_(self): + if self._parameters is None: + return "0" + scalar = self._scalar + if scalar == 1: + s = "" + elif scalar._is_atomic(): + scalar = latex(scalar) + if scalar == "-1": + s = "-" + else: + s = scalar + else: + s = r"\left(%s\right)" % scalar + top = self.top() + bottom = self.bottom() + s += r"\,_{%s} F_{%s} " % (len(top), len(bottom)) + s += r"\left(\begin{matrix} " + s += ",".join(latex(a) for a in top) + s += r"\\" + s += ",".join(latex(b) for b in bottom) + s += r"\end{matrix}; %s \right)" % self.parent().latex_variable_name() + return s + + def base_ring(self): + r""" + Return the ring over which this hypergeometric function is defined. + + EXAMPLES:: + + sage: S. = QQ[] + sage: f = hypergeometric([1/3, 2/3], [1/2], x) + sage: f.base_ring() + Rational Field + + :: + + sage: T. = Qp(5)[] + sage: g = hypergeometric([1/3, 2/3], [1/2], y) + sage: g.base_ring() + 5-adic Field with capped relative precision 20 + + :: + + sage: U. = GF(5)[] + sage: h = hypergeometric([1/3, 2/3], [1/2], z) + sage: h.base_ring() + Finite Field of size 5 + + :: + + sage: V. = CC[] + sage: k = hypergeometric([1/3, 2/3], [1/2], w) + sage: k.base_ring() + Complex Field with 53 bits of precision + """ + return self.parent().base_ring() + + def top(self): + r""" + Return the top parameters of this hypergeometric function. + + EXAMPLES:: + + sage: S. = QQ[] + sage: f = hypergeometric([1/3, 2/3], [1/2], x) + sage: f.top() + (1/3, 2/3) + """ + return self._parameters.top + + def bottom(self): + r""" + Return the bottom parameters of this hypergeometric function (excluding + the extra ``1``). + + EXAMPLES:: + + sage: S. = QQ[] + sage: f = hypergeometric([1/3, 2/3], [1/2], x) + sage: f.bottom() + (1/2,) + """ + return self._parameters.bottom[:-1] + + def scalar(self): + r""" + Return the scalar of this hypergeometric function. + + EXAMPLES:: + + sage: S. = QQ[] + sage: f = hypergeometric([1/3, 2/3], [1/2], x) + sage: f.scalar() + 1 + sage: g = 4*f + sage: g.scalar() + 4 + """ + return self._scalar + + def change_ring(self, R): + r""" + Return this hypergeometric function with changed base ring. + + INPUT: + + - ``R`` -- a commutative ring + + EXAMPLES:: + + sage: S. = QQ[] + sage: f = hypergeometric([1/3, 2/3], [1/2], x) + sage: f.base_ring() + Rational Field + sage: g = f.change_ring(Qp(5)) + sage: g.base_ring() + 5-adic Field with capped relative precision 20 + """ + H = self.parent().change_ring(R) + return H(self._parameters, None, self._scalar) + + def change_variable_name(self, name): + r""" + Return this hypergeometric function with changed variable name + + INPUT: + + - ``name`` -- a string, the new variable name + + EXAMPLES:: + + sage: S. = Qp(5)[] + sage: T. = Qp(5)[] + sage: f = hypergeometric([1/3, 2/3], [1/2], x) + sage: f + hypergeometric((1/3, 2/3), (1/2,), x) + sage: g = f.change_variable_name('y') + sage: g + hypergeometric((1/3, 2/3), (1/2,), y) + """ + H = self.parent().change_variable_name(name) + return H(self._parameters, None, self._scalar) + + def _add_(self, other): + r""" + Return the (formal) sum of the hypergeometric function + and ``other``. + + INPUT: + + - ``other`` -- a hypergeometric function + + EXAMPLES:: + + sage: S. = QQ[] + sage: f = hypergeometric([1/3, 2/3], [1/2], x) + sage: g = 1/2 * hypergeometric([1/3, 2/3], [1/2], x) + sage: h = hypergeometric([1/5, 2/5], [3/5], x) + sage: f + g + 3/2*hypergeometric((1/3, 2/3), (1/2,), x) + sage: f + h + hypergeometric((1/3, 2/3), (1/2,), x) + hypergeometric((1/5, 2/5), (3/5,), x) + + :: + + sage: f + cos(x) + cos(x) + hypergeometric((1/3, 2/3), (1/2,), x) + """ + if self._parameters is None: + return other + if isinstance(other, HypergeometricAlgebraic): + if other._parameters is None: + return self + if self._parameters == other._parameters: + scalar = self._scalar + other._scalar + return self.parent()(self._parameters, scalar=scalar) + return SR(self) + SR(other) + + def _neg_(self): + r""" + Return the negative of this hypergeometric function. + + EXAMPLES:: + + sage: S. = QQ[] + sage: f = 2*hypergeometric([1/3, 2/3], [1/2], x) + sage: -f + -2*hypergeometric((1/3, 2/3), (1/2,), x) + """ + if self._parameters is None: + return self + return self.parent()(self._parameters, scalar=-self._scalar) + + def _sub_(self, other): + r""" + Return the (formal) difference of the hypergeometric function + with ``other``. + + INPUT: + + - ``other`` -- a hypergeometric function or a formal expression + + EXAMPLES:: + + sage: S. = QQ[] + sage: f = hypergeometric([1/3, 2/3], [1/2], x) + sage: g = 1/2 * hypergeometric([1/3, 2/3], [1/2], x) + sage: h = hypergeometric([1/5, 2/5], [3/5], x) + sage: f - g + 1/2*hypergeometric((1/3, 2/3), (1/2,), x) + sage: f - h + hypergeometric((1/3, 2/3), (1/2,), x) - hypergeometric((1/5, 2/5), (3/5,), x) + + :: + + sage: f - sin(x) + hypergeometric((1/3, 2/3), (1/2,), x) - sin(x) + """ + if self._parameters is None: + return other + if isinstance(other, HypergeometricAlgebraic): + if other._parameters is None: + return self + if self._parameters == other._parameters: + scalar = self._scalar - other._scalar + return self.parent()(self._parameters, scalar=scalar) + return SR(self) - SR(other) + + def _mul_(self, other): + r""" + Return the (formal) product of the hypergeometric function + and ``other`` + + INPUT: + + - ``other`` -- a hypergeometric function or a formal expression + + EXAMPLES:: + + sage: S. = QQ[] + sage: f = hypergeometric([1/3, 2/3], [1/2], x) + sage: g = 1/2 * hypergeometric([1/3, 2/3], [1/2], x) + sage: h = hypergeometric([1/5, 2/5], [3/5], x) + sage: f*g + 1/2*hypergeometric((1/3, 2/3), (1/2,), x)^2 + sage: f*h + hypergeometric((1/3, 2/3), (1/2,), x)*hypergeometric((1/5, 2/5), (3/5,), x) + + :: + + sage: sin(x)*f + x + hypergeometric((1/3, 2/3), (1/2,), x)*sin(x) + x + """ + return SR(self) * SR(other) + + def __call__(self, x): + r""" + Return the value of this hypergeometric function at ``x``. + + INPUT: + + - ``x`` -- an element + + EXAMPLES:: + + sage: S. = RR[] + sage: f = hypergeometric([1/3, 2/3], [1/2], x) + sage: f(0.5) + 1.36602540378444 + + :: + + sage: g = 2*f + sage: g(0.2) + 2.20941633798502 + """ + scalar = self._scalar + if scalar == 0: + return self.base_ring().zero() + X = SR('X') + h = hypergeometric(self.top(), self.bottom(), X) + if scalar != 1: + h *= scalar + return h(X=x) + + def _compute_coeffs(self, prec): + r""" + Compute the coefficients of the series representation of this + hypergeometric function up to a given precision, and store + them in ``self._coeffs``. + + INPUT: + + - ``prec`` -- a positive integer + + EXAMPLES:: + + sage: S. = QQ[] + sage: f = hypergeometric([1/3, 2/3], [1/2], x) + sage: f._coeffs + [1] + sage: f._compute_coeffs(3) + sage: f._coeffs + [1, 4/9, 80/243] + """ + coeffs = self._coeffs + start = len(coeffs) - 1 + c = coeffs[-1] + for i in range(start, prec - 1): + for a in self._parameters.top: + c *= a + i + for b in self._parameters.bottom: + c /= b + i + coeffs.append(c) + + def nth_coefficient(self, n): + r""" + Return the ``n``-th coefficient of the series representation of this + hypergeoimetric function. + + INPUT: + + - ``n`` -- a non-negative integer + + EXAMPLES: + + sage: S. = QQ[] + sage: f = hypergeometric([1/3, 2/3], [1/2], x) + sage: f.nth_coefficient(9) + 409541017600/2541865828329 + sage: g = f % 5 + sage: g.nth_coefficient(9) + 0 + """ + self._compute_coeffs(n+1) + S = self.base_ring() + return S(self._coeffs[n]) + + def power_series(self, prec=20): + r""" + Return the power series representation of this hypergeometric + function up to a given precision. + + INPUT: + + - ``prec`` -- a positive integer + + EXAMPLES:: + + sage: S. = QQ[] + sage: f = hypergeometric([1/3, 2/3], [1/2], x) + sage: f.power_series(3) + 1 + 4/9*x + 80/243*x^2 + O(x^3) + """ + S = self.parent().power_series_ring() + self._compute_coeffs(prec) + return S(self._coeffs, prec=prec) + + def shift(self, s): + r""" + Return this hypergeometric function, where each parameter + (including the additional ``1`` as a bottom parameter) is + increased by ``s``. + + INPUT: + + - ``s`` -- a rational number + + EXAMPLES:: + + sage: S. = QQ[] + sage: f = hypergeometric([1/3, 2/3], [1/2], x) + sage: g = f.shift(3/2) + sage: g + hypergeometric((1, 11/6, 13/6), (2, 5/2), x) + """ + return self.parent()(self._parameters.shift(s), scalar=self._scalar) + + @coerce_binop + def hadamard_product(self, other): + r""" + Return the hadamard product of this hypergeometric function + and ``other``. + + INPUT: + + - ``other`` -- a hypergeometric function + + EXAMPLES:: + + sage: S. = QQ[] + sage: f = hypergeometric([1/3, 2/3], [1/2], x) + sage: h = 1/2*hypergeometric([1/5, 2/5], [3/5], x) + sage: f.hadamard_product(h) + 1/2*hypergeometric((1/5, 1/3, 2/5, 2/3), (1/2, 3/5, 1), x) + """ + if self._scalar == 0: + return self + if other._scalar == 0: + return other + top = self.top() + other.top() + bottom = self._parameters.bottom + other.bottom() + scalar = self._scalar * other._scalar + return self.parent()(top, bottom, scalar=scalar) + + def _div_(self, other): + r""" + Return the (formal) quotient of the hypergeometric function + and ``other``. + + INPUT: + + - ``other`` -- a hypergeometric function or a formal expression + + EXAMPLES:: + + sage: S. = QQ[] + sage: f = hypergeometric([1/3, 2/3], [1/2], x) + sage: g = 1/2 * hypergeometric([1/3, 2/3], [1/2], x) + sage: h = hypergeometric([1/5, 2/5], [3/5], x) + sage: f/g + 2 + sage: f/h + hypergeometric((1/3, 2/3), (1/2,), x)/hypergeometric((1/5, 2/5), (3/5,), x) + + :: + + sage: f/sin(x) + x + x + hypergeometric((1/3, 2/3), (1/2,), x)/sin(x) + """ + return SR(self) / SR(other) + + def denominator(self): + r""" + Return the smallest common denominator of the parameters. + + EXAMPLES:: + + sage: S. = QQ[] + sage: f = hypergeometric([1/3, 2/3], [1/2], x) + sage: f.denominator() + 6 + """ + return self._parameters.d + + def differential_operator(self, var='d'): + r""" + Return the hypergeometric differential operator that annihilates + this hypergeometric function as an Ore polynomial in the variable + ``var``. + + INPUT: + + - ``var`` -- a string (default: ``d``), the variable name of + the derivation + + EXAMPLES:: + + sage: S. = QQ[] + sage: f = hypergeometric([1/3, 2/3], [1/2], x) + sage: f.differential_operator(var='D') + (-x^2 + x)*D^2 + (-2*x + 1/2)*D - 2/9 + + Note that this does not necessarily give the minimal differential + operator annihilating this hypergeometric function: in the example + below, this method returns an operator of order `3` where `g` is + solution of a differential equation of order `2`:: + + sage: g = hypergeometric([1/3, 2/3, 6/5], [1/5, 1/2], x) + sage: L = g.differential_operator() + sage: L.degree() + 3 + sage: gs = g.power_series(100) + sage: (72*x^3 - 234*x^2 + 162*x)*gs.derivative(2) + (144*x^2 - 450*x + 81)*gs.derivative() + (16*x - 216)*gs + O(x^99) + """ + S = self.parent().polynomial_ring() + x = S.gen() + D = OrePolynomialRing(S, S.derivation(), names=var) + if self._scalar == 0: + return D.one() + t = x * D.gen() + A = D.one() + for a in self._parameters.top: + A *= t + S(a) + B = D.one() + for b in self._parameters.bottom: + B *= t + S(b-1) + L = B - x*A + return D([c//x for c in L.list()]) + + def derivative(self): + r""" + Return the derivative of this hypergeometric function. + + EXAMPLES:: + + sage: S. = QQ[] + sage: f = hypergeometric([1/3, 2/3], [1/2], x) + sage: f.derivative() + 4/9*hypergeometric((4/3, 5/3), (3/2,), x) + """ + top = [a+1 for a in self.top()] + bottom = [b+1 for b in self.bottom()] + scalar = prod(self._parameters.top) / prod(self._parameters.bottom) + scalar = self.base_ring()(scalar) * self._scalar + return self.parent()(top, bottom, scalar) + + +# Over the rationals + +class HypergeometricAlgebraic_QQ(HypergeometricAlgebraic): + def __mod__(self, p): + r""" + Return the reduction of the hypergeometric function modulo ``p``. + + INPUT: + + - ``p`` -- a prime number. + + EXAMPLES:: + + sage: S. = QQ[] + sage: f = hypergeometric([1/3, 2/3], [1/2], x) + sage: g = f % 5 + sage: g + hypergeometric((1/3, 2/3), (1/2,), x) + sage: g.base_ring() + Finite Field of size 5 + """ + k = FiniteField(p) + val = self._scalar.valuation(p) + if val == 0: + return self.change_ring(k) + h = self.change_ring(Qp(p, 1)) + return h.residue() + + def valuation(self, p, position=False): + r""" + Return the `p`-adic valuation of this hypergeometric function, i.e., the + maximal `s`, such that `p^{-s}` times this hypergeometric function has + p-integral coefficients. + + INPUT: + + - ``p`` -- a prime number + + - ``position`` -- a boolean (default: ``False``); if ``True``, return + also the first index in the series expansion at which the valuation + is attained. + + EXAMPLES:: + + sage: S. = QQ[x] + sage: f = hypergeometric([1/3, 2/3], [1/2], x) + sage: f.valuation(5) + 0 + sage: g = 5*f + sage: g.valuation(5) + 1 + + An example where we ask for the position:: + + sage: h = hypergeometric([1/5, 1/5, 1/5], [1/3, 9/5], x) + sage: h.valuation(3, position=True) + (-1, 1) + + We can check that the coefficient in `x` in the series expansion + has indeed valuation `-1`:: + + sage: s = h.power_series() + sage: s + 1 + 1/75*x + 27/8750*x^2 + ... + O(x^20) + sage: s[1].valuation(3) + -1 + + TESTS:: + + sage: g.valuation(9) + Traceback (most recent call last): + ... + ValueError: p must be a prime number + """ + if not p.is_prime(): + raise ValueError("p must be a prime number") + val, pos = self._parameters.valuation_position(p) + val += self._scalar.valuation(p) + if position: + return val, pos + else: + return val + + def has_good_reduction(self, p): + r""" + Return whether the `p`-adic valuation of this hypergeometric + function is nonnegative, i.e., if its reduction modulo ``p`` + is well-defined. + + INPUT: + + - ``p`` -- a prime number + + EXAMPLES:: + + sage: S. = QQ[x] + sage: f = hypergeometric([1/3, 2/3], [1/2], x) + sage: f.valuation(5) + 0 + sage: f.has_good_reduction(5) + True + sage: g = 1/5*f + sage: g.has_good_reduction(5) + False + """ + return self.valuation(p) >= 0 + + def good_reduction_primes(self): + r""" + Return the set of prime numbers modulo which this hypergeometric + function can be reduced, i.e., the p-adic valuation is nonnegative. + + EXAMPLES:: + + sage: S. = QQ[] + sage: f = hypergeometric([1/3, 2/3], [1/2], x) + sage: f.good_reduction_primes() + Set of all prime numbers with 3 excluded: 2, 5, 7, 11, ... + + ALGORITHM: + + We rely on Christol's criterion ([Chr1986]_, Prop. 1) for globally + bounded hypergeometric function, from which a criterion can be deduced + modulo which primes a hypergeometric function can be reduced + ([CFV2025]_, Thm. 3.1.3). For small primes `p`, we compute the `p`-adic + valuation of the hypergeometric function individually. + """ + params = self._parameters + d = params.d + + # We check the parenthesis criterion for c=1 + if not params.parenthesis_criterion(1): + return Primes(modulus=0) + + # We check the parenthesis criterion for other c + # and derive congruence classes with good reduction + cs = [c for c in range(d) if d.gcd(c) == 1] + goods = {c: None for c in cs} + goods[1] = True + for c in cs: + if goods[c] is not None: + continue + cc = c + goods[c] = True + while cc != 1: + if goods[cc] is False or not params.parenthesis_criterion(cc): + goods[c] = False + break + cc = (cc * c) % d + if goods[c]: + cc = c + while cc != 1: + goods[cc] = True + cc = (cc * c) % d + + # We treat exceptional primes + bound = params.bound + exceptions = {} + for p in Primes(): + if p > bound: + break + if d % p == 0 and self.valuation(p) >= 0: + exceptions[p] = True + if d % p == 0 or not goods[p % d]: + continue + if self.valuation(p) < 0: + exceptions[p] = False + + goods = [c for c, v in goods.items() if v] + return Primes(modulus=d, classes=goods, exceptions=exceptions) + + def is_algebraic(self): + r""" + Return ``True`` if this hypergeometric function is algebraic over + the rational functions, return ``False`` otherwise. + + EXAMPLES:: + + sage: S. = QQ[] + sage: f = hypergeometric([1/3, 2/3], [1/2], x) + sage: f.is_algebraic() + True + sage: g = hypergeometric([1/3, 2/3, 1/4], [5/4, 1/2], x) + sage: g.is_algebraic() + False + + ALGORITHM: + + We rely on the (Christol-)Beukers-Heckmann interlacing criterion + (see [Chr1986]_, p.15, Cor.; [BeukersHeckman]_, Thm. 4.5). For integer + differences between parameters we follow the flowchart in + [FY2024]_, Fig. 1. + """ + if any(a in ZZ and a <= 0 for a in self.top()): + return True + if not self._parameters.is_balanced(): + return False + simplified_parameters = self._parameters.remove_positive_integer_differences() + if simplified_parameters.has_negative_integer_differences(): + return False + d = simplified_parameters.d + return all(simplified_parameters.interlacing_criterion(c) + for c in range(d) if d.gcd(c) == 1) + + def is_globally_bounded(self, include_infinity=True): + r""" + Return whether this hypergeometric function is globally bounded + (if ``include_infinity`` is ``False`` it is not checked whether + the radius of convergence is finite). + + INPUT: + + - ``include_infinity`` -- Boolean (default: ``True``) + + EXAMPLES: + + sage: S. = QQ[] + sage: f = hypergeometric([1/9, 4/9, 5/9], [1/3, 1], x) + sage: f.is_globally_bounded() + True + sage: g = hypergeometric([1/9, 4/9, 5/9], [1/3], x) + sage: g.is_globally_bounded() + False + sage: g.is_globally_bounded(include_infinity=False) + True + + ALGORITHM: + + We rely on Christol's classification of globally bounded hypergeometric + functions (see [Chr1986]_, Prop. 1). + """ + if include_infinity and len(self.top()) > len(self.bottom()) + 1: + return False + d = self.denominator() + for c in range(d): + if d.gcd(c) == 1: + if not self._parameters.parenthesis_criterion(c): + return False + return True + + def p_curvature_ranks(self): + # Should this return the coranks of the p-curvature depending on the + # congruence class of a prime? + # YES! + raise NotImplementedError + + def monodromy(self, x=0, var='z'): + r""" + Return a local monodromy matrix of the hypergeometric differential + equation associated to this hypergeoemtric function at the point + ``x``. + + INPUT: + + - ``x`` -- a complex number (default: ``0``) + + - ``var`` -- a string (default: ``z``), the name of the variable + representing a `d`-th root of unity for `d` being least common + multiple of the parameters. + + EXAMPLES:: + + sage: S. = QQ[] + sage: f = hypergeometric([1/3, 2/3], [1/2], x) + sage: f.monodromy() + [0 1] + [1 0] + + :: + + The bases of the solution space are chosen in a compatible way + across the three singularities of the differential equation: + + sage: g = hypergeometric([1/9, 4/9, 5/9], [1/3, 1], x) + sage: g.monodromy(var='a') + [ -a^3 + 1 1 0] + [2*a^3 + 1 0 1] + [ -a^3 - 1 0 0] + sage: g.monodromy(x=Infinity) * g.monodromy(x=1) * g.monodromy() + [1 0 0] + [0 1 0] + [0 0 1] + + ALGORITHM: + + We use the explicit formulas for the monodromy matrices presented in + [BeukersHeckman]_, Thm. 3.5, attributed to Levelt. + """ + params = self._parameters + if not params.is_balanced(): + raise ValueError("hypergeometric equation is not Fuchsian") + d = params.d + K = CyclotomicField(d, names=var) + z = K.gen() + S = PolynomialRing(K, names='X') + X = S.gen() + if x == 0: + B = prod(X - z**(b*d) for b in params.bottom) + return companion_matrix(B, format='right').inverse() + elif x == 1: + A = prod(X - z**(a*d) for a in params.top) + B = prod(X - z**(b*d) for b in params.bottom) + return companion_matrix(A, format='right').inverse() * companion_matrix(B, format='right') + elif x is infinity: + A = prod(X - z**(a*d) for a in params.top) + return companion_matrix(A, format='right') + else: + n = len(params.top) + return identity_matrix(QQ, n) + + def is_maximum_unipotent_monodromy(self): + r""" + Return whether the hypergeometric differential operator associated + to this hypergeometric function has maximal unipotent monodromy (MUM). + + EXAMPLES:: + + sage: S. = QQ[] + sage: f = hypergeometric([1/3, 2/3], [1/2], x) + sage: f.is_maximum_unipotent_monodromy() + False + sage: g = hypergeometric([1/9, 4/9, 5/9], [1, 2], x) + sage: g.is_maximum_unipotent_monodromy() + True + """ + return all(b in ZZ for b in self.bottom()) + + is_mum = is_maximum_unipotent_monodromy + + +# Over the p-adics + +class HypergeometricAlgebraic_padic(HypergeometricAlgebraic): + def __init__(self, parent, arg1, arg2=None, scalar=None): + r""" + Initialize this hypergeometric function. + + INPUT: + + - ``parent`` -- the parent of this function, which has to be defined + over the p-adics + + - ``arg1``, ``arg2`` -- arguments defining this hypergeometric + function, they can be: + - the top and bottom paramters + - a hypergeometric function and ``None`` + - an instance of the class :class:`HypergeometricParameters` and ``None`` + + - ``scalar`` -- an element in the base ring, the scalar by + which the hypergeometric function is multiplied + + TESTS:: + + sage: S. = Qp(5, 3)[] + sage: h = hypergeometric((1/2, 1/3), (1,), x) + sage: type(h) + + sage: TestSuite(h).run() + """ + HypergeometricAlgebraic.__init__(self, parent, arg1, arg2, scalar) + K = self.base_ring() + self._p = K.prime() + self._e = K.e() + + def residue(self): + r""" + Return the reduction of this hypergeometric function in the residue + field of the p-adics over which this hypergeometric function is + defined. + + EXAMPLES:: + + sage: S. = Qp(5)[] + sage: f = hypergeometric([1/3, 2/3], [1/2], x) + sage: f.parent() + Hypergeometric functions in x over 5-adic Field with capped relative precision 20 + sage: g = f.residue() + sage: g.parent() + Hypergeometric functions in x over Finite Field of size 5 + """ + k = self.base_ring().residue_field() + if self._scalar.valuation() == 0: + return self.change_ring(k) + val, pos = self._parameters.valuation_position(self._p) + if val < 0: + raise ValueError("bad reduction") + if val > 0: + H = self.parent().change_ring(k) + return H(self._parameters, scalar=0) + raise NotImplementedError("the reduction is not a hypergeometric function") + # In fact, it is x^s * h[s] * h, with + # . s = pos + # . h = self.shift(s) + + def dwork_image(self): + r""" + Return the hypergeometric function obtained from this one + by applying the Dwork map to each of its parameters. + + EXAMPLES:: + + sage: S. = QQ[] + sage: f = hypergeometric([1/4, 1/3, 1/2], [2/5, 3/5, 1], x) + sage: f.dwork_image() + hypergeometric((1/3, 1/2, 3/4), (1/5, 4/5, 1), x) + """ + parameters = self._parameters.dwork_image(self._p) + return self.parent()(parameters, scalar=self._scalar) + + def log_radius_of_convergence(self): + r""" + Return the logarithmic p-adic radius of convergence of this + hypergeometric function. + + EXAMPLES:: + + sage: S. = Qp(5)[] + sage: f = hypergeometric([1/3, 2/3], [1/2], x) + sage: f.log_radius_of_convergence() + 0 + sage: g = hypergeometric([1/3, 2/3], [1/5], x) + sage: g.log_radius_of_convergence() + 5/4 + """ + p = self._p + step = self._e / (p - 1) + log_radius = 0 + for a in self._parameters.top: + if a in ZZ and a <= 0: + return infinity + v = a.valuation(p) + if v < 0: + log_radius += v + else: + log_radius += step + for b in self._parameters.bottom: + v = b.valuation(p) + if v < 0: + log_radius -= v + else: + log_radius -= step + return log_radius + + def valuation(self, log_radius=0, position=False): + r""" + Return the p-adic valuation of this hypergeometric function on the + disk of logarithmic radius ``log_radius``, and, if ``position`` is + ``True`` the index of the first coefficient of the series that + attains this valuation. + + INPUT: + + - ``log_radius`` -- a rational number + + - ``position`` -- a boolean (default: ``False``), if ``True`` the + index of the first coefficient attaining the valuation is also + returned + + EXAMPLES:: + + sage: S. = Qp(5)[] + sage: f = hypergeometric([1/3, 2/3], [1/2], x) + sage: f.valuation() + 0 + + :: + + sage: S. = Qp(5)[] + sage: g = 1/5 * hypergeometric([1/3, 2/3], [5^3/3], x) + sage: g.valuation(-1, position=True) + (-2, 1) + """ + drift = -log_radius / self._e + val, pos = self._parameters.valuation_position(self._p, drift) + if position: + return val, pos + else: + return val + + def newton_polygon(self, log_radius): + raise NotImplementedError + + def _truncation_bound(self, log_radius, prec): + convergence = self.log_radius_of_convergence() + margin = convergence - log_radius + if margin <= 0: + raise ValueError("outside the domain of convergence") + val = self.valuation(convergence) + if val is not -infinity: + lr = convergence + else: + # We choose an intermediate log_radius + # It can be anything between convergence and log_radius + # but it seems that the following works well (in the sense + # that it gives good bounds at the end). + lr = convergence - margin / max(prec, 2) + val = self.valuation(lr) + # Now, we know that + # val(h_k) >= -lr*k + val + # and we want to find k such that + # val(h_k) >= -log_radius*k + prec + # So we just solve the equation. + k = (prec - val) / (lr - log_radius) + return 1 + max(0, floor(k)) + + def tate_series(self, log_radius, prec=None): + K = self.base_ring() + name = self.parent().variable_name() + S = TateAlgebra(K, log_radii=[log_radius], names=name) + if prec is None: + prec = self.base_ring().precision_cap() + trunc = self._truncation_bound(log_radius, prec) + self._compute_coeffs(trunc) + coeffs = {(i,): self._coeffs[i] for i in range(trunc)} + return self._scalar * S(coeffs, prec) + + def __call__(self, x): + K = self.base_ring() + x = K(x) + val = min(x.valuation(), x.precision_absolute()) + if val is infinity: + return K.one() + w = self.valuation(-val) + prec = w + K.precision_cap() + trunc = self._truncation_bound(-val, prec) + self._compute_coeffs(trunc) + ans = sum(self._coeffs[i] * x**i for i in range(trunc)) + ans = ans.add_bigoh(prec) + return self._scalar * ans + + +# Over prime finite fields + +class HypergeometricAlgebraic_GFp(HypergeometricAlgebraic): + def __init__(self, parent, arg1, arg2=None, scalar=None): + # TODO: do we want to simplify automatically if the + # hypergeometric series is a polynomial? + HypergeometricAlgebraic.__init__(self, parent, arg1, arg2, scalar) + self._p = p = self.base_ring().cardinality() + self._coeffs = [Qp(p, 1)(self._scalar)] + + def power_series(self, prec=20): + S = self.parent().power_series_ring() + self._compute_coeffs(prec) + try: + f = S(self._coeffs, prec=prec) + except ValueError: + raise ValueError("denominator appears in the series at the required precision") + return f + + def is_almost_defined(self): + p = self._char + d = self.denominator() + if d.gcd(p) > 1: + return False + u = 1 + if not self._parameters.parenthesis_criterion(u): + return False + u = p % d + while u != 1: + if not self._parameters.parenthesis_criterion(u): + return False + u = p*u % d + return True + + def is_defined(self): + p = self._char + if not self.is_almost_defined(): + return False + bound = self._parameters.bound + if bound < p: + return True + prec = 1 + p ** ceil(log(self._parameters.bound, p)) + try: + self.series(prec) + except ValueError: + return False + return True + + def is_defined_conjectural(self): + p = self._char + if not self.is_almost_defined(): + return False + bound = self._parameters.bound + q = p + while q <= bound: + if not self._parameters.q_parenthesis_criterion(q): + return False + q *= p + return True + + def __call__(self, x): + return self.polynomial()(x) + + def is_polynomial(self): + raise NotImplementedError + + def degree(self): + raise NotImplementedError + + def polynomial(self): + raise NotImplementedError + + def is_algebraic(self): + return True + + def p_curvature(self): + L = self.differential_operator() + K = L.base_ring().fraction_field() + S = OrePolynomialRing(K, L.parent().twisting_derivation().extend_to_fraction_field(), names='d') + L = S(L.list()) + d = S.gen() + p = self._char + rows = [ ] + n = L.degree() + for i in range(p, p + n): + Li = d**i % L + rows.append([Li[j] for j in range(n)]) + return matrix(rows) + + def p_curvature_corank(self): # maybe p_curvature_rank is preferable? + # TODO: check if it is also correct when the parameters are not balanced + return self._parameters.q_interlacing_number(self._char) + + def dwork_relation(self): + r""" + Return (P1, h1), ..., (Ps, hs) such that + + self = P1*h1^p + ... + Ps*hs^p + """ + parameters = self._parameters + if not parameters.is_balanced(): + raise ValueError("the hypergeometric function must be a pFq with q = p-1") + p = self._char + H = self.parent() + F = H.base_ring() + Hp = H.change_ring(Qp(p, 1)) + x = H.polynomial_ring().gen() + coeffs = self._coeffs + Ps = {} + for r in range(p): + params = parameters.shift(r).dwork_image(p) + _, s = params.valuation_position(p) + h = H(params.shift(s)) + e = s*p + r + if e >= len(coeffs): + self._compute_coeffs(e + 1) + c = F(coeffs[e]) + if c: + if h in Ps: + Ps[h] += c * x**e + else: + Ps[h] = c * x**e + return Ps + + def annihilating_ore_polynomial(self, var='Frob'): + # QUESTION: does this method actually return the + # minimal Ore polynomial annihilating self? + # Probably not :-( + parameters = self._parameters + if not parameters.is_balanced(): + raise NotImplementedError("the hypergeometric function is not a pFq with q = p-1") + + p = self._char + S = self.parent().polynomial_ring() + zero = S.zero() + Frob = S.frobenius_endomorphism() + Ore = OrePolynomialRing(S, Frob, names=var) + + # We remove the scalar + if self._scalar == 0: + return Ore.one() + self = self.parent()(parameters) + + order = parameters.frobenius_order(p) + bound = self.p_curvature_corank() + + rows = [{self: S.one()}] + # If row is the i-th item of rows, we have: + # self = sum_g row[g] * g**(p**i) + q = 1 + while True: + row = {} + previous_row = rows[-1] + for _ in range(order): + row = {} + for g, P in previous_row.items(): + for h, Q in g.dwork_relation().items(): + # here g = sum(Q * h^p) + if h in row: + row[h] += P * insert_zeroes(Q, q) + else: + row[h] = P * insert_zeroes(Q, q) + previous_row = row + q *= p # q = p**i + rows.append(row) + + i = len(rows) + Mrows = [] + Mqo = 1 + columns = {} + for j in range(i-1, max(-1, i-2-bound), -1): + for col in rows[j]: + columns[col] = None + for j in range(i-1, max(-1, i-2-bound), -1): + Mrow = [] + for col in columns: + Mrow.append(insert_zeroes(rows[j].get(col, zero), Mqo)) + Mrows.append(Mrow) + Mqo *= p ** order + M = matrix(S, Mrows) + + ker = kernel(M) + if ker is not None: + return insert_zeroes(Ore(ker), order) + + def is_lucas(self): + p = self._char + if self._parameters.frobenius_order(p) > 1: + # TODO: check this + return False + S = self.parent().polynomial_ring() + K = S.fraction_field() + Ore = OrePolynomialRing(K, K.frobenius_endomorphism(), names='F') + Z = Ore(self.annihilating_ore_polynomial()) + Ap = self.series(p).polynomial() + F = Ap * Ore.gen() - 1 + return (Z % F).is_zero() + + +# Parent +######## + +class HypergeometricToSR(Map): + def _call_(self, h): + return h.scalar() * hypergeometric(h.top(), h.bottom(), SR.var(h.parent().variable_name())) + + +class ScalarMultiplication(Action): + def _act_(self, scalar, h): + return h.parent()(h, scalar=scalar) + + +class HypergeometricFunctions(Parent, UniqueRepresentation): + def __init__(self, base, name, category=None): + self._name = normalize_names(1, name)[0] + self._latex_name = latex_variable_name(self._name) + self._char = char = base.characteristic() + if char == 0: + base = pushout(base, QQ) + if base in FiniteFields() and base.is_prime_field(): + self.Element = HypergeometricAlgebraic_GFp + elif base is QQ: + self.Element = HypergeometricAlgebraic_QQ + elif isinstance(base, pAdicGeneric): + self.Element = HypergeometricAlgebraic_padic + else: + self.Element = HypergeometricAlgebraic + Parent.__init__(self, base, category=category) + self.register_action(ScalarMultiplication(base, self, False, operator.mul)) + self.register_action(ScalarMultiplication(base, self, True, operator.mul)) + if char == 0: + SR.register_coercion(HypergeometricToSR(self.Hom(SR))) + + def _repr_(self): + return "Hypergeometric functions in %s over %s" % (self._name, self._base) + + def _element_constructor_(self, *args, **kwds): + return self.element_class(self, *args, **kwds) + + def _coerce_map_from_(self, other): + if (isinstance(other, HypergeometricFunctions) + and other.has_coerce_map_from(self)): + return True + + def _pushout_(self, other): + if isinstance(other, HypergeometricFunctions) and self._name == other._name: + base = pushout(self.base_ring(), other.base_ring()) + if base is not None: + return HypergeometricFunctions(base, self._name) + if SR.has_coerce_map_from(other): + return SR + + def base_ring(self): + return self._base + + def variable_name(self): + return self._name + + def latex_variable_name(self): + return self._latex_name + + def change_ring(self, R): + return HypergeometricFunctions(R, self._name) + + def change_variable_name(self, name): + return HypergeometricFunctions(self._base, name) + + def polynomial_ring(self): + return PolynomialRing(self.base_ring(), self._name) + + def power_series_ring(self, default_prec=None): + return PowerSeriesRing(self.base_ring(), self._name, default_prec=default_prec) + + +# Helper functions +################## + +def insert_zeroes(P, n): + cs = P.list() + coeffs = n * len(cs) * [0] + for i in range(len(cs)): + coeffs[n*i] = cs[i] + return P.parent()(coeffs) + + +def kernel(M, repeat=2): + n = M.nrows() + m = M.ncols() + if n > m + 1: + raise RuntimeError + if n <= m: + K = M.base_ring().base_ring() + for _ in range(repeat): + a = K.random_element() + Me = matrix(n, m, [f(a) for f in M.list()]) + if Me.rank() == n: + return + for J in Subsets(range(m), n-1): + MJ = M.matrix_from_columns(J) + minor = MJ.delete_rows([0]).determinant() + if minor.is_zero(): + continue + ker = [minor] + for i in range(1, n): + minor = MJ.delete_rows([i]).determinant() + ker.append((-1)**i * minor) + Z = matrix(ker) * M + if not Z.is_zero(): + return + g = ker[0].leading_coefficient() * gcd(ker) + ker = [c//g for c in ker] + return ker diff --git a/src/sage/functions/hypergeometric_parameters.py b/src/sage/functions/hypergeometric_parameters.py new file mode 100644 index 00000000000..310fc59eb7c --- /dev/null +++ b/src/sage/functions/hypergeometric_parameters.py @@ -0,0 +1,730 @@ +r""" +Parameters for hypergeometric functions + +AUTHORS: + +- Xavier Caruso, Florian Fürnsinn (2025-10): initial version +""" + +# *************************************************************************** +# Copyright (C) 2025 Xavier Caruso +# Florian Fürnsinn +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 2 of the License, or +# (at your option) any later version. +# https://www.gnu.org/licenses/ +# *************************************************************************** + +from sage.misc.cachefunc import cached_method + +from sage.functions.other import ceil +from sage.arith.functions import lcm + +from sage.rings.infinity import infinity +from sage.rings.integer_ring import ZZ +from sage.rings.rational_field import QQ +from sage.rings.finite_rings.integer_mod_ring import IntegerModRing +from sage.rings.semirings.tropical_semiring import TropicalSemiring + +from sage.matrix.constructor import matrix +from sage.matrix.special import identity_matrix + + +# HypergeometricParameters of hypergeometric functions +######################################## + +class HypergeometricParameters(): + r""" + Class for parameters of hypergeometric functions. + """ + def __init__(self, top, bottom, add_one=True): + r""" + Initialize this set of parameters. + + INPUT: + + - ``top`` -- list of top parameters + + - ``bottom`` -- list of bottom parameters + + - ``add_one`` -- boolean (default: ``True``), + if ``True``, add an additional one to the bottom + parameters. + + EXAMPLES:: + + sage: from sage.functions.hypergeometric_algebraic import HypergeometricParameters + sage: pa = HypergeometricParameters([1/4, 1/3, 1/2], [2/5, 3/5]) + sage: pa + ((1/4, 1/3, 1/2), (2/5, 3/5, 1)) + sage: type(pa) + + + By default, parameters are sorted, duplicates are removed and + a trailing `1` is added to the bottom parameters:: + + sage: HypergeometricParameters([1/2, 1/3, 2/3], [2/3]) + ((1/3, 1/2), (1,)) + + We can avoid adding the trailing `1` by passing ``add_one=False``:: + + sage: HypergeometricParameters([1/2, 1/3, 2/3], [2/3], add_one=False) + ((1/3, 1/2, 1), (1,)) + """ + try: + top = sorted([QQ(a) for a in top if a is not None]) + bottom = sorted([QQ(b) for b in bottom if b is not None]) + except TypeError: + raise NotImplementedError("parameters must be rational numbers") + i = j = 0 + while i < len(top) and j < len(bottom): + if top[i] == bottom[j]: + del top[i] + del bottom[j] + elif top[i] > bottom[j]: + j += 1 + else: + i += 1 + if add_one: + bottom.append(QQ(1)) + else: + try: + i = bottom.index(QQ(1)) + bottom.append(QQ(1)) + del bottom[i] + except ValueError: + bottom.append(QQ(1)) + top.append(QQ(1)) + top.sort() + self.top = tuple(top) + self.bottom = tuple(bottom) + if len(top) == 0 and len(bottom) == 0: + self.d = 1 + self.bound = 1 + else: + self.d = lcm([ a.denominator() for a in top ] + + [ b.denominator() for b in bottom ]) + self.bound = 2 * self.d * max(top + bottom) + 1 + + def __repr__(self): + r""" + Return a string representation of these parameters. + + EXAMPLES:: + + sage: from sage.functions.hypergeometric_algebraic import HypergeometricParameters + sage: pa = HypergeometricParameters([1/4, 1/3, 1/2], [2/5, 3/5]) + sage: pa # indirect doctest + ((1/4, 1/3, 1/2), (2/5, 3/5, 1)) + """ + return "(%s, %s)" % (self.top, self.bottom) + + def __hash__(self): + return hash((self.top, self.bottom)) + + def __eq__(self, other): + return (isinstance(other, HypergeometricParameters) + and self.top == other.top and self.bottom == other.bottom) + + def is_balanced(self): + r""" + Return ``True`` if there are as many top parameters as bottom + parameters; ``False`` otherwise. + + EXAMPLES:: + + sage: from sage.functions.hypergeometric_algebraic import HypergeometricParameters + sage: pa = HypergeometricParameters([1/4, 1/3, 1/2], [2/5, 3/5]) + sage: pa + ((1/4, 1/3, 1/2), (2/5, 3/5, 1)) + sage: pa.is_balanced() + True + + :: + + sage: pa = HypergeometricParameters([1/4, 1/3, 1/2], [2/5, 3/5], add_one=False) + sage: pa + ((1/4, 1/3, 1/2, 1), (2/5, 3/5, 1)) + sage: pa.is_balanced() + False + """ + return len(self.top) == len(self.bottom) + + @cached_method + def christol_sorting(self, c=1): + r""" + Return a sorted list of triples, where each triple is associated to one + of the parameters a, and consists of the decimal part of d*c*a (where + integers are assigned 1 instead of 0), the negative value of a, and a + sign (plus or minus 1), where top parameters are assigned -1 and bottom + parameters +1. Sorting the list lexecographically according to the + first two entries of the tuples sorts the corresponing parameters + according to the total ordering (defined on p.6 in [Chr1986]_). + + INPUT: + + - ``c`` -- an integer (default: ``1``) + + EXAMPLES:: + + sage: from sage.functions.hypergeometric_algebraic import HypergeometricParameters + sage: pa = HypergeometricParameters([1/4, 1/3, 1/2], [2/5, 3/5]) + sage: pa + ((1/4, 1/3, 1/2), (2/5, 3/5, 1)) + sage: pa.christol_sorting(7) + [(12, -3/5, 1), + (20, -1/3, -1), + (30, -1/2, -1), + (45, -1/4, -1), + (48, -2/5, 1), + (60, -1, 1)] + """ + d = self.d + A = [(d - (-d*c*a) % d, -a, -1) for a in self.top] + B = [(d - (-d*c*b) % d, -b, 1) for b in self.bottom] + return sorted(A + B) + + def parenthesis_criterion(self, c): + r""" + Return ``True`` if in each prefix of the list + ``self.christol_sorting(c)`` there are at least as many triples with + third entry -1 as triples with third entry +1. Return ``False`` + otherwise. + + INPUT: + + - ``c`` -- an integer + + EXAMPLES:: + + sage: from sage.functions.hypergeometric_algebraic import HypergeometricParameters + sage: pa = HypergeometricParameters([1/4, 1/3, 1/2], [2/5, 3/5]) + sage: pa + ((1/4, 1/3, 1/2), (2/5, 3/5, 1)) + sage: pa.christol_sorting(7) + [(12, -3/5, 1), + (20, -1/3, -1), + (30, -1/2, -1), + (45, -1/4, -1), + (48, -2/5, 1), + (60, -1, 1)] + sage: pa.parenthesis_criterion(7) + False + sage: pa.christol_sorting(1) + [(15, -1/4, -1), + (20, -1/3, -1), + (24, -2/5, 1), + (30, -1/2, -1), + (36, -3/5, 1), + (60, -1, 1)] + sage: pa.parenthesis_criterion(1) + True + """ + parenthesis = 0 + for _, _, paren in self.christol_sorting(c): + parenthesis += paren + if parenthesis > 0: + return False + return parenthesis <= 0 + + def interlacing_criterion(self, c): + r""" + Return ``True`` if the sorted lists of the decimal parts (where integers + are assigned 1 instead of 0) of c*a and c*b for a in the top parameters + and b in the bottom parameters interlace, i.e., the entries in the sorted + union of the two lists alternate between entries from the first and from + the second list. Used to determine algebraicity of the hypergeometric + function with these parameters with the Beukers-Heckman criterion. + + INPUT: + + - ``c`` -- an integer + + EXAMPLES:: + + sage: from sage.functions.hypergeometric_algebraic import HypergeometricParameters + sage: pa = HypergeometricParameters([1/3, 2/3], [1/2]) + sage: pa + ((1/3, 2/3), (1/2, 1)) + sage: pa.interlacing_criterion(1) + True + sage: pa.interlacing_criterion(5) + True + + :: + + sage: from sage.functions.hypergeometric_algebraic import HypergeometricParameters + sage: pa = HypergeometricParameters([1/8, 3/8, 5/8], [1/4, 1/2]) + sage: pa + ((1/8, 3/8, 5/8), (1/4, 1/2, 1)) + sage: pa.interlacing_criterion(1) + True + sage: pa.interlacing_criterion(3) + False + """ + previous_paren = 1 + for _, _, paren in self.christol_sorting(c): + if paren == previous_paren: + return False + previous_paren = paren + return True + + def q_christol_sorting(self, q): + r""" + Return a sorted list of pairs, one associated to each top parameter a, + and one associated to each bottom parameter b where the pair is either + (1/2 + (-a) % q, -1) or (1 + (-b) % q, 1). + + INPUT: + + - ``q`` -- an integer + + EXAMPLES:: + + sage: from sage.functions.hypergeometric_algebraic import HypergeometricParameters + sage: pa = HypergeometricParameters([1/4, 1/3, 1/2], [2/5, 3/5]) + sage: pa + ((1/4, 1/3, 1/2), (2/5, 3/5, 1)) + sage: pa.q_christol_sorting(7) + [(2, 1), (2.5, -1), (3.5, -1), (5.5, -1), (6, 1), (7, 1)] + """ + A = [(1/2 + (-a) % q, -1) for a in self.top] + B = [(1 + (-b) % q, 1) for b in self.bottom] + return sorted(A + B) + + def q_parenthesis(self, q): + r""" + Return maximal value of the sum of all the second entries of the pairs + in a prefix of ``self.q_christol_sorting(q)`` and the first entry of + the last pair in the prefix of smallest length where this value is + attained. + + INPUT: + + - ``q`` -- an integer + + EXAMPLES:: + + sage: from sage.functions.hypergeometric_algebraic import HypergeometricParameters + sage: pa = HypergeometricParameters([1/4, 1/3, 1/2], [2/5, 3/5]) + sage: pa + ((1/4, 1/3, 1/2), (2/5, 3/5, 1)) + sage: pa.q_christol_sorting(7) + [(2, 1), (2.5, -1), (3.5, -1), (5.5, -1), (6, 1), (7, 1)] + sage: pa.q_parenthesis(7) + (2, 1) + """ + parenthesis = maximum = shift = 0 + for s, paren in self.q_christol_sorting(q): + parenthesis += paren + if parenthesis > maximum: + maximum = parenthesis + shift = s + return shift, maximum + + def q_parenthesis_criterion(self, q): + r""" + Return ``True`` if in each prefix of the list + ``self.q_christol_sorting(q)`` there are at least as many pairs with + second entry -1 as pairs with second entry +1. Return ``False`` + otherwise. + + INPUT: + + - ``q`` -- an integer + + EXAMPLES:: + + sage: from sage.functions.hypergeometric_algebraic import HypergeometricParameters + sage: pa = HypergeometricParameters([1/4, 1/3, 1/2], [2/5, 3/5]) + sage: pa + ((1/4, 1/3, 1/2), (2/5, 3/5, 1)) + sage: pa.q_christol_sorting(7) + [(2, 1), (2.5, -1), (3.5, -1), (5.5, -1), (6, 1), (7, 1)] + sage: pa.q_parenthesis_criterion(7) + False + sage: pa.q_christol_sorting(61) + [(15.5, -1), (20.5, -1), (25, 1), (30.5, -1), (37, 1), (61, 1)] + sage: pa.q_parenthesis_criterion(61) + True + """ + parenthesis = 0 + for _, paren in self.q_christol_sorting(q): + parenthesis += paren + if parenthesis > 0: + return False + return parenthesis <= 0 + + def q_interlacing_number(self, q): + r""" + Return the number of pairs in the list ``self.q_christol_sorting(q)`` + with second entry 1, that were preceded by a pair with second entry + -1. + + INPUT: + + - ``q`` -- an integer. + + EXAMPLES:: + + sage: from sage.functions.hypergeometric_algebraic import HypergeometricParameters + sage: pa = HypergeometricParameters([1/4, 1/3, 1/2], [2/5, 3/5]) + sage: pa + ((1/4, 1/3, 1/2), (2/5, 3/5, 1)) + sage: pa.q_christol_sorting(7) + [(2, 1), (2.5, -1), (3.5, -1), (5.5, -1), (6, 1), (7, 1)] + sage: pa.q_interlacing_number(7) + 1 + """ + interlacing = 0 + previous_paren = 1 + for _, paren in self.q_christol_sorting(q): + if paren == 1 and previous_paren == -1: + interlacing += 1 + previous_paren = paren + return interlacing + + def remove_positive_integer_differences(self): + r""" + Return parameters, where pairs consisting of a top parameter + and a bottom parameter with positive integer differences are + removed, starting with pairs of minimal positive integer + difference. + + EXAMPLES:: + + sage: from sage.functions.hypergeometric_algebraic import HypergeometricParameters + sage: pa = HypergeometricParameters([5/2, -1/2, 5/3], [3/2, 1/3]) + sage: pa + ((-1/2, 5/3, 5/2), (1/3, 3/2, 1)) + sage: pa.remove_positive_integer_differences() + ((-1/2, 5/3), (1/3, 1)) + + The choice of which pair with integer differences to remove first + is important:: + + sage: pa = HypergeometricParameters([4, 2, 1/2], [1, 3]) + sage: pa + ((1/2, 2, 4), (1, 3, 1)) + sage: pa.remove_positive_integer_differences() + ((1/2,), (1,)) + """ + differences = [] + top = list(self.top) + bottom = list(self.bottom) + for i in range(len(top)): + for j in range(len(bottom)): + diff = top[i] - bottom[j] + if diff in ZZ and diff > 0: + differences.append((diff, i, j)) + for _, i, j in sorted(differences): + if top[i] is not None and bottom[j] is not None: + top[i] = None + bottom[j] = None + return HypergeometricParameters(top, bottom, add_one=False) + + def has_negative_integer_differences(self): + r""" + Return ``True`` if there exists a pair of a top parameter and a bottom + parameter, such that the top one minus the bottom one is a negative integer; + return ``False`` otherwise. + + EXAMPLES:: + + sage: from sage.functions.hypergeometric_algebraic import HypergeometricParameters + sage: pa = HypergeometricParameters([1/4, 1/3, 1/2], [2/5, 3/5]) + sage: pa + ((1/4, 1/3, 1/2), (2/5, 3/5, 1)) + sage: pa.has_negative_integer_differences() + False + + :: + + sage: pa = HypergeometricParameters([1/4, 1/3, 1/2], [2/5, 3/2]) + sage: pa + ((1/4, 1/3, 1/2), (2/5, 3/2, 1)) + sage: pa.has_negative_integer_differences() + True + """ + return any(a - b in ZZ and a < b for a in self.top for b in self.bottom) + + def shift(self, s): + r""" + Return the parameters obtained by adding s to each of them. + + INPUT: + + - ``s`` -- a rational number + + EXAMPLES:: + + sage: from sage.functions.hypergeometric_algebraic import HypergeometricParameters + sage: pa = HypergeometricParameters([1/4, 1/3, 1/2], [2/5, 3/5]) + sage: pa + ((1/4, 1/3, 1/2), (2/5, 3/5, 1)) + sage: pa.shift(2) + ((1, 9/4, 7/3, 5/2), (12/5, 13/5, 3, 1)) + """ + top = [a+s for a in self.top] + bottom = [b+s for b in self.bottom] + return HypergeometricParameters(top, bottom, add_one=False) + + def decimal_part(self): + r""" + Return the parameters obtained by taking the decimal part of each of + the parameters, where integers are assigned 1 instead of 0. + + EXAMPLES:: + + sage: from sage.functions.hypergeometric_algebraic import HypergeometricParameters + sage: pa = HypergeometricParameters([5/4, 1/3, 2], [2/5, -2/5]) + sage: pa + ((1/3, 5/4, 2), (-2/5, 2/5, 1)) + sage: pa.decimal_part() + ((1/4, 1/3, 1), (2/5, 3/5, 1)) + """ + top = [1 + a - ceil(a) for a in self.top] + bottom = [1 + b - ceil(b) for b in self.bottom] + return HypergeometricParameters(top, bottom, add_one=False) + + def valuation_position(self, p, drift=0): + r""" + If the `h_k`s are the coefficients of the hypergeometric + series corresponding to these parameters and `\delta` is + the drift, return the smallest value of + + .. MATH:: + + \text{val}_p(h_k) + \delta k + + and the first index `k` where this minimum is reached. + + INPUT: + + - ``p`` -- a prime number + + - ``drift`` -- a rational number (default: ``0``) + + EXAMPLES:: + + sage: from sage.functions.hypergeometric_algebraic import HypergeometricParameters + sage: pa = HypergeometricParameters([1/5, 1/5, 1/5], [1/3, 3^10/5]) + sage: pa.valuation_position(3) + (-9, 1) + + When the relevant sequence is not bounded from below, the + tuple ``(-Infinity, None)`` is returned:: + + sage: pa.valuation_position(5) + (-Infinity, None) + + An example with a drift:: + + sage: pa.valuation_position(3, drift=-7/5) + (-54/5, 7) + """ + # We treat the case of parameters having p in the denominator + # When it happens, we remove the parameter and update the drift accordingly + params = {} + for a in self.top: + if a in ZZ and a <= 0: + raise NotImplementedError # TODO + v = a.valuation(p) + if v < 0: + drift += v + else: + params[a] = params.get(a, 0) + 1 + for b in self.bottom: + v = b.valuation(p) + if v < 0: + drift -= v + else: + params[b] = params.get(b, 0) - 1 + params = [(pa, dw) for pa, dw in params.items() if dw != 0] + + # We check that we are inside the disk of convergence + diff = sum(dw for _, dw in params) + growth = (p-1)*drift + diff + if growth < 0: + return -infinity, None + + # Main part: computation of the valuation + # We use Christol's formula (see Lemma 3.1.10 of [CFVM2025]) + # with modification in order to take the drift into account + n = len(params) + 1 + TSR = TropicalSemiring(QQ) + TM = identity_matrix(TSR, n) + + d = lcm(pa.denominator() for pa, _ in params) + order = IntegerModRing(d)(p).multiplicative_order() + bound = 1 + max(pa.abs() for pa, _ in params) + thresold = d * sum(dw for _, dw in params if dw > 0) + + valuation = position = ZZ(0) + valfinal = signature_prev = None + indices = {} + count = 0 + q = 1 + while True: + # We take into account the contribution of V({k/p^r}, p^r). + # We represent the partial sum until r by the list signature. + # Its entries are triples (valuation, position, parameter): + # - parameter is the parameter corresponding to a point of + # discontinuity of the last summand V({k/p^r}, p^r) + # - valuation is the minimum of the partial sum on the + # range starting at this discontinuity point (included) + # and ending at the next one (excluded) + # - position is the first position where the minimum is reached + # (The dictionary indices allows for finding rapidly + # an entry in signature with a given parameter.) + # The list signature_prev and the dictionary indices_prev + # correspond to the same data for r-1. + + pq = p * q + + # We compute the points of discontinuity of V({k/p^r}, p^r) + # and store them in the list jumps + # Each entry of jumps has the form (x, dw, parameter) where: + # - x is the position of the discontinuity point + # - dw is the jump of V({k/p^r}, p^r) at this point + # - parameter is the corresponding parameter + jumps = [(1 + (-pa) % pq, dw, pa) for pa, dw in params] + jumps = [(0, 0, 0)] + sorted(jumps) + [(pq, None, 0)] + + # We compute the signature + signature = [] + indices = {} + w = 0 + TMstep = matrix(TSR, n) + for i in range(n): + x, dw, param = jumps[i] # discontinuity point + y, _, right = jumps[i+1] # next discontinuity point + w += dw # the value of V({k/p^r}, p^r) on this interval + indices[param] = len(signature) + if x == y: + # Case of empty interval + val = infinity + pos = None + elif signature_prev is None: + # Case r = 1 + if drift < 0: + pos = y - 1 + else: + pos = x + val = drift * pos + else: + # Case r > 1 + # The variable complete stores whether the interval + # [x,y] covers all [0, p^(r-1)) modulo p^(r-1) + complete = (y - x >= q) + if complete and drift < 0: + interval = ((y-1) // q) - 1 + j = j0 = indices_prev[right] + else: + interval = max(0, (x-1) // q) + j = j0 = indices_prev[param] + val = infinity + while True: + valj, posj, paramj = signature_prev[j] + valj += drift * interval + TMstep[i,j] = TSR(drift*interval + w) + if valj < val: + val = valj + pos = posj + q * interval + j += 1 + if j >= n: + j = 0 + interval += 1 + if j == j0 or (not complete and signature_prev[j][2] == right): + break + signature.append((val + w, pos, param)) + + # The halting criterion + minimum = min(signature) + valuation, position, _ = minimum + if q > bound: + if drift > 2*thresold and all(signature[i][0] > valuation + thresold for i in range(1, n)): + return valuation, position + if growth == 0: + if count < order: + TM = TMstep * TM + count += 1 + if count == order: + try: + TM = TM.weak_transitive_closure() + except ValueError: + return -infinity, None + valfinal = min(TM[i,j].lift() + signature[j][0] + for i in range(n) for j in range(n)) + if valuation == valfinal: + return valuation, position + + # We update the values for the next r + q = pq + drift = p*drift + diff + signature_prev = signature + indices_prev = indices + + def dwork_image(self, p): + r""" + Return the parameters obtained by applying the Dwork map to each of + the parameters. The Dwork map `D_p(x)` of a p-adic integer x is defined + as the unique p-adic integer such that `p D_p(x) - x` is a nonnegative + integer smaller than p. + + INPUT: + + - ``p`` -- a prime number + + EXAMPLE:: + + sage: from sage.functions.hypergeometric_algebraic import HypergeometricParameters + sage: pa = HypergeometricParameters([1/4, 1/3, 1/2], [2/5, 3/5]) + sage: pa + ((1/4, 1/3, 1/2), (2/5, 3/5, 1)) + sage: pa.dwork_image(7) + ((1/3, 1/2, 3/4), (1/5, 4/5, 1)) + + If `p` is not coprime to the common denominators of the parameters, + a ``ValueError`` is raised:: + + sage: pa.dwork_image(3) + Traceback (most recent call last): + ... + ValueError: denominators of parameters are not coprime to p + """ + try: + top = [(a + (-a) % p) / p for a in self.top] + bottom = [(b + (-b) % p) / p for b in self.bottom] + except ZeroDivisionError: + raise ValueError("denominators of parameters are not coprime to p") + return HypergeometricParameters(top, bottom, add_one=False) + + def frobenius_order(self, p): + r""" + Return the Frobenius order of the hypergeometric function with this set + of parameters, that is the order of the Dwork map acting on the decimal + parts of the parameters. + + INPUT: + + - ``p`` -- a prime number + + EXAMPLES:: + + sage: from sage.functions.hypergeometric_algebraic import HypergeometricParameters + sage: pa = HypergeometricParameters([1/4, 1/3, 1/2], [2/5, 3/5]) + sage: pa + ((1/4, 1/3, 1/2), (2/5, 3/5, 1)) + sage: pa.frobenius_order(7) + 2 + """ + param = self.decimal_part() + iter = param.dwork_image(p) + i = 1 + while param != iter: + iter = iter.dwork_image(p) + i += 1 + return i diff --git a/src/sage/functions/meson.build b/src/sage/functions/meson.build index c37ec96e9ff..4b781dadfcf 100644 --- a/src/sage/functions/meson.build +++ b/src/sage/functions/meson.build @@ -9,6 +9,8 @@ py.install_sources( 'generalized.py', 'hyperbolic.py', 'hypergeometric.py', + 'hypergeometric_algebraic.py', + 'hypergeometric_parameters.py', 'jacobi.py', 'log.py', 'min_max.py', diff --git a/src/sage/matrix/action.pyx b/src/sage/matrix/action.pyx index 732a1312ca1..f40cf212066 100644 --- a/src/sage/matrix/action.pyx +++ b/src/sage/matrix/action.pyx @@ -279,7 +279,7 @@ cdef class MatrixMatrixAction(MatrixMulAction): B = B.dense_matrix() else: A = A.dense_matrix() - assert type(A) is type(B), (type(A), type(B)) + # assert type(A) is type(B), (type(A), type(B)) prod = A._matrix_times_matrix_(B) if A._subdivisions is not None or B._subdivisions is not None: Asubs = A.subdivisions() diff --git a/src/sage/matrix/matrix_space.py b/src/sage/matrix/matrix_space.py index eafab99b0e9..ed27f7a41bb 100644 --- a/src/sage/matrix/matrix_space.py +++ b/src/sage/matrix/matrix_space.py @@ -50,6 +50,8 @@ from sage.misc.lazy_attribute import lazy_attribute from sage.misc.superseded import deprecated_function_alias from sage.misc.persist import register_unpickle_override +from sage.categories.sets_cat import Sets +from sage.categories.semirings import Semirings from sage.categories.rings import Rings from sage.categories.fields import Fields from sage.categories.enumerated_sets import EnumeratedSets @@ -60,7 +62,6 @@ feature=Meataxe()) lazy_import('sage.groups.matrix_gps.matrix_group', ['MatrixGroup_base']) -_Rings = Rings() _Fields = Fields() @@ -320,6 +321,15 @@ def get_matrix_class(R, nrows, ncols, sparse, implementation): else: return matrix_laurent_mpolynomial_dense.Matrix_laurent_mpolynomial_dense + try: + from sage.rings.semirings.tropical_semiring import TropicalSemiring + except ImportError: + pass + else: + if isinstance(R, TropicalSemiring): + from sage.rings.semirings import tropical_matrix + return tropical_matrix.Matrix_tropical_dense + # The fallback from sage.matrix.matrix_generic_dense import Matrix_generic_dense return Matrix_generic_dense @@ -725,8 +735,8 @@ def __classcall__(cls, base_ring, sage: MS2._my_option False """ - if base_ring not in _Rings: - raise TypeError("base_ring (=%s) must be a ring" % base_ring) + if base_ring not in Semirings(): + raise TypeError("base_ring (=%s) must be a ring or a semiring" % base_ring) if ncols_or_column_keys is not None: try: @@ -898,12 +908,17 @@ def __init__(self, base_ring, nrows, ncols, sparse, implementation): from sage.categories.modules import Modules from sage.categories.algebras import Algebras - if nrows == ncols: - category = Algebras(base_ring.category()) + if base_ring in Rings(): + if nrows == ncols: + category = Algebras(base_ring.category()) + else: + category = Modules(base_ring.category()) + category = category.WithBasis().FiniteDimensional() else: - category = Modules(base_ring.category()) - - category = category.WithBasis().FiniteDimensional() + if nrows == ncols: + category = Semirings() + else: + category = Sets() if not self.__nrows or not self.__ncols: is_finite = True @@ -2024,8 +2039,9 @@ def identity_matrix(self): if self.__nrows != self.__ncols: raise TypeError("identity matrix must be square") A = self.zero_matrix().__copy__() + one = self.base_ring().one() for i in range(self.__nrows): - A[i, i] = 1 + A[i, i] = one A.set_immutable() return A diff --git a/src/sage/matrix/special.py b/src/sage/matrix/special.py index 7182c27662f..f5eb173479c 100644 --- a/src/sage/matrix/special.py +++ b/src/sage/matrix/special.py @@ -936,11 +936,19 @@ def identity_matrix(ring, n=0, sparse=False): Full MatrixSpace of 3 by 3 sparse matrices over Integer Ring sage: M.is_mutable() True + + :: + + sage: T = TropicalSemiring(QQ) + sage: identity_matrix(T, 3) + [ 0 +infinity +infinity] + [+infinity 0 +infinity] + [+infinity +infinity 0] """ if isinstance(ring, (Integer, int)): n = ring ring = ZZ - return matrix_space.MatrixSpace(ring, n, n, sparse)(1) + return matrix_space.MatrixSpace(ring, n, n, sparse)(ring.one()) @matrix_method diff --git a/src/sage/misc/sagedoc.py b/src/sage/misc/sagedoc.py index ab5ff639dd6..7ff39a8a7ee 100644 --- a/src/sage/misc/sagedoc.py +++ b/src/sage/misc/sagedoc.py @@ -1454,7 +1454,7 @@ class _sage_doc: sage: browse_sage_doc._open("reference", testing=True)[0] # needs sagemath_doc_html 'http://localhost:8000/doc/live/reference/index.html' - sage: browse_sage_doc(identity_matrix, 'rst')[-107:-47] # needs sage.modules + sage: browse_sage_doc(identity_matrix, 'rst')[-311:-251] # needs sage.modules '...Full MatrixSpace of 3 by 3 sparse matrices...' """ def __init__(self): diff --git a/src/sage/rings/semirings/meson.build b/src/sage/rings/semirings/meson.build index b92ea837850..0a9762b756d 100644 --- a/src/sage/rings/semirings/meson.build +++ b/src/sage/rings/semirings/meson.build @@ -2,6 +2,7 @@ py.install_sources( '__init__.py', 'all.py', 'non_negative_integer_semiring.py', + 'tropical_matrix.py', 'tropical_mpolynomial.py', 'tropical_polynomial.py', 'tropical_semiring.pyx', diff --git a/src/sage/rings/semirings/tropical_matrix.py b/src/sage/rings/semirings/tropical_matrix.py new file mode 100644 index 00000000000..12761978673 --- /dev/null +++ b/src/sage/rings/semirings/tropical_matrix.py @@ -0,0 +1,199 @@ +r""" +Matrices over tropical semirings + +AUTHORS: + +- Xavier Caruso (2025-11): initial version +""" + +# **************************************************************************** +# Copyright (C) 2025 Xavier Caruso +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 2 of the License, or +# (at your option) any later version. +# https://www.gnu.org/licenses/ +# **************************************************************************** + +from sage.matrix.constructor import matrix +from sage.matrix.matrix_generic_dense import Matrix_generic_dense +from sage.rings.infinity import infinity + + +class Matrix_tropical_dense(Matrix_generic_dense): + r""" + A class for dense matrices over a tropical semiring. + + EXAMPLES:: + + sage: from sage.rings.semirings.tropical_matrix import Matrix_tropical_dense + sage: T = TropicalSemiring(QQ) + sage: M = matrix(T, [[1, 2], [3, 4]]) + sage: isinstance(M, Matrix_tropical_dense) + True + """ + def extremum_cycle_mean(self): + r""" + Return the extremal (that is, minimal if the addition is max + and maximum is the addition is min) mean weight of this matrix + It is also the smallest/largest eigenvalue of this matrix. + + ALGORITHM: + + We implement Karp's algorithm described in []_, Section 1.6.1. + + EXAMPLES:: + + sage: T = TropicalSemiring(QQ, use_min=False) + sage: M = matrix(T, [[-2, 1, -3], + ....: [ 3, 0, 3], + ....: [ 5, 2, 1]]) + sage: M.extremum_cycle_mean() + 3 + + :: + + sage: T = TropicalSemiring(QQ) + sage: z = T.zero() + sage: M = matrix(T, [[z, 1, 10, z], + ....: [z, z, 3, z], + ....: [z, z, z, 2], + ....: [8, 0, z, z]]) + sage: M.extremum_cycle_mean() + 5/3 + """ + T = self.base_ring() + n = self.ncols() + if self.nrows() != n: + raise TypeError("matrix must be square") + v = matrix(1, n, n*[T.one()]) + vs = [v] + for _ in range(n): + v = v * self + vs.append(v) + w = [vs[n][0,j].lift() for j in range(n)] + if T._use_min: + return min(max((w[j] - vs[k][0,j].lift()) / (n-k) for k in range(n)) + for j in range(n) if w[j] is not infinity) + else: + return max(min((w[j] - vs[k][0,j].lift()) / (n-k) for k in range(n)) + for j in range(n) if w[j] is not infinity) + + def weak_transitive_closure(self): + r""" + Return the weak transitive closure of this matrix `M`, + that is, by definition + + .. MATH:: + + A \oplus A^2 \oplus A^3 \oplus A^4 \oplus \cdots + + or raise an error if this sum does not converge. + + ALGORITHM: + + We implement the Floyd-Warshall algorithm described in + [But2010]_, Algorithm 1.6.21. + + EXAMPLES:: + + sage: T = TropicalSemiring(QQ) + sage: z = T.zero() + sage: M = matrix(T, [[z, 1, 10, z], + ....: [z, z, 3, z], + ....: [z, z, z, 2], + ....: [8, 0, z, z]]) + sage: M.weak_transitive_closure() + [14 1 4 6] + [13 5 3 5] + [10 2 5 2] + [ 8 0 3 5] + + We check that the minimal cycle mean of `M` is the largest + value `a` such that `(-a) \otimes M` has a weak transitive + closure:: + + sage: M.extremum_cycle_mean() + 5/3 + sage: aM = T(-5/3) * M + sage: aM.weak_transitive_closure() + [22/3 -2/3 2/3 1] + [ 8 0 4/3 5/3] + [20/3 -4/3 0 1/3] + [19/3 -5/3 -1/3 0] + sage: bM = T(-2) * M + sage: bM.weak_transitive_closure() + Traceback (most recent call last): + ... + ValueError: negative cycle exists + + .. SEEALSO:: + + :meth:`strong_transitive_closure` + """ + T = self.base_ring() + n = self.ncols() + if self.nrows() != n: + raise TypeError("matrix must be square") + G = self.__copy__() + for p in range(n): + for i in range(n): + if i == p: + continue + for j in range(n): + if j == p: + continue + G[i,j] += G[i,p] * G[p,j] + if i == j: + if T._use_min and G[i,i].lift() < 0: + raise ValueError("negative cycle exists") + if not T._use_min and G[i,i].lift() > 0: + raise ValueError("positive cycle exists") + return G + + def strong_transitive_closure(self): + r""" + Return the string transitive closure of this matrix `M`, + that is, by definition + + .. MATH:: + + I \oplus A \oplus A^2 \oplus A^3 \oplus A^4 \oplus \cdots + + or raise an error if this sum does not converge. + + ALGORITHM: + + We implement the Floyd-Warshall algorithm described in + [But2010]_, Algorithm 1.6.21. + + EXAMPLES:: + + sage: T = TropicalSemiring(QQ, use_min=False) + sage: M = matrix(T, [[-5, -2, -6], + ....: [ 0, -3, 0], + ....: [ 2, -1, -2]]) + sage: M.strong_transitive_closure() + [ 0 -2 -2] + [ 2 0 0] + [ 2 0 0] + + :: + + sage: T = TropicalSemiring(QQ) + sage: M = matrix(T, [[-5, -2, -6], + ....: [ 0, -3, 0], + ....: [ 2, -1, -2]]) + sage: M.strong_transitive_closure() + Traceback (most recent call last): + ... + ValueError: negative cycle exists + + .. SEEALSO:: + + :meth:`weak_transitive_closure` + """ + return self.parent().identity_matrix() + self.weak_transitive_closure() + + kleene_star = strong_transitive_closure diff --git a/src/sage/rings/tate_algebra_element.pyx b/src/sage/rings/tate_algebra_element.pyx index ade9ea418ad..c109c2a967c 100644 --- a/src/sage/rings/tate_algebra_element.pyx +++ b/src/sage/rings/tate_algebra_element.pyx @@ -1123,7 +1123,7 @@ cdef class TateAlgebraElement(CommutativeAlgebraElement): sage: S. = TateAlgebra(Qp(5), log_radii=(1,0)) sage: f = 5*x sage: f.add_bigoh(1) - (5 + O(5^2))*x + O(5 * ) + (5 + O(5^2))*x + O(5 * <5*x, y>) """ self._is_normalized = True if self._prec is Infinity: @@ -1169,6 +1169,15 @@ cdef class TateAlgebraElement(CommutativeAlgebraElement): sage: A(x + 2*x^2 + x^3, prec=5) ...00001*x^3 + ...00001*x + ...00010*x^2 + O(2^5 * ) + + TESTS:: + + sage: S. = TateAlgebra(R, log_radii=[-1]) + sage: S(x, 5) + ...0001*x + O(2^5 * ) + sage: S. = TateAlgebra(R, log_radii=[1]) + sage: S(x, 5) + ...000001*x + O(2^5 * <2*x>) """ vars = self._parent.variable_names() s = "" @@ -1191,14 +1200,14 @@ cdef class TateAlgebraElement(CommutativeAlgebraElement): for i in range(len(vars)): if lr[i] == 0: sv.append(vars[i]) - elif lr[i] == -1: - sv.append("%s*%s" % (su, vars[i])) elif lr[i] == 1: + sv.append("%s*%s" % (su, vars[i])) + elif lr[i] == -1: sv.append("%s/%s" % (vars[i], su)) - elif lr[i] < 0: - sv.append("%s^%s*%s" % (su, -lr[i], vars[i])) + elif lr[i] > 0: + sv.append("%s^%s*%s" % (su, lr[i], vars[i])) else: - sv.append("%s/%s^%s" % (vars[i], su, lr[i])) + sv.append("%s/%s^%s" % (vars[i], su, -lr[i])) sv = ", ".join(sv) if self._prec == 0: s += "O(<%s>)" % sv @@ -2534,8 +2543,8 @@ cdef class TateAlgebraElement(CommutativeAlgebraElement): However `\log(1+x)` converges on a smaller disk:: sage: f.restriction(-1).log() - ...000000001*x + ...0000000.1*x^3 + ...111111*x^2 + ... - + O(3^10 * <3*x, 3*y>) + ...000000001*x + ...0000000.1*x^3 + ...11111111*x^2 + ... + + O(3^10 * ) TESTS:: @@ -2692,8 +2701,8 @@ cdef class TateAlgebraElement(CommutativeAlgebraElement): However `\exp(x)` converges on a smaller disk:: sage: f.restriction(-1).exp() - ...0000000001 + ...000000001*x + ...1111111.2*x^3 + ...111112*x^2 - + ... + O(3^10 * <3*x, 3*y>) + ...0000000001 + ...000000001*x + ...1111111.2*x^3 + ...11111112*x^2 + + ... + O(3^10 * ) TESTS::