Source code for flatsurf.geometry.subfield

r"""
Various helper functions for computations in number fields.

References

.. [KoRoTr2015]
    Pierre-Vincent Koseleff, Fabrice Rouillier, Cuong Tran.
    On the Sign of a Trigonometric Expression.
    ISSAC 2015. Proceedings of the 2015 ACM International
    Symposium on Symbolic and Algebraic Computation.
    Pages 259--266.
    DOI: http://dx.doi.org/10.1145/2755996.2756664
"""
######################################################################
#  This file is part of sage-flatsurf.
#
#        Copyright (C) 2020 Samuel Lelièvre
#                      2020 Vincent Delecroix
#                      2023 Julian Rüth
#
#  sage-flatsurf 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.
#
#  sage-flatsurf is distributed in the hope that it will be useful,
#  but WITHOUT ANY WARRANTY; without even the implied warranty of
#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#  GNU General Public License for more details.
#
#  You should have received a copy of the GNU General Public License
#  along with sage-flatsurf. If not, see <https://www.gnu.org/licenses/>.
######################################################################
from sage.misc.misc_c import prod
from sage.rings.all import ZZ, QQ, AA, NumberField, polygen
from sage.structure.element import parent
from sage.modules.free_module import VectorSpace
from sage.matrix.constructor import matrix
from sage.categories.homset import Hom
from sage.categories.fields import Fields
from sage.rings.qqbar import do_polred


