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