Coverage for /usr/lib/python3/dist-packages/sympy/polys/numberfields/subfield.py: 14%
174 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1r"""
2Functions in ``polys.numberfields.subfield`` solve the "Subfield Problem" and
3allied problems, for algebraic number fields.
5Following Cohen (see [Cohen93]_ Section 4.5), we can define the main problem as
6follows:
8* **Subfield Problem:**
10 Given two number fields $\mathbb{Q}(\alpha)$, $\mathbb{Q}(\beta)$
11 via the minimal polynomials for their generators $\alpha$ and $\beta$, decide
12 whether one field is isomorphic to a subfield of the other.
14From a solution to this problem flow solutions to the following problems as
15well:
17* **Primitive Element Problem:**
19 Given several algebraic numbers
20 $\alpha_1, \ldots, \alpha_m$, compute a single algebraic number $\theta$
21 such that $\mathbb{Q}(\alpha_1, \ldots, \alpha_m) = \mathbb{Q}(\theta)$.
23* **Field Isomorphism Problem:**
25 Decide whether two number fields
26 $\mathbb{Q}(\alpha)$, $\mathbb{Q}(\beta)$ are isomorphic.
28* **Field Membership Problem:**
30 Given two algebraic numbers $\alpha$,
31 $\beta$, decide whether $\alpha \in \mathbb{Q}(\beta)$, and if so write
32 $\alpha = f(\beta)$ for some $f(x) \in \mathbb{Q}[x]$.
33"""
35from sympy.core.add import Add
36from sympy.core.numbers import AlgebraicNumber
37from sympy.core.singleton import S
38from sympy.core.symbol import Dummy
39from sympy.core.sympify import sympify, _sympify
40from sympy.ntheory import sieve
41from sympy.polys.densetools import dup_eval
42from sympy.polys.domains import QQ
43from sympy.polys.numberfields.minpoly import _choose_factor, minimal_polynomial
44from sympy.polys.polyerrors import IsomorphismFailed
45from sympy.polys.polytools import Poly, PurePoly, factor_list
46from sympy.utilities import public
48from mpmath import MPContext
51def is_isomorphism_possible(a, b):
52 """Necessary but not sufficient test for isomorphism. """
53 n = a.minpoly.degree()
54 m = b.minpoly.degree()
56 if m % n != 0:
57 return False
59 if n == m:
60 return True
62 da = a.minpoly.discriminant()
63 db = b.minpoly.discriminant()
65 i, k, half = 1, m//n, db//2
67 while True:
68 p = sieve[i]
69 P = p**k
71 if P > half:
72 break
74 if ((da % p) % 2) and not (db % P):
75 return False
77 i += 1
79 return True
82def field_isomorphism_pslq(a, b):
83 """Construct field isomorphism using PSLQ algorithm. """
84 if not a.root.is_real or not b.root.is_real:
85 raise NotImplementedError("PSLQ doesn't support complex coefficients")
87 f = a.minpoly
88 g = b.minpoly.replace(f.gen)
90 n, m, prev = 100, b.minpoly.degree(), None
91 ctx = MPContext()
93 for i in range(1, 5):
94 A = a.root.evalf(n)
95 B = b.root.evalf(n)
97 basis = [1, B] + [ B**i for i in range(2, m) ] + [-A]
99 ctx.dps = n
100 coeffs = ctx.pslq(basis, maxcoeff=10**10, maxsteps=1000)
102 if coeffs is None:
103 # PSLQ can't find an integer linear combination. Give up.
104 break
106 if coeffs != prev:
107 prev = coeffs
108 else:
109 # Increasing precision didn't produce anything new. Give up.
110 break
112 # We have
113 # c0 + c1*B + c2*B^2 + ... + cm-1*B^(m-1) - cm*A ~ 0.
114 # So bring cm*A to the other side, and divide through by cm,
115 # for an approximate representation of A as a polynomial in B.
116 # (We know cm != 0 since `b.minpoly` is irreducible.)
117 coeffs = [S(c)/coeffs[-1] for c in coeffs[:-1]]
119 # Throw away leading zeros.
120 while not coeffs[-1]:
121 coeffs.pop()
123 coeffs = list(reversed(coeffs))
124 h = Poly(coeffs, f.gen, domain='QQ')
126 # We only have A ~ h(B). We must check whether the relation is exact.
127 if f.compose(h).rem(g).is_zero:
128 # Now we know that h(b) is in fact equal to _some conjugate of_ a.
129 # But from the very precise approximation A ~ h(B) we can assume
130 # the conjugate is a itself.
131 return coeffs
132 else:
133 n *= 2
135 return None
138def field_isomorphism_factor(a, b):
139 """Construct field isomorphism via factorization. """
140 _, factors = factor_list(a.minpoly, extension=b)
141 for f, _ in factors:
142 if f.degree() == 1:
143 # Any linear factor f(x) represents some conjugate of a in QQ(b).
144 # We want to know whether this linear factor represents a itself.
145 # Let f = x - c
146 c = -f.rep.TC()
147 # Write c as polynomial in b
148 coeffs = c.to_sympy_list()
149 d, terms = len(coeffs) - 1, []
150 for i, coeff in enumerate(coeffs):
151 terms.append(coeff*b.root**(d - i))
152 r = Add(*terms)
153 # Check whether we got the number a
154 if a.minpoly.same_root(r, a):
155 return coeffs
157 # If none of the linear factors represented a in QQ(b), then in fact a is
158 # not an element of QQ(b).
159 return None
162@public
163def field_isomorphism(a, b, *, fast=True):
164 r"""
165 Find an embedding of one number field into another.
167 Explanation
168 ===========
170 This function looks for an isomorphism from $\mathbb{Q}(a)$ onto some
171 subfield of $\mathbb{Q}(b)$. Thus, it solves the Subfield Problem.
173 Examples
174 ========
176 >>> from sympy import sqrt, field_isomorphism, I
177 >>> print(field_isomorphism(3, sqrt(2))) # doctest: +SKIP
178 [3]
179 >>> print(field_isomorphism( I*sqrt(3), I*sqrt(3)/2)) # doctest: +SKIP
180 [2, 0]
182 Parameters
183 ==========
185 a : :py:class:`~.Expr`
186 Any expression representing an algebraic number.
187 b : :py:class:`~.Expr`
188 Any expression representing an algebraic number.
189 fast : boolean, optional (default=True)
190 If ``True``, we first attempt a potentially faster way of computing the
191 isomorphism, falling back on a slower method if this fails. If
192 ``False``, we go directly to the slower method, which is guaranteed to
193 return a result.
195 Returns
196 =======
198 List of rational numbers, or None
199 If $\mathbb{Q}(a)$ is not isomorphic to some subfield of
200 $\mathbb{Q}(b)$, then return ``None``. Otherwise, return a list of
201 rational numbers representing an element of $\mathbb{Q}(b)$ to which
202 $a$ may be mapped, in order to define a monomorphism, i.e. an
203 isomorphism from $\mathbb{Q}(a)$ to some subfield of $\mathbb{Q}(b)$.
204 The elements of the list are the coefficients of falling powers of $b$.
206 """
207 a, b = sympify(a), sympify(b)
209 if not a.is_AlgebraicNumber:
210 a = AlgebraicNumber(a)
212 if not b.is_AlgebraicNumber:
213 b = AlgebraicNumber(b)
215 a = a.to_primitive_element()
216 b = b.to_primitive_element()
218 if a == b:
219 return a.coeffs()
221 n = a.minpoly.degree()
222 m = b.minpoly.degree()
224 if n == 1:
225 return [a.root]
227 if m % n != 0:
228 return None
230 if fast:
231 try:
232 result = field_isomorphism_pslq(a, b)
234 if result is not None:
235 return result
236 except NotImplementedError:
237 pass
239 return field_isomorphism_factor(a, b)
242def _switch_domain(g, K):
243 # An algebraic relation f(a, b) = 0 over Q can also be written
244 # g(b) = 0 where g is in Q(a)[x] and h(a) = 0 where h is in Q(b)[x].
245 # This function transforms g into h where Q(b) = K.
246 frep = g.rep.inject()
247 hrep = frep.eject(K, front=True)
249 return g.new(hrep, g.gens[0])
252def _linsolve(p):
253 # Compute root of linear polynomial.
254 c, d = p.rep.rep
255 return -d/c
258@public
259def primitive_element(extension, x=None, *, ex=False, polys=False):
260 r"""
261 Find a single generator for a number field given by several generators.
263 Explanation
264 ===========
266 The basic problem is this: Given several algebraic numbers
267 $\alpha_1, \alpha_2, \ldots, \alpha_n$, find a single algebraic number
268 $\theta$ such that
269 $\mathbb{Q}(\alpha_1, \alpha_2, \ldots, \alpha_n) = \mathbb{Q}(\theta)$.
271 This function actually guarantees that $\theta$ will be a linear
272 combination of the $\alpha_i$, with non-negative integer coefficients.
274 Furthermore, if desired, this function will tell you how to express each
275 $\alpha_i$ as a $\mathbb{Q}$-linear combination of the powers of $\theta$.
277 Examples
278 ========
280 >>> from sympy import primitive_element, sqrt, S, minpoly, simplify
281 >>> from sympy.abc import x
282 >>> f, lincomb, reps = primitive_element([sqrt(2), sqrt(3)], x, ex=True)
284 Then ``lincomb`` tells us the primitive element as a linear combination of
285 the given generators ``sqrt(2)`` and ``sqrt(3)``.
287 >>> print(lincomb)
288 [1, 1]
290 This means the primtiive element is $\sqrt{2} + \sqrt{3}$.
291 Meanwhile ``f`` is the minimal polynomial for this primitive element.
293 >>> print(f)
294 x**4 - 10*x**2 + 1
295 >>> print(minpoly(sqrt(2) + sqrt(3), x))
296 x**4 - 10*x**2 + 1
298 Finally, ``reps`` (which was returned only because we set keyword arg
299 ``ex=True``) tells us how to recover each of the generators $\sqrt{2}$ and
300 $\sqrt{3}$ as $\mathbb{Q}$-linear combinations of the powers of the
301 primitive element $\sqrt{2} + \sqrt{3}$.
303 >>> print([S(r) for r in reps[0]])
304 [1/2, 0, -9/2, 0]
305 >>> theta = sqrt(2) + sqrt(3)
306 >>> print(simplify(theta**3/2 - 9*theta/2))
307 sqrt(2)
308 >>> print([S(r) for r in reps[1]])
309 [-1/2, 0, 11/2, 0]
310 >>> print(simplify(-theta**3/2 + 11*theta/2))
311 sqrt(3)
313 Parameters
314 ==========
316 extension : list of :py:class:`~.Expr`
317 Each expression must represent an algebraic number $\alpha_i$.
318 x : :py:class:`~.Symbol`, optional (default=None)
319 The desired symbol to appear in the computed minimal polynomial for the
320 primitive element $\theta$. If ``None``, we use a dummy symbol.
321 ex : boolean, optional (default=False)
322 If and only if ``True``, compute the representation of each $\alpha_i$
323 as a $\mathbb{Q}$-linear combination over the powers of $\theta$.
324 polys : boolean, optional (default=False)
325 If ``True``, return the minimal polynomial as a :py:class:`~.Poly`.
326 Otherwise return it as an :py:class:`~.Expr`.
328 Returns
329 =======
331 Pair (f, coeffs) or triple (f, coeffs, reps), where:
332 ``f`` is the minimal polynomial for the primitive element.
333 ``coeffs`` gives the primitive element as a linear combination of the
334 given generators.
335 ``reps`` is present if and only if argument ``ex=True`` was passed,
336 and is a list of lists of rational numbers. Each list gives the
337 coefficients of falling powers of the primitive element, to recover
338 one of the original, given generators.
340 """
341 if not extension:
342 raise ValueError("Cannot compute primitive element for empty extension")
343 extension = [_sympify(ext) for ext in extension]
345 if x is not None:
346 x, cls = sympify(x), Poly
347 else:
348 x, cls = Dummy('x'), PurePoly
350 if not ex:
351 gen, coeffs = extension[0], [1]
352 g = minimal_polynomial(gen, x, polys=True)
353 for ext in extension[1:]:
354 if ext.is_Rational:
355 coeffs.append(0)
356 continue
357 _, factors = factor_list(g, extension=ext)
358 g = _choose_factor(factors, x, gen)
359 s, _, g = g.sqf_norm()
360 gen += s*ext
361 coeffs.append(s)
363 if not polys:
364 return g.as_expr(), coeffs
365 else:
366 return cls(g), coeffs
368 gen, coeffs = extension[0], [1]
369 f = minimal_polynomial(gen, x, polys=True)
370 K = QQ.algebraic_field((f, gen)) # incrementally constructed field
371 reps = [K.unit] # representations of extension elements in K
372 for ext in extension[1:]:
373 if ext.is_Rational:
374 coeffs.append(0) # rational ext is not included in the expression of a primitive element
375 reps.append(K.convert(ext)) # but it is included in reps
376 continue
377 p = minimal_polynomial(ext, x, polys=True)
378 L = QQ.algebraic_field((p, ext))
379 _, factors = factor_list(f, domain=L)
380 f = _choose_factor(factors, x, gen)
381 s, g, f = f.sqf_norm()
382 gen += s*ext
383 coeffs.append(s)
384 K = QQ.algebraic_field((f, gen))
385 h = _switch_domain(g, K)
386 erep = _linsolve(h.gcd(p)) # ext as element of K
387 ogen = K.unit - s*erep # old gen as element of K
388 reps = [dup_eval(_.rep, ogen, K) for _ in reps] + [erep]
390 if K.ext.root.is_Rational: # all extensions are rational
391 H = [K.convert(_).rep for _ in extension]
392 coeffs = [0]*len(extension)
393 f = cls(x, domain=QQ)
394 else:
395 H = [_.rep for _ in reps]
396 if not polys:
397 return f.as_expr(), coeffs, H
398 else:
399 return f, coeffs, H
402@public
403def to_number_field(extension, theta=None, *, gen=None, alias=None):
404 r"""
405 Express one algebraic number in the field generated by another.
407 Explanation
408 ===========
410 Given two algebraic numbers $\eta, \theta$, this function either expresses
411 $\eta$ as an element of $\mathbb{Q}(\theta)$, or else raises an exception
412 if $\eta \not\in \mathbb{Q}(\theta)$.
414 This function is essentially just a convenience, utilizing
415 :py:func:`~.field_isomorphism` (our solution of the Subfield Problem) to
416 solve this, the Field Membership Problem.
418 As an additional convenience, this function allows you to pass a list of
419 algebraic numbers $\alpha_1, \alpha_2, \ldots, \alpha_n$ instead of $\eta$.
420 It then computes $\eta$ for you, as a solution of the Primitive Element
421 Problem, using :py:func:`~.primitive_element` on the list of $\alpha_i$.
423 Examples
424 ========
426 >>> from sympy import sqrt, to_number_field
427 >>> eta = sqrt(2)
428 >>> theta = sqrt(2) + sqrt(3)
429 >>> a = to_number_field(eta, theta)
430 >>> print(type(a))
431 <class 'sympy.core.numbers.AlgebraicNumber'>
432 >>> a.root
433 sqrt(2) + sqrt(3)
434 >>> print(a)
435 sqrt(2)
436 >>> a.coeffs()
437 [1/2, 0, -9/2, 0]
439 We get an :py:class:`~.AlgebraicNumber`, whose ``.root`` is $\theta$, whose
440 value is $\eta$, and whose ``.coeffs()`` show how to write $\eta$ as a
441 $\mathbb{Q}$-linear combination in falling powers of $\theta$.
443 Parameters
444 ==========
446 extension : :py:class:`~.Expr` or list of :py:class:`~.Expr`
447 Either the algebraic number that is to be expressed in the other field,
448 or else a list of algebraic numbers, a primitive element for which is
449 to be expressed in the other field.
450 theta : :py:class:`~.Expr`, None, optional (default=None)
451 If an :py:class:`~.Expr` representing an algebraic number, behavior is
452 as described under **Explanation**. If ``None``, then this function
453 reduces to a shorthand for calling :py:func:`~.primitive_element` on
454 ``extension`` and turning the computed primitive element into an
455 :py:class:`~.AlgebraicNumber`.
456 gen : :py:class:`~.Symbol`, None, optional (default=None)
457 If provided, this will be used as the generator symbol for the minimal
458 polynomial in the returned :py:class:`~.AlgebraicNumber`.
459 alias : str, :py:class:`~.Symbol`, None, optional (default=None)
460 If provided, this will be used as the alias symbol for the returned
461 :py:class:`~.AlgebraicNumber`.
463 Returns
464 =======
466 AlgebraicNumber
467 Belonging to $\mathbb{Q}(\theta)$ and equaling $\eta$.
469 Raises
470 ======
472 IsomorphismFailed
473 If $\eta \not\in \mathbb{Q}(\theta)$.
475 See Also
476 ========
478 field_isomorphism
479 primitive_element
481 """
482 if hasattr(extension, '__iter__'):
483 extension = list(extension)
484 else:
485 extension = [extension]
487 if len(extension) == 1 and isinstance(extension[0], tuple):
488 return AlgebraicNumber(extension[0], alias=alias)
490 minpoly, coeffs = primitive_element(extension, gen, polys=True)
491 root = sum([ coeff*ext for coeff, ext in zip(coeffs, extension) ])
493 if theta is None:
494 return AlgebraicNumber((minpoly, root), alias=alias)
495 else:
496 theta = sympify(theta)
498 if not theta.is_AlgebraicNumber:
499 theta = AlgebraicNumber(theta, gen=gen, alias=alias)
501 coeffs = field_isomorphism(root, theta)
503 if coeffs is not None:
504 return AlgebraicNumber(theta, coeffs, alias=alias)
505 else:
506 raise IsomorphismFailed(
507 "%s is not in a subfield of %s" % (root, theta.root))