[docs] def number_field_elements_from_algebraics(elts, name="a"): r""" The native Sage function ``number_field_elements_from_algebraics`` currently returns number field *without* embedding. This function return field with embedding! EXAMPLES:: sage: from flatsurf.geometry.subfield import number_field_elements_from_algebraics sage: z = QQbar.zeta(5) sage: c = z.real() sage: s = z.imag() sage: number_field_elements_from_algebraics((c, s)) (Number Field in a with defining polynomial y^4 - 5*y^2 + 5 with a = 1.902113032590308?, [1/2*a^2 - 3/2, 1/2*a]) sage: number_field_elements_from_algebraics([AA(1), AA(2/3)]) (Rational Field, [1, 2/3]) """ # case when all elements are rationals if all(x in QQ for x in elts): return QQ, [QQ(x) for x in elts] # general case from sage.rings.qqbar import number_field_elements_from_algebraics field, elts, phi = number_field_elements_from_algebraics(elts, minimal=True) polys = [x.polynomial() for x in elts] K = NumberField(field.polynomial(), name, embedding=AA(phi(field.gen()))) gen = K.gen() return K, [x.polynomial()(gen) for x in elts]
[docs] def subfield_from_elements(self, alpha, name=None, polred=True, threshold=None): r""" Return the subfield generated by the elements ``alpha``. .. NOTE:: See https://trac.sagemath.org/ticket/29331 INPUT: - ``alpha`` - list of elements in this number field - ``name`` - a name for the generator of the new number field - ``polred`` (boolean, default ``True``) - whether to optimize the generator of the newly created field - ``threshold`` (positive number, default ``None``) - threshold to be passed to the ``do_polred`` function EXAMPLES:: sage: from flatsurf.geometry.subfield import subfield_from_elements sage: x = polygen(QQ) sage: poly = x^4 - 4*x^2 + 1 sage: emb = AA.polynomial_root(poly, RIF(0.51, 0.52)) sage: K.<a> = NumberField(poly, embedding=emb) sage: sqrt2 = -a^3 + 3*a sage: sqrt3 = -a^2 + 2 sage: assert sqrt2 ** 2 == 2 and sqrt3 ** 2 == 3 sage: L, elts, phi = subfield_from_elements(K, [sqrt2, 1 - sqrt2/3]) sage: L Number Field in a0 with defining polynomial x^2 - 2 with a0 = 1.414213562373095? sage: elts [a0, -1/3*a0 + 1] sage: phi Ring morphism: From: Number Field in a0 with defining polynomial x^2 - 2 with a0 = 1.414213562373095? To: Number Field in a with defining polynomial x^4 - 4*x^2 + 1 with a = 0.5176380902050415? Defn: a0 |--> -a^3 + 3*a sage: assert phi(elts[0]) == sqrt2 sage: assert phi(elts[1]) == 1 - sqrt2/3 sage: L, elts, phi = subfield_from_elements(K, [1, sqrt3]) sage: assert phi(elts[0]) == 1 sage: assert phi(elts[1]) == sqrt3 TESTS:: sage: from flatsurf.geometry.subfield import subfield_from_elements sage: x = polygen(QQ) sage: p = x^8 - 12*x^6 + 23*x^4 - 12*x^2 + 1 sage: K.<a> = NumberField(p) sage: sqrt2 = 6/7*a^7 - 71/7*a^5 + 125/7*a^3 - 43/7*a sage: sqrt3 = 3/7*a^6 - 32/7*a^4 + 24/7*a^2 + 10/7 sage: sqrt5 = 8/7*a^6 - 90/7*a^4 + 120/7*a^2 - 27/7 sage: assert sqrt2**2 == 2 and sqrt3**2 == 3 and sqrt5**2 == 5 sage: L, elts, phi = subfield_from_elements(K, [sqrt2, sqrt3], name='b') sage: assert phi(elts[0]) == sqrt2 sage: assert phi(elts[1]) == sqrt3 sage: L, elts, phi = subfield_from_elements(K, [sqrt2, sqrt5], name='b') sage: assert phi(elts[0]) == sqrt2 sage: assert phi(elts[1]) == sqrt5 sage: L, elts, phi = subfield_from_elements(K, [sqrt3, sqrt5], name='b') sage: assert phi(elts[0]) == sqrt3 sage: assert phi(elts[1]) == sqrt5 sage: L, elts, phi = subfield_from_elements(K, [-149582/214245 + 21423/5581*sqrt2], name='b') sage: assert L.polynomial() == x^2 - 2 sage: L, elts, phi = subfield_from_elements(K, [131490/777 - 1359/22*sqrt3], name='b') sage: assert L.polynomial() == x^2 - 3 sage: L, elts, phi = subfield_from_elements(K, [12241829/177 - 321121/22459 * sqrt5], name='b') sage: assert L.polynomial() == x^2 - x - 1 sage: from sage.rings.qqbar import number_field_elements_from_algebraics sage: R.<x> = QQ[] sage: p1 = x^3 - x - 1 sage: roots1 = p1.roots(QQbar, False) sage: for _ in range(10): # long time (1.5s) ....: p2 = R.random_element(degree=2) ....: while not p2.is_irreducible(): p2 = R.random_element(degree=2) ....: roots2 = p2.roots(QQbar, False) ....: K, (a1,b1,c1,a2,b2), phi = number_field_elements_from_algebraics(roots1 + roots2) ....: u1 = 1 - a1/17 + 3/7*a1**2 ....: u2 = 2 + 33/35 * a1 ....: L, (v1,v2), phi = subfield_from_elements(K, [u1, u2], threshold=100) ....: assert L.polynomial() == p1 ....: assert phi(v1) == u1 and phi(v2) == u2 """ V = VectorSpace(QQ, self.degree()) alpha = [self(a) for a in alpha] # Rational case (degree 0) if all(a in QQ for a in alpha): return (QQ, [QQ(a) for a in alpha], self.coerce_map_from(QQ)) # Trivial maximal case (an element generating the field) if any(a.minpoly().degree() == self.degree() for a in alpha): return (self, alpha, Hom(self, self, Fields()).identity()) # Saturate with multiplication vecs = [(a * a.denominator()).vector() for a in alpha] U = V.subspace(vecs) modified = True while modified: modified = False d = U.dimension() if d == self.degree(): return (self, alpha, Hom(self, self, Fields()).identity()) B = U.basis() new_vecs = [ (self(B[i]) * self(B[j])).vector() for i in range(d) for j in range(i, d) ] if any(vv not in U for vv in new_vecs): U = V.subspace(list(B) + new_vecs) modified = True # Strict subfield, find a generator vgen = None for b in U.basis(): if self(b).minpoly().degree() == d: vgen = b break if vgen is None: s = 1 while True: vgen = U.random_element(proba=0.5, x=-s, y=s) if self(vgen).minpoly().degree() == d: break s *= 2 # Minimize the generator via PARI polred gen = self(vgen) p = gen.minpoly() if polred: if threshold: fwd, back, q = do_polred(p, threshold) else: fwd, back, q = do_polred(p) else: q = p fwd = back = self.polynomial_ring().gen() new_gen = fwd(gen) assert new_gen.minpoly() == q K, hom = self.subfield(new_gen, name=name) # need to express the elements in the basis 1, a, a^2, ..., a^(d-1) M = matrix(QQ, [(new_gen**i).vector() for i in range(d)]) new_alpha = [K(M.solve_left(elt.vector())) for elt in alpha] return (K, new_alpha, hom)
[docs] def is_embedded_subfield(K, L, certificate=False): r""" Return whether there exists a field morphism from ``K`` to ``L`` compatible with the embeddings. EXAMPLES:: sage: from flatsurf.geometry.subfield import is_embedded_subfield sage: x = polygen(QQ) sage: y = polygen(QQ, 'y') sage: z = polygen(QQ, 'z') sage: p = x^4 - 4*x^2 + 1 sage: emb = AA.polynomial_root(p, RIF(0.5, 0.6)) sage: L = NumberField(p, 'a', embedding=emb) sage: K1 = NumberField(y^2 - 2, 'sqrt2', embedding=AA(2)**(1/2)) sage: K2 = NumberField(z^2 - 2, 'msqrt2', embedding=-AA(2)**(1/2)) sage: is_embedded_subfield(K1, L) True sage: is_embedded_subfield(K1, L, True) (True, Generic morphism: From: Number Field in sqrt2 with defining polynomial y^2 - 2 with sqrt2 = 1.414213562373095? To: Number Field in a with defining polynomial x^4 - 4*x^2 + 1 with a = 0.5176380902050415? Defn: sqrt2 -> -a^3 + 3*a) sage: is_embedded_subfield(K2, L, True) (True, Generic morphism: From: Number Field in msqrt2 with defining polynomial z^2 - 2 with msqrt2 = -1.414213562373095? To: Number Field in a with defining polynomial x^4 - 4*x^2 + 1 with a = 0.5176380902050415? Defn: msqrt2 -> a^3 - 3*a) sage: is_embedded_subfield(K1, K2, True) (True, Generic morphism: From: Number Field in sqrt2 with defining polynomial y^2 - 2 with sqrt2 = 1.414213562373095? To: Number Field in msqrt2 with defining polynomial z^2 - 2 with msqrt2 = -1.414213562373095? Defn: sqrt2 -> -msqrt2) """ assert K.coerce_embedding() is not None assert L.coerce_embedding() is not None assert K.coerce_embedding().codomain() == L.coerce_embedding().codomain() pK = K.polynomial() pL = L.polynomial() emb = K.coerce_embedding().gen_image() P = pL.parent() pK = P(pK) phi = L.coerce_embedding() for r in pK.roots(L, multiplicities=False): if phi(r) == emb: return (True, K.hom(L, [r])) if certificate else True return (False, None) if certificate else False
[docs] def chebyshev_T(n, c): r""" Return the Chebyshev polynomial T_n so that T_n(2 cos(x)) = 2 cos(n x) .. NOTE:: We use the polynomial as defined in [KoRoTr2015]_, not the Sage version which uses a different normalization. EXAMPLES:: sage: from flatsurf.geometry.subfield import chebyshev_T, cos_minpoly sage: x = polygen(ZZ) sage: chebyshev_T(1, x) - 1 == cos_minpoly(3, x) True sage: chebyshev_T(2, x) - chebyshev_T(1, x) + 1 == cos_minpoly(5, x) True sage: chebyshev_T(3, x) - chebyshev_T(2, x) + chebyshev_T(1, x) - 1 == cos_minpoly(7, x) True sage: cos_minpoly(5, x)(chebyshev_T(3, x)) == cos_minpoly(15, x) * cos_minpoly(5, x) True sage: cos_minpoly(3, x)(chebyshev_T(5, x)) == cos_minpoly(15, x) * cos_minpoly(3, x) True """ # T_0(x) = 2 # T_1(x) = x # and T_{n+1}(x) = x T_n(x) - T_{n-1}(x) T0 = parent(c)(2) if n == 0: return T0 T1 = c for i in range(n - 1): T0, T1 = T1, c * T1 - T0 return T1
[docs] def cos_minpoly_odd_prime(p, var): r""" (-1)^k + (-1)^{k-1} T1 + ... + T_k EXAMPLES:: sage: from flatsurf.geometry.subfield import cos_minpoly_odd_prime sage: x = polygen(ZZ) sage: cos_minpoly_odd_prime(3, x) x - 1 sage: cos_minpoly_odd_prime(5, x) x^2 - x - 1 sage: cos_minpoly_odd_prime(7, x) x^3 - x^2 - 2*x + 1 sage: cos_minpoly_odd_prime(11, x) x^5 - x^4 - 4*x^3 + 3*x^2 + 3*x - 1 """ T0 = 2 T1 = var k = (p - 1) // 2 s = (-1) ** k minpoly = s * (1 - T1) for i in range(k - 1): T0, T1 = T1, var * T1 - T0 minpoly += s * T1 s *= -1 return minpoly
[docs] def cos_minpoly(n, var="x"): r""" Return the minimal polynomial of 2 cos pi/n Done via [KoRoTr2015]_ Algorithm 3. EXAMPLES:: sage: from flatsurf.geometry.subfield import cos_minpoly sage: cos_minpoly(1) x + 2 sage: cos_minpoly(2) x sage: cos_minpoly(3) x - 1 sage: cos_minpoly(4) x^2 - 2 sage: cos_minpoly(5) x^2 - x - 1 sage: cos_minpoly(6) x^2 - 3 sage: cos_minpoly(8) x^4 - 4*x^2 + 2 sage: cos_minpoly(9) x^3 - 3*x - 1 sage: cos_minpoly(10) x^4 - 5*x^2 + 5 sage: cos_minpoly(90) x^24 - 24*x^22 + 252*x^20 - 1519*x^18 + 5796*x^16 - 14553*x^14 + 24206*x^12 - 26169*x^10 + 17523*x^8 - 6623*x^6 + 1182*x^4 - 72*x^2 + 1 sage: cos_minpoly(148) x^72 - 72*x^70 + 2483*x^68 - 54604*x^66 + 860081*x^64 ... - 3511656*x^6 + 77691*x^4 - 684*x^2 + 1 """ n = ZZ(n) facs = list(n.factor()) if isinstance(var, str): var = polygen(ZZ, var) if not facs: # 2 cos (pi) = -2 return var + 2 if facs[0][0] == 2: # n is even k = facs[0][1] facs.pop(0) else: k = 0 if not facs: # 0. n = 2^k # ([KoRoTr2015] Lemma 12) return chebyshev_T(2 ** (k - 1), var) # 1. Compute M_{n0} = M_{p1 ... ps} # ([KoRoTr2015] Lemma 14 and Lemma 15) M = cos_minpoly_odd_prime(facs[0][0], var) for i in range(1, len(facs)): p = facs[i][0] M, r = M(chebyshev_T(p, var)).quo_rem(M) assert r.is_zero() # 2. Compute M_{2^k p1^{a1} ... ps^{as}} # ([KoRoTr2015] Lemma 12) nn = 2**k * prod(p ** (a - 1) for p, a in facs) if nn != 1: M = M(chebyshev_T(nn, var)) return M