Coverage for /usr/lib/python3/dist-packages/sympy/polys/numberfields/modules.py: 27%
580 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1r"""Modules in number fields.
3The classes defined here allow us to work with finitely generated, free
4modules, whose generators are algebraic numbers.
6There is an abstract base class called :py:class:`~.Module`, which has two
7concrete subclasses, :py:class:`~.PowerBasis` and :py:class:`~.Submodule`.
9Every module is defined by its basis, or set of generators:
11* For a :py:class:`~.PowerBasis`, the generators are the first $n$ powers
12 (starting with the zeroth) of an algebraic integer $\theta$ of degree $n$.
13 The :py:class:`~.PowerBasis` is constructed by passing either the minimal
14 polynomial of $\theta$, or an :py:class:`~.AlgebraicField` having $\theta$
15 as its primitive element.
17* For a :py:class:`~.Submodule`, the generators are a set of
18 $\mathbb{Q}$-linear combinations of the generators of another module. That
19 other module is then the "parent" of the :py:class:`~.Submodule`. The
20 coefficients of the $\mathbb{Q}$-linear combinations may be given by an
21 integer matrix, and a positive integer denominator. Each column of the matrix
22 defines a generator.
24>>> from sympy.polys import Poly, cyclotomic_poly, ZZ
25>>> from sympy.abc import x
26>>> from sympy.polys.matrices import DomainMatrix, DM
27>>> from sympy.polys.numberfields.modules import PowerBasis
28>>> T = Poly(cyclotomic_poly(5, x))
29>>> A = PowerBasis(T)
30>>> print(A)
31PowerBasis(x**4 + x**3 + x**2 + x + 1)
32>>> B = A.submodule_from_matrix(2 * DomainMatrix.eye(4, ZZ), denom=3)
33>>> print(B)
34Submodule[[2, 0, 0, 0], [0, 2, 0, 0], [0, 0, 2, 0], [0, 0, 0, 2]]/3
35>>> print(B.parent)
36PowerBasis(x**4 + x**3 + x**2 + x + 1)
38Thus, every module is either a :py:class:`~.PowerBasis`,
39or a :py:class:`~.Submodule`, some ancestor of which is a
40:py:class:`~.PowerBasis`. (If ``S`` is a :py:class:`~.Submodule`, then its
41ancestors are ``S.parent``, ``S.parent.parent``, and so on).
43The :py:class:`~.ModuleElement` class represents a linear combination of the
44generators of any module. Critically, the coefficients of this linear
45combination are not restricted to be integers, but may be any rational
46numbers. This is necessary so that any and all algebraic integers be
47representable, starting from the power basis in a primitive element $\theta$
48for the number field in question. For example, in a quadratic field
49$\mathbb{Q}(\sqrt{d})$ where $d \equiv 1 \mod{4}$, a denominator of $2$ is
50needed.
52A :py:class:`~.ModuleElement` can be constructed from an integer column vector
53and a denominator:
55>>> U = Poly(x**2 - 5)
56>>> M = PowerBasis(U)
57>>> e = M(DM([[1], [1]], ZZ), denom=2)
58>>> print(e)
59[1, 1]/2
60>>> print(e.module)
61PowerBasis(x**2 - 5)
63The :py:class:`~.PowerBasisElement` class is a subclass of
64:py:class:`~.ModuleElement` that represents elements of a
65:py:class:`~.PowerBasis`, and adds functionality pertinent to elements
66represented directly over powers of the primitive element $\theta$.
69Arithmetic with module elements
70===============================
72While a :py:class:`~.ModuleElement` represents a linear combination over the
73generators of a particular module, recall that every module is either a
74:py:class:`~.PowerBasis` or a descendant (along a chain of
75:py:class:`~.Submodule` objects) thereof, so that in fact every
76:py:class:`~.ModuleElement` represents an algebraic number in some field
77$\mathbb{Q}(\theta)$, where $\theta$ is the defining element of some
78:py:class:`~.PowerBasis`. It thus makes sense to talk about the number field
79to which a given :py:class:`~.ModuleElement` belongs.
81This means that any two :py:class:`~.ModuleElement` instances can be added,
82subtracted, multiplied, or divided, provided they belong to the same number
83field. Similarly, since $\mathbb{Q}$ is a subfield of every number field,
84any :py:class:`~.ModuleElement` may be added, multiplied, etc. by any
85rational number.
87>>> from sympy import QQ
88>>> from sympy.polys.numberfields.modules import to_col
89>>> T = Poly(cyclotomic_poly(5))
90>>> A = PowerBasis(T)
91>>> C = A.submodule_from_matrix(3 * DomainMatrix.eye(4, ZZ))
92>>> e = A(to_col([0, 2, 0, 0]), denom=3)
93>>> f = A(to_col([0, 0, 0, 7]), denom=5)
94>>> g = C(to_col([1, 1, 1, 1]))
95>>> e + f
96[0, 10, 0, 21]/15
97>>> e - f
98[0, 10, 0, -21]/15
99>>> e - g
100[-9, -7, -9, -9]/3
101>>> e + QQ(7, 10)
102[21, 20, 0, 0]/30
103>>> e * f
104[-14, -14, -14, -14]/15
105>>> e ** 2
106[0, 0, 4, 0]/9
107>>> f // g
108[7, 7, 7, 7]/15
109>>> f * QQ(2, 3)
110[0, 0, 0, 14]/15
112However, care must be taken with arithmetic operations on
113:py:class:`~.ModuleElement`, because the module $C$ to which the result will
114belong will be the nearest common ancestor (NCA) of the modules $A$, $B$ to
115which the two operands belong, and $C$ may be different from either or both
116of $A$ and $B$.
118>>> A = PowerBasis(T)
119>>> B = A.submodule_from_matrix(2 * DomainMatrix.eye(4, ZZ))
120>>> C = A.submodule_from_matrix(3 * DomainMatrix.eye(4, ZZ))
121>>> print((B(0) * C(0)).module == A)
122True
124Before the arithmetic operation is performed, copies of the two operands are
125automatically converted into elements of the NCA (the operands themselves are
126not modified). This upward conversion along an ancestor chain is easy: it just
127requires the successive multiplication by the defining matrix of each
128:py:class:`~.Submodule`.
130Conversely, downward conversion, i.e. representing a given
131:py:class:`~.ModuleElement` in a submodule, is also supported -- namely by
132the :py:meth:`~sympy.polys.numberfields.modules.Submodule.represent` method
133-- but is not guaranteed to succeed in general, since the given element may
134not belong to the submodule. The main circumstance in which this issue tends
135to arise is with multiplication, since modules, while closed under addition,
136need not be closed under multiplication.
139Multiplication
140--------------
142Generally speaking, a module need not be closed under multiplication, i.e. need
143not form a ring. However, many of the modules we work with in the context of
144number fields are in fact rings, and our classes do support multiplication.
146Specifically, any :py:class:`~.Module` can attempt to compute its own
147multiplication table, but this does not happen unless an attempt is made to
148multiply two :py:class:`~.ModuleElement` instances belonging to it.
150>>> A = PowerBasis(T)
151>>> print(A._mult_tab is None)
152True
153>>> a = A(0)*A(1)
154>>> print(A._mult_tab is None)
155False
157Every :py:class:`~.PowerBasis` is, by its nature, closed under multiplication,
158so instances of :py:class:`~.PowerBasis` can always successfully compute their
159multiplication table.
161When a :py:class:`~.Submodule` attempts to compute its multiplication table,
162it converts each of its own generators into elements of its parent module,
163multiplies them there, in every possible pairing, and then tries to
164represent the results in itself, i.e. as $\mathbb{Z}$-linear combinations
165over its own generators. This will succeed if and only if the submodule is
166in fact closed under multiplication.
169Module Homomorphisms
170====================
172Many important number theoretic algorithms require the calculation of the
173kernel of one or more module homomorphisms. Accordingly we have several
174lightweight classes, :py:class:`~.ModuleHomomorphism`,
175:py:class:`~.ModuleEndomorphism`, :py:class:`~.InnerEndomorphism`, and
176:py:class:`~.EndomorphismRing`, which provide the minimal necessary machinery
177to support this.
179"""
181from sympy.core.numbers import igcd, ilcm
182from sympy.core.symbol import Dummy
183from sympy.polys.polyclasses import ANP
184from sympy.polys.polytools import Poly
185from sympy.polys.densetools import dup_clear_denoms
186from sympy.polys.domains.algebraicfield import AlgebraicField
187from sympy.polys.domains.finitefield import FF
188from sympy.polys.domains.rationalfield import QQ
189from sympy.polys.domains.integerring import ZZ
190from sympy.polys.matrices.domainmatrix import DomainMatrix
191from sympy.polys.matrices.exceptions import DMBadInputError
192from sympy.polys.matrices.normalforms import hermite_normal_form
193from sympy.polys.polyerrors import CoercionFailed, UnificationFailed
194from sympy.polys.polyutils import IntegerPowerable
195from .exceptions import ClosureFailure, MissingUnityError, StructureError
196from .utilities import AlgIntPowers, is_rat, get_num_denom
199def to_col(coeffs):
200 r"""Transform a list of integer coefficients into a column vector."""
201 return DomainMatrix([[ZZ(c) for c in coeffs]], (1, len(coeffs)), ZZ).transpose()
204class Module:
205 """
206 Generic finitely-generated module.
208 This is an abstract base class, and should not be instantiated directly.
209 The two concrete subclasses are :py:class:`~.PowerBasis` and
210 :py:class:`~.Submodule`.
212 Every :py:class:`~.Submodule` is derived from another module, referenced
213 by its ``parent`` attribute. If ``S`` is a submodule, then we refer to
214 ``S.parent``, ``S.parent.parent``, and so on, as the "ancestors" of
215 ``S``. Thus, every :py:class:`~.Module` is either a
216 :py:class:`~.PowerBasis` or a :py:class:`~.Submodule`, some ancestor of
217 which is a :py:class:`~.PowerBasis`.
218 """
220 @property
221 def n(self):
222 """The number of generators of this module."""
223 raise NotImplementedError
225 def mult_tab(self):
226 """
227 Get the multiplication table for this module (if closed under mult).
229 Explanation
230 ===========
232 Computes a dictionary ``M`` of dictionaries of lists, representing the
233 upper triangular half of the multiplication table.
235 In other words, if ``0 <= i <= j < self.n``, then ``M[i][j]`` is the
236 list ``c`` of coefficients such that
237 ``g[i] * g[j] == sum(c[k]*g[k], k in range(self.n))``,
238 where ``g`` is the list of generators of this module.
240 If ``j < i`` then ``M[i][j]`` is undefined.
242 Examples
243 ========
245 >>> from sympy.polys import Poly, cyclotomic_poly
246 >>> from sympy.polys.numberfields.modules import PowerBasis
247 >>> T = Poly(cyclotomic_poly(5))
248 >>> A = PowerBasis(T)
249 >>> print(A.mult_tab()) # doctest: +SKIP
250 {0: {0: [1, 0, 0, 0], 1: [0, 1, 0, 0], 2: [0, 0, 1, 0], 3: [0, 0, 0, 1]},
251 1: {1: [0, 0, 1, 0], 2: [0, 0, 0, 1], 3: [-1, -1, -1, -1]},
252 2: {2: [-1, -1, -1, -1], 3: [1, 0, 0, 0]},
253 3: {3: [0, 1, 0, 0]}}
255 Returns
256 =======
258 dict of dict of lists
260 Raises
261 ======
263 ClosureFailure
264 If the module is not closed under multiplication.
266 """
267 raise NotImplementedError
269 @property
270 def parent(self):
271 """
272 The parent module, if any, for this module.
274 Explanation
275 ===========
277 For a :py:class:`~.Submodule` this is its ``parent`` attribute; for a
278 :py:class:`~.PowerBasis` this is ``None``.
280 Returns
281 =======
283 :py:class:`~.Module`, ``None``
285 See Also
286 ========
288 Module
290 """
291 return None
293 def represent(self, elt):
294 r"""
295 Represent a module element as an integer-linear combination over the
296 generators of this module.
298 Explanation
299 ===========
301 In our system, to "represent" always means to write a
302 :py:class:`~.ModuleElement` as a :ref:`ZZ`-linear combination over the
303 generators of the present :py:class:`~.Module`. Furthermore, the
304 incoming :py:class:`~.ModuleElement` must belong to an ancestor of
305 the present :py:class:`~.Module` (or to the present
306 :py:class:`~.Module` itself).
308 The most common application is to represent a
309 :py:class:`~.ModuleElement` in a :py:class:`~.Submodule`. For example,
310 this is involved in computing multiplication tables.
312 On the other hand, representing in a :py:class:`~.PowerBasis` is an
313 odd case, and one which tends not to arise in practice, except for
314 example when using a :py:class:`~.ModuleEndomorphism` on a
315 :py:class:`~.PowerBasis`.
317 In such a case, (1) the incoming :py:class:`~.ModuleElement` must
318 belong to the :py:class:`~.PowerBasis` itself (since the latter has no
319 proper ancestors) and (2) it is "representable" iff it belongs to
320 $\mathbb{Z}[\theta]$ (although generally a
321 :py:class:`~.PowerBasisElement` may represent any element of
322 $\mathbb{Q}(\theta)$, i.e. any algebraic number).
324 Examples
325 ========
327 >>> from sympy import Poly, cyclotomic_poly
328 >>> from sympy.polys.numberfields.modules import PowerBasis, to_col
329 >>> from sympy.abc import zeta
330 >>> T = Poly(cyclotomic_poly(5))
331 >>> A = PowerBasis(T)
332 >>> a = A(to_col([2, 4, 6, 8]))
334 The :py:class:`~.ModuleElement` ``a`` has all even coefficients.
335 If we represent ``a`` in the submodule ``B = 2*A``, the coefficients in
336 the column vector will be halved:
338 >>> B = A.submodule_from_gens([2*A(i) for i in range(4)])
339 >>> b = B.represent(a)
340 >>> print(b.transpose()) # doctest: +SKIP
341 DomainMatrix([[1, 2, 3, 4]], (1, 4), ZZ)
343 However, the element of ``B`` so defined still represents the same
344 algebraic number:
346 >>> print(a.poly(zeta).as_expr())
347 8*zeta**3 + 6*zeta**2 + 4*zeta + 2
348 >>> print(B(b).over_power_basis().poly(zeta).as_expr())
349 8*zeta**3 + 6*zeta**2 + 4*zeta + 2
351 Parameters
352 ==========
354 elt : :py:class:`~.ModuleElement`
355 The module element to be represented. Must belong to some ancestor
356 module of this module (including this module itself).
358 Returns
359 =======
361 :py:class:`~.DomainMatrix` over :ref:`ZZ`
362 This will be a column vector, representing the coefficients of a
363 linear combination of this module's generators, which equals the
364 given element.
366 Raises
367 ======
369 ClosureFailure
370 If the given element cannot be represented as a :ref:`ZZ`-linear
371 combination over this module.
373 See Also
374 ========
376 .Submodule.represent
377 .PowerBasis.represent
379 """
380 raise NotImplementedError
382 def ancestors(self, include_self=False):
383 """
384 Return the list of ancestor modules of this module, from the
385 foundational :py:class:`~.PowerBasis` downward, optionally including
386 ``self``.
388 See Also
389 ========
391 Module
393 """
394 c = self.parent
395 a = [] if c is None else c.ancestors(include_self=True)
396 if include_self:
397 a.append(self)
398 return a
400 def power_basis_ancestor(self):
401 """
402 Return the :py:class:`~.PowerBasis` that is an ancestor of this module.
404 See Also
405 ========
407 Module
409 """
410 if isinstance(self, PowerBasis):
411 return self
412 c = self.parent
413 if c is not None:
414 return c.power_basis_ancestor()
415 return None
417 def nearest_common_ancestor(self, other):
418 """
419 Locate the nearest common ancestor of this module and another.
421 Returns
422 =======
424 :py:class:`~.Module`, ``None``
426 See Also
427 ========
429 Module
431 """
432 sA = self.ancestors(include_self=True)
433 oA = other.ancestors(include_self=True)
434 nca = None
435 for sa, oa in zip(sA, oA):
436 if sa == oa:
437 nca = sa
438 else:
439 break
440 return nca
442 @property
443 def number_field(self):
444 r"""
445 Return the associated :py:class:`~.AlgebraicField`, if any.
447 Explanation
448 ===========
450 A :py:class:`~.PowerBasis` can be constructed on a :py:class:`~.Poly`
451 $f$ or on an :py:class:`~.AlgebraicField` $K$. In the latter case, the
452 :py:class:`~.PowerBasis` and all its descendant modules will return $K$
453 as their ``.number_field`` property, while in the former case they will
454 all return ``None``.
456 Returns
457 =======
459 :py:class:`~.AlgebraicField`, ``None``
461 """
462 return self.power_basis_ancestor().number_field
464 def is_compat_col(self, col):
465 """Say whether *col* is a suitable column vector for this module."""
466 return isinstance(col, DomainMatrix) and col.shape == (self.n, 1) and col.domain.is_ZZ
468 def __call__(self, spec, denom=1):
469 r"""
470 Generate a :py:class:`~.ModuleElement` belonging to this module.
472 Examples
473 ========
475 >>> from sympy.polys import Poly, cyclotomic_poly
476 >>> from sympy.polys.numberfields.modules import PowerBasis, to_col
477 >>> T = Poly(cyclotomic_poly(5))
478 >>> A = PowerBasis(T)
479 >>> e = A(to_col([1, 2, 3, 4]), denom=3)
480 >>> print(e) # doctest: +SKIP
481 [1, 2, 3, 4]/3
482 >>> f = A(2)
483 >>> print(f) # doctest: +SKIP
484 [0, 0, 1, 0]
486 Parameters
487 ==========
489 spec : :py:class:`~.DomainMatrix`, int
490 Specifies the numerators of the coefficients of the
491 :py:class:`~.ModuleElement`. Can be either a column vector over
492 :ref:`ZZ`, whose length must equal the number $n$ of generators of
493 this module, or else an integer ``j``, $0 \leq j < n$, which is a
494 shorthand for column $j$ of $I_n$, the $n \times n$ identity
495 matrix.
496 denom : int, optional (default=1)
497 Denominator for the coefficients of the
498 :py:class:`~.ModuleElement`.
500 Returns
501 =======
503 :py:class:`~.ModuleElement`
504 The coefficients are the entries of the *spec* vector, divided by
505 *denom*.
507 """
508 if isinstance(spec, int) and 0 <= spec < self.n:
509 spec = DomainMatrix.eye(self.n, ZZ)[:, spec].to_dense()
510 if not self.is_compat_col(spec):
511 raise ValueError('Compatible column vector required.')
512 return make_mod_elt(self, spec, denom=denom)
514 def starts_with_unity(self):
515 """Say whether the module's first generator equals unity."""
516 raise NotImplementedError
518 def basis_elements(self):
519 """
520 Get list of :py:class:`~.ModuleElement` being the generators of this
521 module.
522 """
523 return [self(j) for j in range(self.n)]
525 def zero(self):
526 """Return a :py:class:`~.ModuleElement` representing zero."""
527 return self(0) * 0
529 def one(self):
530 """
531 Return a :py:class:`~.ModuleElement` representing unity,
532 and belonging to the first ancestor of this module (including
533 itself) that starts with unity.
534 """
535 return self.element_from_rational(1)
537 def element_from_rational(self, a):
538 """
539 Return a :py:class:`~.ModuleElement` representing a rational number.
541 Explanation
542 ===========
544 The returned :py:class:`~.ModuleElement` will belong to the first
545 module on this module's ancestor chain (including this module
546 itself) that starts with unity.
548 Examples
549 ========
551 >>> from sympy.polys import Poly, cyclotomic_poly, QQ
552 >>> from sympy.polys.numberfields.modules import PowerBasis
553 >>> T = Poly(cyclotomic_poly(5))
554 >>> A = PowerBasis(T)
555 >>> a = A.element_from_rational(QQ(2, 3))
556 >>> print(a) # doctest: +SKIP
557 [2, 0, 0, 0]/3
559 Parameters
560 ==========
562 a : int, :ref:`ZZ`, :ref:`QQ`
564 Returns
565 =======
567 :py:class:`~.ModuleElement`
569 """
570 raise NotImplementedError
572 def submodule_from_gens(self, gens, hnf=True, hnf_modulus=None):
573 """
574 Form the submodule generated by a list of :py:class:`~.ModuleElement`
575 belonging to this module.
577 Examples
578 ========
580 >>> from sympy.polys import Poly, cyclotomic_poly
581 >>> from sympy.polys.numberfields.modules import PowerBasis
582 >>> T = Poly(cyclotomic_poly(5))
583 >>> A = PowerBasis(T)
584 >>> gens = [A(0), 2*A(1), 3*A(2), 4*A(3)//5]
585 >>> B = A.submodule_from_gens(gens)
586 >>> print(B) # doctest: +SKIP
587 Submodule[[5, 0, 0, 0], [0, 10, 0, 0], [0, 0, 15, 0], [0, 0, 0, 4]]/5
589 Parameters
590 ==========
592 gens : list of :py:class:`~.ModuleElement` belonging to this module.
593 hnf : boolean, optional (default=True)
594 If True, we will reduce the matrix into Hermite Normal Form before
595 forming the :py:class:`~.Submodule`.
596 hnf_modulus : int, None, optional (default=None)
597 Modulus for use in the HNF reduction algorithm. See
598 :py:func:`~sympy.polys.matrices.normalforms.hermite_normal_form`.
600 Returns
601 =======
603 :py:class:`~.Submodule`
605 See Also
606 ========
608 submodule_from_matrix
610 """
611 if not all(g.module == self for g in gens):
612 raise ValueError('Generators must belong to this module.')
613 n = len(gens)
614 if n == 0:
615 raise ValueError('Need at least one generator.')
616 m = gens[0].n
617 d = gens[0].denom if n == 1 else ilcm(*[g.denom for g in gens])
618 B = DomainMatrix.zeros((m, 0), ZZ).hstack(*[(d // g.denom) * g.col for g in gens])
619 if hnf:
620 B = hermite_normal_form(B, D=hnf_modulus)
621 return self.submodule_from_matrix(B, denom=d)
623 def submodule_from_matrix(self, B, denom=1):
624 """
625 Form the submodule generated by the elements of this module indicated
626 by the columns of a matrix, with an optional denominator.
628 Examples
629 ========
631 >>> from sympy.polys import Poly, cyclotomic_poly, ZZ
632 >>> from sympy.polys.matrices import DM
633 >>> from sympy.polys.numberfields.modules import PowerBasis
634 >>> T = Poly(cyclotomic_poly(5))
635 >>> A = PowerBasis(T)
636 >>> B = A.submodule_from_matrix(DM([
637 ... [0, 10, 0, 0],
638 ... [0, 0, 7, 0],
639 ... ], ZZ).transpose(), denom=15)
640 >>> print(B) # doctest: +SKIP
641 Submodule[[0, 10, 0, 0], [0, 0, 7, 0]]/15
643 Parameters
644 ==========
646 B : :py:class:`~.DomainMatrix` over :ref:`ZZ`
647 Each column gives the numerators of the coefficients of one
648 generator of the submodule. Thus, the number of rows of *B* must
649 equal the number of generators of the present module.
650 denom : int, optional (default=1)
651 Common denominator for all generators of the submodule.
653 Returns
654 =======
656 :py:class:`~.Submodule`
658 Raises
659 ======
661 ValueError
662 If the given matrix *B* is not over :ref:`ZZ` or its number of rows
663 does not equal the number of generators of the present module.
665 See Also
666 ========
668 submodule_from_gens
670 """
671 m, n = B.shape
672 if not B.domain.is_ZZ:
673 raise ValueError('Matrix must be over ZZ.')
674 if not m == self.n:
675 raise ValueError('Matrix row count must match base module.')
676 return Submodule(self, B, denom=denom)
678 def whole_submodule(self):
679 """
680 Return a submodule equal to this entire module.
682 Explanation
683 ===========
685 This is useful when you have a :py:class:`~.PowerBasis` and want to
686 turn it into a :py:class:`~.Submodule` (in order to use methods
687 belonging to the latter).
689 """
690 B = DomainMatrix.eye(self.n, ZZ)
691 return self.submodule_from_matrix(B)
693 def endomorphism_ring(self):
694 """Form the :py:class:`~.EndomorphismRing` for this module."""
695 return EndomorphismRing(self)
698class PowerBasis(Module):
699 """The module generated by the powers of an algebraic integer."""
701 def __init__(self, T):
702 """
703 Parameters
704 ==========
706 T : :py:class:`~.Poly`, :py:class:`~.AlgebraicField`
707 Either (1) the monic, irreducible, univariate polynomial over
708 :ref:`ZZ`, a root of which is the generator of the power basis,
709 or (2) an :py:class:`~.AlgebraicField` whose primitive element
710 is the generator of the power basis.
712 """
713 K = None
714 if isinstance(T, AlgebraicField):
715 K, T = T, T.ext.minpoly_of_element()
716 # Sometimes incoming Polys are formally over QQ, although all their
717 # coeffs are integral. We want them to be formally over ZZ.
718 T = T.set_domain(ZZ)
719 self.K = K
720 self.T = T
721 self._n = T.degree()
722 self._mult_tab = None
724 @property
725 def number_field(self):
726 return self.K
728 def __repr__(self):
729 return f'PowerBasis({self.T.as_expr()})'
731 def __eq__(self, other):
732 if isinstance(other, PowerBasis):
733 return self.T == other.T
734 return NotImplemented
736 @property
737 def n(self):
738 return self._n
740 def mult_tab(self):
741 if self._mult_tab is None:
742 self.compute_mult_tab()
743 return self._mult_tab
745 def compute_mult_tab(self):
746 theta_pow = AlgIntPowers(self.T)
747 M = {}
748 n = self.n
749 for u in range(n):
750 M[u] = {}
751 for v in range(u, n):
752 M[u][v] = theta_pow[u + v]
753 self._mult_tab = M
755 def represent(self, elt):
756 r"""
757 Represent a module element as an integer-linear combination over the
758 generators of this module.
760 See Also
761 ========
763 .Module.represent
764 .Submodule.represent
766 """
767 if elt.module == self and elt.denom == 1:
768 return elt.column()
769 else:
770 raise ClosureFailure('Element not representable in ZZ[theta].')
772 def starts_with_unity(self):
773 return True
775 def element_from_rational(self, a):
776 return self(0) * a
778 def element_from_poly(self, f):
779 """
780 Produce an element of this module, representing *f* after reduction mod
781 our defining minimal polynomial.
783 Parameters
784 ==========
786 f : :py:class:`~.Poly` over :ref:`ZZ` in same var as our defining poly.
788 Returns
789 =======
791 :py:class:`~.PowerBasisElement`
793 """
794 n, k = self.n, f.degree()
795 if k >= n:
796 f = f % self.T
797 if f == 0:
798 return self.zero()
799 d, c = dup_clear_denoms(f.rep.rep, QQ, convert=True)
800 c = list(reversed(c))
801 ell = len(c)
802 z = [ZZ(0)] * (n - ell)
803 col = to_col(c + z)
804 return self(col, denom=d)
806 def _element_from_rep_and_mod(self, rep, mod):
807 """
808 Produce a PowerBasisElement representing a given algebraic number.
810 Parameters
811 ==========
813 rep : list of coeffs
814 Represents the number as polynomial in the primitive element of the
815 field.
817 mod : list of coeffs
818 Represents the minimal polynomial of the primitive element of the
819 field.
821 Returns
822 =======
824 :py:class:`~.PowerBasisElement`
826 """
827 if mod != self.T.rep.rep:
828 raise UnificationFailed('Element does not appear to be in the same field.')
829 return self.element_from_poly(Poly(rep, self.T.gen))
831 def element_from_ANP(self, a):
832 """Convert an ANP into a PowerBasisElement. """
833 return self._element_from_rep_and_mod(a.rep, a.mod)
835 def element_from_alg_num(self, a):
836 """Convert an AlgebraicNumber into a PowerBasisElement. """
837 return self._element_from_rep_and_mod(a.rep.rep, a.minpoly.rep.rep)
840class Submodule(Module, IntegerPowerable):
841 """A submodule of another module."""
843 def __init__(self, parent, matrix, denom=1, mult_tab=None):
844 """
845 Parameters
846 ==========
848 parent : :py:class:`~.Module`
849 The module from which this one is derived.
850 matrix : :py:class:`~.DomainMatrix` over :ref:`ZZ`
851 The matrix whose columns define this submodule's generators as
852 linear combinations over the parent's generators.
853 denom : int, optional (default=1)
854 Denominator for the coefficients given by the matrix.
855 mult_tab : dict, ``None``, optional
856 If already known, the multiplication table for this module may be
857 supplied.
859 """
860 self._parent = parent
861 self._matrix = matrix
862 self._denom = denom
863 self._mult_tab = mult_tab
864 self._n = matrix.shape[1]
865 self._QQ_matrix = None
866 self._starts_with_unity = None
867 self._is_sq_maxrank_HNF = None
869 def __repr__(self):
870 r = 'Submodule' + repr(self.matrix.transpose().to_Matrix().tolist())
871 if self.denom > 1:
872 r += f'/{self.denom}'
873 return r
875 def reduced(self):
876 """
877 Produce a reduced version of this submodule.
879 Explanation
880 ===========
882 In the reduced version, it is guaranteed that 1 is the only positive
883 integer dividing both the submodule's denominator, and every entry in
884 the submodule's matrix.
886 Returns
887 =======
889 :py:class:`~.Submodule`
891 """
892 if self.denom == 1:
893 return self
894 g = igcd(self.denom, *self.coeffs)
895 if g == 1:
896 return self
897 return type(self)(self.parent, (self.matrix / g).convert_to(ZZ), denom=self.denom // g, mult_tab=self._mult_tab)
899 def discard_before(self, r):
900 """
901 Produce a new module by discarding all generators before a given
902 index *r*.
903 """
904 W = self.matrix[:, r:]
905 s = self.n - r
906 M = None
907 mt = self._mult_tab
908 if mt is not None:
909 M = {}
910 for u in range(s):
911 M[u] = {}
912 for v in range(u, s):
913 M[u][v] = mt[r + u][r + v][r:]
914 return Submodule(self.parent, W, denom=self.denom, mult_tab=M)
916 @property
917 def n(self):
918 return self._n
920 def mult_tab(self):
921 if self._mult_tab is None:
922 self.compute_mult_tab()
923 return self._mult_tab
925 def compute_mult_tab(self):
926 gens = self.basis_element_pullbacks()
927 M = {}
928 n = self.n
929 for u in range(n):
930 M[u] = {}
931 for v in range(u, n):
932 M[u][v] = self.represent(gens[u] * gens[v]).flat()
933 self._mult_tab = M
935 @property
936 def parent(self):
937 return self._parent
939 @property
940 def matrix(self):
941 return self._matrix
943 @property
944 def coeffs(self):
945 return self.matrix.flat()
947 @property
948 def denom(self):
949 return self._denom
951 @property
952 def QQ_matrix(self):
953 """
954 :py:class:`~.DomainMatrix` over :ref:`QQ`, equal to
955 ``self.matrix / self.denom``, and guaranteed to be dense.
957 Explanation
958 ===========
960 Depending on how it is formed, a :py:class:`~.DomainMatrix` may have
961 an internal representation that is sparse or dense. We guarantee a
962 dense representation here, so that tests for equivalence of submodules
963 always come out as expected.
965 Examples
966 ========
968 >>> from sympy.polys import Poly, cyclotomic_poly, ZZ
969 >>> from sympy.abc import x
970 >>> from sympy.polys.matrices import DomainMatrix
971 >>> from sympy.polys.numberfields.modules import PowerBasis
972 >>> T = Poly(cyclotomic_poly(5, x))
973 >>> A = PowerBasis(T)
974 >>> B = A.submodule_from_matrix(3*DomainMatrix.eye(4, ZZ), denom=6)
975 >>> C = A.submodule_from_matrix(DomainMatrix.eye(4, ZZ), denom=2)
976 >>> print(B.QQ_matrix == C.QQ_matrix)
977 True
979 Returns
980 =======
982 :py:class:`~.DomainMatrix` over :ref:`QQ`
984 """
985 if self._QQ_matrix is None:
986 self._QQ_matrix = (self.matrix / self.denom).to_dense()
987 return self._QQ_matrix
989 def starts_with_unity(self):
990 if self._starts_with_unity is None:
991 self._starts_with_unity = self(0).equiv(1)
992 return self._starts_with_unity
994 def is_sq_maxrank_HNF(self):
995 if self._is_sq_maxrank_HNF is None:
996 self._is_sq_maxrank_HNF = is_sq_maxrank_HNF(self._matrix)
997 return self._is_sq_maxrank_HNF
999 def is_power_basis_submodule(self):
1000 return isinstance(self.parent, PowerBasis)
1002 def element_from_rational(self, a):
1003 if self.starts_with_unity():
1004 return self(0) * a
1005 else:
1006 return self.parent.element_from_rational(a)
1008 def basis_element_pullbacks(self):
1009 """
1010 Return list of this submodule's basis elements as elements of the
1011 submodule's parent module.
1012 """
1013 return [e.to_parent() for e in self.basis_elements()]
1015 def represent(self, elt):
1016 """
1017 Represent a module element as an integer-linear combination over the
1018 generators of this module.
1020 See Also
1021 ========
1023 .Module.represent
1024 .PowerBasis.represent
1026 """
1027 if elt.module == self:
1028 return elt.column()
1029 elif elt.module == self.parent:
1030 try:
1031 # The given element should be a ZZ-linear combination over our
1032 # basis vectors; however, due to the presence of denominators,
1033 # we need to solve over QQ.
1034 A = self.QQ_matrix
1035 b = elt.QQ_col
1036 x = A._solve(b)[0].transpose()
1037 x = x.convert_to(ZZ)
1038 except DMBadInputError:
1039 raise ClosureFailure('Element outside QQ-span of this basis.')
1040 except CoercionFailed:
1041 raise ClosureFailure('Element in QQ-span but not ZZ-span of this basis.')
1042 return x
1043 elif isinstance(self.parent, Submodule):
1044 coeffs_in_parent = self.parent.represent(elt)
1045 parent_element = self.parent(coeffs_in_parent)
1046 return self.represent(parent_element)
1047 else:
1048 raise ClosureFailure('Element outside ancestor chain of this module.')
1050 def is_compat_submodule(self, other):
1051 return isinstance(other, Submodule) and other.parent == self.parent
1053 def __eq__(self, other):
1054 if self.is_compat_submodule(other):
1055 return other.QQ_matrix == self.QQ_matrix
1056 return NotImplemented
1058 def add(self, other, hnf=True, hnf_modulus=None):
1059 """
1060 Add this :py:class:`~.Submodule` to another.
1062 Explanation
1063 ===========
1065 This represents the module generated by the union of the two modules'
1066 sets of generators.
1068 Parameters
1069 ==========
1071 other : :py:class:`~.Submodule`
1072 hnf : boolean, optional (default=True)
1073 If ``True``, reduce the matrix of the combined module to its
1074 Hermite Normal Form.
1075 hnf_modulus : :ref:`ZZ`, None, optional
1076 If a positive integer is provided, use this as modulus in the
1077 HNF reduction. See
1078 :py:func:`~sympy.polys.matrices.normalforms.hermite_normal_form`.
1080 Returns
1081 =======
1083 :py:class:`~.Submodule`
1085 """
1086 d, e = self.denom, other.denom
1087 m = ilcm(d, e)
1088 a, b = m // d, m // e
1089 B = (a * self.matrix).hstack(b * other.matrix)
1090 if hnf:
1091 B = hermite_normal_form(B, D=hnf_modulus)
1092 return self.parent.submodule_from_matrix(B, denom=m)
1094 def __add__(self, other):
1095 if self.is_compat_submodule(other):
1096 return self.add(other)
1097 return NotImplemented
1099 __radd__ = __add__
1101 def mul(self, other, hnf=True, hnf_modulus=None):
1102 """
1103 Multiply this :py:class:`~.Submodule` by a rational number, a
1104 :py:class:`~.ModuleElement`, or another :py:class:`~.Submodule`.
1106 Explanation
1107 ===========
1109 To multiply by a rational number or :py:class:`~.ModuleElement` means
1110 to form the submodule whose generators are the products of this
1111 quantity with all the generators of the present submodule.
1113 To multiply by another :py:class:`~.Submodule` means to form the
1114 submodule whose generators are all the products of one generator from
1115 the one submodule, and one generator from the other.
1117 Parameters
1118 ==========
1120 other : int, :ref:`ZZ`, :ref:`QQ`, :py:class:`~.ModuleElement`, :py:class:`~.Submodule`
1121 hnf : boolean, optional (default=True)
1122 If ``True``, reduce the matrix of the product module to its
1123 Hermite Normal Form.
1124 hnf_modulus : :ref:`ZZ`, None, optional
1125 If a positive integer is provided, use this as modulus in the
1126 HNF reduction. See
1127 :py:func:`~sympy.polys.matrices.normalforms.hermite_normal_form`.
1129 Returns
1130 =======
1132 :py:class:`~.Submodule`
1134 """
1135 if is_rat(other):
1136 a, b = get_num_denom(other)
1137 if a == b == 1:
1138 return self
1139 else:
1140 return Submodule(self.parent,
1141 self.matrix * a, denom=self.denom * b,
1142 mult_tab=None).reduced()
1143 elif isinstance(other, ModuleElement) and other.module == self.parent:
1144 # The submodule is multiplied by an element of the parent module.
1145 # We presume this means we want a new submodule of the parent module.
1146 gens = [other * e for e in self.basis_element_pullbacks()]
1147 return self.parent.submodule_from_gens(gens, hnf=hnf, hnf_modulus=hnf_modulus)
1148 elif self.is_compat_submodule(other):
1149 # This case usually means you're multiplying ideals, and want another
1150 # ideal, i.e. another submodule of the same parent module.
1151 alphas, betas = self.basis_element_pullbacks(), other.basis_element_pullbacks()
1152 gens = [a * b for a in alphas for b in betas]
1153 return self.parent.submodule_from_gens(gens, hnf=hnf, hnf_modulus=hnf_modulus)
1154 return NotImplemented
1156 def __mul__(self, other):
1157 return self.mul(other)
1159 __rmul__ = __mul__
1161 def _first_power(self):
1162 return self
1164 def reduce_element(self, elt):
1165 r"""
1166 If this submodule $B$ has defining matrix $W$ in square, maximal-rank
1167 Hermite normal form, then, given an element $x$ of the parent module
1168 $A$, we produce an element $y \in A$ such that $x - y \in B$, and the
1169 $i$th coordinate of $y$ satisfies $0 \leq y_i < w_{i,i}$. This
1170 representative $y$ is unique, in the sense that every element of
1171 the coset $x + B$ reduces to it under this procedure.
1173 Explanation
1174 ===========
1176 In the special case where $A$ is a power basis for a number field $K$,
1177 and $B$ is a submodule representing an ideal $I$, this operation
1178 represents one of a few important ways of reducing an element of $K$
1179 modulo $I$ to obtain a "small" representative. See [Cohen00]_ Section
1180 1.4.3.
1182 Examples
1183 ========
1185 >>> from sympy import QQ, Poly, symbols
1186 >>> t = symbols('t')
1187 >>> k = QQ.alg_field_from_poly(Poly(t**3 + t**2 - 2*t + 8))
1188 >>> Zk = k.maximal_order()
1189 >>> A = Zk.parent
1190 >>> B = (A(2) - 3*A(0))*Zk
1191 >>> B.reduce_element(A(2))
1192 [3, 0, 0]
1194 Parameters
1195 ==========
1197 elt : :py:class:`~.ModuleElement`
1198 An element of this submodule's parent module.
1200 Returns
1201 =======
1203 elt : :py:class:`~.ModuleElement`
1204 An element of this submodule's parent module.
1206 Raises
1207 ======
1209 NotImplementedError
1210 If the given :py:class:`~.ModuleElement` does not belong to this
1211 submodule's parent module.
1212 StructureError
1213 If this submodule's defining matrix is not in square, maximal-rank
1214 Hermite normal form.
1216 References
1217 ==========
1219 .. [Cohen00] Cohen, H. *Advanced Topics in Computational Number
1220 Theory.*
1222 """
1223 if not elt.module == self.parent:
1224 raise NotImplementedError
1225 if not self.is_sq_maxrank_HNF():
1226 msg = "Reduction not implemented unless matrix square max-rank HNF"
1227 raise StructureError(msg)
1228 B = self.basis_element_pullbacks()
1229 a = elt
1230 for i in range(self.n - 1, -1, -1):
1231 b = B[i]
1232 q = a.coeffs[i]*b.denom // (b.coeffs[i]*a.denom)
1233 a -= q*b
1234 return a
1237def is_sq_maxrank_HNF(dm):
1238 r"""
1239 Say whether a :py:class:`~.DomainMatrix` is in that special case of Hermite
1240 Normal Form, in which the matrix is also square and of maximal rank.
1242 Explanation
1243 ===========
1245 We commonly work with :py:class:`~.Submodule` instances whose matrix is in
1246 this form, and it can be useful to be able to check that this condition is
1247 satisfied.
1249 For example this is the case with the :py:class:`~.Submodule` ``ZK``
1250 returned by :py:func:`~sympy.polys.numberfields.basis.round_two`, which
1251 represents the maximal order in a number field, and with ideals formed
1252 therefrom, such as ``2 * ZK``.
1254 """
1255 if dm.domain.is_ZZ and dm.is_square and dm.is_upper:
1256 n = dm.shape[0]
1257 for i in range(n):
1258 d = dm[i, i].element
1259 if d <= 0:
1260 return False
1261 for j in range(i + 1, n):
1262 if not (0 <= dm[i, j].element < d):
1263 return False
1264 return True
1265 return False
1268def make_mod_elt(module, col, denom=1):
1269 r"""
1270 Factory function which builds a :py:class:`~.ModuleElement`, but ensures
1271 that it is a :py:class:`~.PowerBasisElement` if the module is a
1272 :py:class:`~.PowerBasis`.
1273 """
1274 if isinstance(module, PowerBasis):
1275 return PowerBasisElement(module, col, denom=denom)
1276 else:
1277 return ModuleElement(module, col, denom=denom)
1280class ModuleElement(IntegerPowerable):
1281 r"""
1282 Represents an element of a :py:class:`~.Module`.
1284 NOTE: Should not be constructed directly. Use the
1285 :py:meth:`~.Module.__call__` method or the :py:func:`make_mod_elt()`
1286 factory function instead.
1287 """
1289 def __init__(self, module, col, denom=1):
1290 """
1291 Parameters
1292 ==========
1294 module : :py:class:`~.Module`
1295 The module to which this element belongs.
1296 col : :py:class:`~.DomainMatrix` over :ref:`ZZ`
1297 Column vector giving the numerators of the coefficients of this
1298 element.
1299 denom : int, optional (default=1)
1300 Denominator for the coefficients of this element.
1302 """
1303 self.module = module
1304 self.col = col
1305 self.denom = denom
1306 self._QQ_col = None
1308 def __repr__(self):
1309 r = str([int(c) for c in self.col.flat()])
1310 if self.denom > 1:
1311 r += f'/{self.denom}'
1312 return r
1314 def reduced(self):
1315 """
1316 Produce a reduced version of this ModuleElement, i.e. one in which the
1317 gcd of the denominator together with all numerator coefficients is 1.
1318 """
1319 if self.denom == 1:
1320 return self
1321 g = igcd(self.denom, *self.coeffs)
1322 if g == 1:
1323 return self
1324 return type(self)(self.module,
1325 (self.col / g).convert_to(ZZ),
1326 denom=self.denom // g)
1328 def reduced_mod_p(self, p):
1329 """
1330 Produce a version of this :py:class:`~.ModuleElement` in which all
1331 numerator coefficients have been reduced mod *p*.
1332 """
1333 return make_mod_elt(self.module,
1334 self.col.convert_to(FF(p)).convert_to(ZZ),
1335 denom=self.denom)
1337 @classmethod
1338 def from_int_list(cls, module, coeffs, denom=1):
1339 """
1340 Make a :py:class:`~.ModuleElement` from a list of ints (instead of a
1341 column vector).
1342 """
1343 col = to_col(coeffs)
1344 return cls(module, col, denom=denom)
1346 @property
1347 def n(self):
1348 """The length of this element's column."""
1349 return self.module.n
1351 def __len__(self):
1352 return self.n
1354 def column(self, domain=None):
1355 """
1356 Get a copy of this element's column, optionally converting to a domain.
1357 """
1358 return self.col.convert_to(domain)
1360 @property
1361 def coeffs(self):
1362 return self.col.flat()
1364 @property
1365 def QQ_col(self):
1366 """
1367 :py:class:`~.DomainMatrix` over :ref:`QQ`, equal to
1368 ``self.col / self.denom``, and guaranteed to be dense.
1370 See Also
1371 ========
1373 .Submodule.QQ_matrix
1375 """
1376 if self._QQ_col is None:
1377 self._QQ_col = (self.col / self.denom).to_dense()
1378 return self._QQ_col
1380 def to_parent(self):
1381 """
1382 Transform into a :py:class:`~.ModuleElement` belonging to the parent of
1383 this element's module.
1384 """
1385 if not isinstance(self.module, Submodule):
1386 raise ValueError('Not an element of a Submodule.')
1387 return make_mod_elt(
1388 self.module.parent, self.module.matrix * self.col,
1389 denom=self.module.denom * self.denom)
1391 def to_ancestor(self, anc):
1392 """
1393 Transform into a :py:class:`~.ModuleElement` belonging to a given
1394 ancestor of this element's module.
1396 Parameters
1397 ==========
1399 anc : :py:class:`~.Module`
1401 """
1402 if anc == self.module:
1403 return self
1404 else:
1405 return self.to_parent().to_ancestor(anc)
1407 def over_power_basis(self):
1408 """
1409 Transform into a :py:class:`~.PowerBasisElement` over our
1410 :py:class:`~.PowerBasis` ancestor.
1411 """
1412 e = self
1413 while not isinstance(e.module, PowerBasis):
1414 e = e.to_parent()
1415 return e
1417 def is_compat(self, other):
1418 """
1419 Test whether other is another :py:class:`~.ModuleElement` with same
1420 module.
1421 """
1422 return isinstance(other, ModuleElement) and other.module == self.module
1424 def unify(self, other):
1425 """
1426 Try to make a compatible pair of :py:class:`~.ModuleElement`, one
1427 equivalent to this one, and one equivalent to the other.
1429 Explanation
1430 ===========
1432 We search for the nearest common ancestor module for the pair of
1433 elements, and represent each one there.
1435 Returns
1436 =======
1438 Pair ``(e1, e2)``
1439 Each ``ei`` is a :py:class:`~.ModuleElement`, they belong to the
1440 same :py:class:`~.Module`, ``e1`` is equivalent to ``self``, and
1441 ``e2`` is equivalent to ``other``.
1443 Raises
1444 ======
1446 UnificationFailed
1447 If ``self`` and ``other`` have no common ancestor module.
1449 """
1450 if self.module == other.module:
1451 return self, other
1452 nca = self.module.nearest_common_ancestor(other.module)
1453 if nca is not None:
1454 return self.to_ancestor(nca), other.to_ancestor(nca)
1455 raise UnificationFailed(f"Cannot unify {self} with {other}")
1457 def __eq__(self, other):
1458 if self.is_compat(other):
1459 return self.QQ_col == other.QQ_col
1460 return NotImplemented
1462 def equiv(self, other):
1463 """
1464 A :py:class:`~.ModuleElement` may test as equivalent to a rational
1465 number or another :py:class:`~.ModuleElement`, if they represent the
1466 same algebraic number.
1468 Explanation
1469 ===========
1471 This method is intended to check equivalence only in those cases in
1472 which it is easy to test; namely, when *other* is either a
1473 :py:class:`~.ModuleElement` that can be unified with this one (i.e. one
1474 which shares a common :py:class:`~.PowerBasis` ancestor), or else a
1475 rational number (which is easy because every :py:class:`~.PowerBasis`
1476 represents every rational number).
1478 Parameters
1479 ==========
1481 other : int, :ref:`ZZ`, :ref:`QQ`, :py:class:`~.ModuleElement`
1483 Returns
1484 =======
1486 bool
1488 Raises
1489 ======
1491 UnificationFailed
1492 If ``self`` and ``other`` do not share a common
1493 :py:class:`~.PowerBasis` ancestor.
1495 """
1496 if self == other:
1497 return True
1498 elif isinstance(other, ModuleElement):
1499 a, b = self.unify(other)
1500 return a == b
1501 elif is_rat(other):
1502 if isinstance(self, PowerBasisElement):
1503 return self == self.module(0) * other
1504 else:
1505 return self.over_power_basis().equiv(other)
1506 return False
1508 def __add__(self, other):
1509 """
1510 A :py:class:`~.ModuleElement` can be added to a rational number, or to
1511 another :py:class:`~.ModuleElement`.
1513 Explanation
1514 ===========
1516 When the other summand is a rational number, it will be converted into
1517 a :py:class:`~.ModuleElement` (belonging to the first ancestor of this
1518 module that starts with unity).
1520 In all cases, the sum belongs to the nearest common ancestor (NCA) of
1521 the modules of the two summands. If the NCA does not exist, we return
1522 ``NotImplemented``.
1523 """
1524 if self.is_compat(other):
1525 d, e = self.denom, other.denom
1526 m = ilcm(d, e)
1527 u, v = m // d, m // e
1528 col = to_col([u * a + v * b for a, b in zip(self.coeffs, other.coeffs)])
1529 return type(self)(self.module, col, denom=m).reduced()
1530 elif isinstance(other, ModuleElement):
1531 try:
1532 a, b = self.unify(other)
1533 except UnificationFailed:
1534 return NotImplemented
1535 return a + b
1536 elif is_rat(other):
1537 return self + self.module.element_from_rational(other)
1538 return NotImplemented
1540 __radd__ = __add__
1542 def __neg__(self):
1543 return self * -1
1545 def __sub__(self, other):
1546 return self + (-other)
1548 def __rsub__(self, other):
1549 return -self + other
1551 def __mul__(self, other):
1552 """
1553 A :py:class:`~.ModuleElement` can be multiplied by a rational number,
1554 or by another :py:class:`~.ModuleElement`.
1556 Explanation
1557 ===========
1559 When the multiplier is a rational number, the product is computed by
1560 operating directly on the coefficients of this
1561 :py:class:`~.ModuleElement`.
1563 When the multiplier is another :py:class:`~.ModuleElement`, the product
1564 will belong to the nearest common ancestor (NCA) of the modules of the
1565 two operands, and that NCA must have a multiplication table. If the NCA
1566 does not exist, we return ``NotImplemented``. If the NCA does not have
1567 a mult. table, ``ClosureFailure`` will be raised.
1568 """
1569 if self.is_compat(other):
1570 M = self.module.mult_tab()
1571 A, B = self.col.flat(), other.col.flat()
1572 n = self.n
1573 C = [0] * n
1574 for u in range(n):
1575 for v in range(u, n):
1576 c = A[u] * B[v]
1577 if v > u:
1578 c += A[v] * B[u]
1579 if c != 0:
1580 R = M[u][v]
1581 for k in range(n):
1582 C[k] += c * R[k]
1583 d = self.denom * other.denom
1584 return self.from_int_list(self.module, C, denom=d)
1585 elif isinstance(other, ModuleElement):
1586 try:
1587 a, b = self.unify(other)
1588 except UnificationFailed:
1589 return NotImplemented
1590 return a * b
1591 elif is_rat(other):
1592 a, b = get_num_denom(other)
1593 if a == b == 1:
1594 return self
1595 else:
1596 return make_mod_elt(self.module,
1597 self.col * a, denom=self.denom * b).reduced()
1598 return NotImplemented
1600 __rmul__ = __mul__
1602 def _zeroth_power(self):
1603 return self.module.one()
1605 def _first_power(self):
1606 return self
1608 def __floordiv__(self, a):
1609 if is_rat(a):
1610 a = QQ(a)
1611 return self * (1/a)
1612 elif isinstance(a, ModuleElement):
1613 return self * (1//a)
1614 return NotImplemented
1616 def __rfloordiv__(self, a):
1617 return a // self.over_power_basis()
1619 def __mod__(self, m):
1620 r"""
1621 Reduce this :py:class:`~.ModuleElement` mod a :py:class:`~.Submodule`.
1623 Parameters
1624 ==========
1626 m : int, :ref:`ZZ`, :ref:`QQ`, :py:class:`~.Submodule`
1627 If a :py:class:`~.Submodule`, reduce ``self`` relative to this.
1628 If an integer or rational, reduce relative to the
1629 :py:class:`~.Submodule` that is our own module times this constant.
1631 See Also
1632 ========
1634 .Submodule.reduce_element
1636 """
1637 if is_rat(m):
1638 m = m * self.module.whole_submodule()
1639 if isinstance(m, Submodule) and m.parent == self.module:
1640 return m.reduce_element(self)
1641 return NotImplemented
1644class PowerBasisElement(ModuleElement):
1645 r"""
1646 Subclass for :py:class:`~.ModuleElement` instances whose module is a
1647 :py:class:`~.PowerBasis`.
1648 """
1650 @property
1651 def T(self):
1652 """Access the defining polynomial of the :py:class:`~.PowerBasis`."""
1653 return self.module.T
1655 def numerator(self, x=None):
1656 """Obtain the numerator as a polynomial over :ref:`ZZ`."""
1657 x = x or self.T.gen
1658 return Poly(reversed(self.coeffs), x, domain=ZZ)
1660 def poly(self, x=None):
1661 """Obtain the number as a polynomial over :ref:`QQ`."""
1662 return self.numerator(x=x) // self.denom
1664 @property
1665 def is_rational(self):
1666 """Say whether this element represents a rational number."""
1667 return self.col[1:, :].is_zero_matrix
1669 @property
1670 def generator(self):
1671 """
1672 Return a :py:class:`~.Symbol` to be used when expressing this element
1673 as a polynomial.
1675 If we have an associated :py:class:`~.AlgebraicField` whose primitive
1676 element has an alias symbol, we use that. Otherwise we use the variable
1677 of the minimal polynomial defining the power basis to which we belong.
1678 """
1679 K = self.module.number_field
1680 return K.ext.alias if K and K.ext.is_aliased else self.T.gen
1682 def as_expr(self, x=None):
1683 """Create a Basic expression from ``self``. """
1684 return self.poly(x or self.generator).as_expr()
1686 def norm(self, T=None):
1687 """Compute the norm of this number."""
1688 T = T or self.T
1689 x = T.gen
1690 A = self.numerator(x=x)
1691 return T.resultant(A) // self.denom ** self.n
1693 def inverse(self):
1694 f = self.poly()
1695 f_inv = f.invert(self.T)
1696 return self.module.element_from_poly(f_inv)
1698 def __rfloordiv__(self, a):
1699 return self.inverse() * a
1701 def _negative_power(self, e, modulo=None):
1702 return self.inverse() ** abs(e)
1704 def to_ANP(self):
1705 """Convert to an equivalent :py:class:`~.ANP`. """
1706 return ANP(list(reversed(self.QQ_col.flat())), QQ.map(self.T.rep.rep), QQ)
1708 def to_alg_num(self):
1709 """
1710 Try to convert to an equivalent :py:class:`~.AlgebraicNumber`.
1712 Explanation
1713 ===========
1715 In general, the conversion from an :py:class:`~.AlgebraicNumber` to a
1716 :py:class:`~.PowerBasisElement` throws away information, because an
1717 :py:class:`~.AlgebraicNumber` specifies a complex embedding, while a
1718 :py:class:`~.PowerBasisElement` does not. However, in some cases it is
1719 possible to convert a :py:class:`~.PowerBasisElement` back into an
1720 :py:class:`~.AlgebraicNumber`, namely when the associated
1721 :py:class:`~.PowerBasis` has a reference to an
1722 :py:class:`~.AlgebraicField`.
1724 Returns
1725 =======
1727 :py:class:`~.AlgebraicNumber`
1729 Raises
1730 ======
1732 StructureError
1733 If the :py:class:`~.PowerBasis` to which this element belongs does
1734 not have an associated :py:class:`~.AlgebraicField`.
1736 """
1737 K = self.module.number_field
1738 if K:
1739 return K.to_alg_num(self.to_ANP())
1740 raise StructureError("No associated AlgebraicField")
1743class ModuleHomomorphism:
1744 r"""A homomorphism from one module to another."""
1746 def __init__(self, domain, codomain, mapping):
1747 r"""
1748 Parameters
1749 ==========
1751 domain : :py:class:`~.Module`
1752 The domain of the mapping.
1754 codomain : :py:class:`~.Module`
1755 The codomain of the mapping.
1757 mapping : callable
1758 An arbitrary callable is accepted, but should be chosen so as
1759 to represent an actual module homomorphism. In particular, should
1760 accept elements of *domain* and return elements of *codomain*.
1762 Examples
1763 ========
1765 >>> from sympy import Poly, cyclotomic_poly
1766 >>> from sympy.polys.numberfields.modules import PowerBasis, ModuleHomomorphism
1767 >>> T = Poly(cyclotomic_poly(5))
1768 >>> A = PowerBasis(T)
1769 >>> B = A.submodule_from_gens([2*A(j) for j in range(4)])
1770 >>> phi = ModuleHomomorphism(A, B, lambda x: 6*x)
1771 >>> print(phi.matrix()) # doctest: +SKIP
1772 DomainMatrix([[3, 0, 0, 0], [0, 3, 0, 0], [0, 0, 3, 0], [0, 0, 0, 3]], (4, 4), ZZ)
1774 """
1775 self.domain = domain
1776 self.codomain = codomain
1777 self.mapping = mapping
1779 def matrix(self, modulus=None):
1780 r"""
1781 Compute the matrix of this homomorphism.
1783 Parameters
1784 ==========
1786 modulus : int, optional
1787 A positive prime number $p$ if the matrix should be reduced mod
1788 $p$.
1790 Returns
1791 =======
1793 :py:class:`~.DomainMatrix`
1794 The matrix is over :ref:`ZZ`, or else over :ref:`GF(p)` if a
1795 modulus was given.
1797 """
1798 basis = self.domain.basis_elements()
1799 cols = [self.codomain.represent(self.mapping(elt)) for elt in basis]
1800 if not cols:
1801 return DomainMatrix.zeros((self.codomain.n, 0), ZZ).to_dense()
1802 M = cols[0].hstack(*cols[1:])
1803 if modulus:
1804 M = M.convert_to(FF(modulus))
1805 return M
1807 def kernel(self, modulus=None):
1808 r"""
1809 Compute a Submodule representing the kernel of this homomorphism.
1811 Parameters
1812 ==========
1814 modulus : int, optional
1815 A positive prime number $p$ if the kernel should be computed mod
1816 $p$.
1818 Returns
1819 =======
1821 :py:class:`~.Submodule`
1822 This submodule's generators span the kernel of this
1823 homomorphism over :ref:`ZZ`, or else over :ref:`GF(p)` if a
1824 modulus was given.
1826 """
1827 M = self.matrix(modulus=modulus)
1828 if modulus is None:
1829 M = M.convert_to(QQ)
1830 # Note: Even when working over a finite field, what we want here is
1831 # the pullback into the integers, so in this case the conversion to ZZ
1832 # below is appropriate. When working over ZZ, the kernel should be a
1833 # ZZ-submodule, so, while the conversion to QQ above was required in
1834 # order for the nullspace calculation to work, conversion back to ZZ
1835 # afterward should always work.
1836 # TODO:
1837 # Watch <https://github.com/sympy/sympy/issues/21834>, which calls
1838 # for fraction-free algorithms. If this is implemented, we can skip
1839 # the conversion to `QQ` above.
1840 K = M.nullspace().convert_to(ZZ).transpose()
1841 return self.domain.submodule_from_matrix(K)
1844class ModuleEndomorphism(ModuleHomomorphism):
1845 r"""A homomorphism from one module to itself."""
1847 def __init__(self, domain, mapping):
1848 r"""
1849 Parameters
1850 ==========
1852 domain : :py:class:`~.Module`
1853 The common domain and codomain of the mapping.
1855 mapping : callable
1856 An arbitrary callable is accepted, but should be chosen so as
1857 to represent an actual module endomorphism. In particular, should
1858 accept and return elements of *domain*.
1860 """
1861 super().__init__(domain, domain, mapping)
1864class InnerEndomorphism(ModuleEndomorphism):
1865 r"""
1866 An inner endomorphism on a module, i.e. the endomorphism corresponding to
1867 multiplication by a fixed element.
1868 """
1870 def __init__(self, domain, multiplier):
1871 r"""
1872 Parameters
1873 ==========
1875 domain : :py:class:`~.Module`
1876 The domain and codomain of the endomorphism.
1878 multiplier : :py:class:`~.ModuleElement`
1879 The element $a$ defining the mapping as $x \mapsto a x$.
1881 """
1882 super().__init__(domain, lambda x: multiplier * x)
1883 self.multiplier = multiplier
1886class EndomorphismRing:
1887 r"""The ring of endomorphisms on a module."""
1889 def __init__(self, domain):
1890 """
1891 Parameters
1892 ==========
1894 domain : :py:class:`~.Module`
1895 The domain and codomain of the endomorphisms.
1897 """
1898 self.domain = domain
1900 def inner_endomorphism(self, multiplier):
1901 r"""
1902 Form an inner endomorphism belonging to this endomorphism ring.
1904 Parameters
1905 ==========
1907 multiplier : :py:class:`~.ModuleElement`
1908 Element $a$ defining the inner endomorphism $x \mapsto a x$.
1910 Returns
1911 =======
1913 :py:class:`~.InnerEndomorphism`
1915 """
1916 return InnerEndomorphism(self.domain, multiplier)
1918 def represent(self, element):
1919 r"""
1920 Represent an element of this endomorphism ring, as a single column
1921 vector.
1923 Explanation
1924 ===========
1926 Let $M$ be a module, and $E$ its ring of endomorphisms. Let $N$ be
1927 another module, and consider a homomorphism $\varphi: N \rightarrow E$.
1928 In the event that $\varphi$ is to be represented by a matrix $A$, each
1929 column of $A$ must represent an element of $E$. This is possible when
1930 the elements of $E$ are themselves representable as matrices, by
1931 stacking the columns of such a matrix into a single column.
1933 This method supports calculating such matrices $A$, by representing
1934 an element of this endomorphism ring first as a matrix, and then
1935 stacking that matrix's columns into a single column.
1937 Examples
1938 ========
1940 Note that in these examples we print matrix transposes, to make their
1941 columns easier to inspect.
1943 >>> from sympy import Poly, cyclotomic_poly
1944 >>> from sympy.polys.numberfields.modules import PowerBasis
1945 >>> from sympy.polys.numberfields.modules import ModuleHomomorphism
1946 >>> T = Poly(cyclotomic_poly(5))
1947 >>> M = PowerBasis(T)
1948 >>> E = M.endomorphism_ring()
1950 Let $\zeta$ be a primitive 5th root of unity, a generator of our field,
1951 and consider the inner endomorphism $\tau$ on the ring of integers,
1952 induced by $\zeta$:
1954 >>> zeta = M(1)
1955 >>> tau = E.inner_endomorphism(zeta)
1956 >>> tau.matrix().transpose() # doctest: +SKIP
1957 DomainMatrix(
1958 [[0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1], [-1, -1, -1, -1]],
1959 (4, 4), ZZ)
1961 The matrix representation of $\tau$ is as expected. The first column
1962 shows that multiplying by $\zeta$ carries $1$ to $\zeta$, the second
1963 column that it carries $\zeta$ to $\zeta^2$, and so forth.
1965 The ``represent`` method of the endomorphism ring ``E`` stacks these
1966 into a single column:
1968 >>> E.represent(tau).transpose() # doctest: +SKIP
1969 DomainMatrix(
1970 [[0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, -1, -1, -1, -1]],
1971 (1, 16), ZZ)
1973 This is useful when we want to consider a homomorphism $\varphi$ having
1974 ``E`` as codomain:
1976 >>> phi = ModuleHomomorphism(M, E, lambda x: E.inner_endomorphism(x))
1978 and we want to compute the matrix of such a homomorphism:
1980 >>> phi.matrix().transpose() # doctest: +SKIP
1981 DomainMatrix(
1982 [[1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1],
1983 [0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, -1, -1, -1, -1],
1984 [0, 0, 1, 0, 0, 0, 0, 1, -1, -1, -1, -1, 1, 0, 0, 0],
1985 [0, 0, 0, 1, -1, -1, -1, -1, 1, 0, 0, 0, 0, 1, 0, 0]],
1986 (4, 16), ZZ)
1988 Note that the stacked matrix of $\tau$ occurs as the second column in
1989 this example. This is because $\zeta$ is the second basis element of
1990 ``M``, and $\varphi(\zeta) = \tau$.
1992 Parameters
1993 ==========
1995 element : :py:class:`~.ModuleEndomorphism` belonging to this ring.
1997 Returns
1998 =======
2000 :py:class:`~.DomainMatrix`
2001 Column vector equalling the vertical stacking of all the columns
2002 of the matrix that represents the given *element* as a mapping.
2004 """
2005 if isinstance(element, ModuleEndomorphism) and element.domain == self.domain:
2006 M = element.matrix()
2007 # Transform the matrix into a single column, which should reproduce
2008 # the original columns, one after another.
2009 m, n = M.shape
2010 if n == 0:
2011 return M
2012 return M[:, 0].vstack(*[M[:, j] for j in range(1, n)])
2013 raise NotImplementedError
2016def find_min_poly(alpha, domain, x=None, powers=None):
2017 r"""
2018 Find a polynomial of least degree (not necessarily irreducible) satisfied
2019 by an element of a finitely-generated ring with unity.
2021 Examples
2022 ========
2024 For the $n$th cyclotomic field, $n$ an odd prime, consider the quadratic
2025 equation whose roots are the two periods of length $(n-1)/2$. Article 356
2026 of Gauss tells us that we should get $x^2 + x - (n-1)/4$ or
2027 $x^2 + x + (n+1)/4$ according to whether $n$ is 1 or 3 mod 4, respectively.
2029 >>> from sympy import Poly, cyclotomic_poly, primitive_root, QQ
2030 >>> from sympy.abc import x
2031 >>> from sympy.polys.numberfields.modules import PowerBasis, find_min_poly
2032 >>> n = 13
2033 >>> g = primitive_root(n)
2034 >>> C = PowerBasis(Poly(cyclotomic_poly(n, x)))
2035 >>> ee = [g**(2*k+1) % n for k in range((n-1)//2)]
2036 >>> eta = sum(C(e) for e in ee)
2037 >>> print(find_min_poly(eta, QQ, x=x).as_expr())
2038 x**2 + x - 3
2039 >>> n = 19
2040 >>> g = primitive_root(n)
2041 >>> C = PowerBasis(Poly(cyclotomic_poly(n, x)))
2042 >>> ee = [g**(2*k+2) % n for k in range((n-1)//2)]
2043 >>> eta = sum(C(e) for e in ee)
2044 >>> print(find_min_poly(eta, QQ, x=x).as_expr())
2045 x**2 + x + 5
2047 Parameters
2048 ==========
2050 alpha : :py:class:`~.ModuleElement`
2051 The element whose min poly is to be found, and whose module has
2052 multiplication and starts with unity.
2054 domain : :py:class:`~.Domain`
2055 The desired domain of the polynomial.
2057 x : :py:class:`~.Symbol`, optional
2058 The desired variable for the polynomial.
2060 powers : list, optional
2061 If desired, pass an empty list. The powers of *alpha* (as
2062 :py:class:`~.ModuleElement` instances) from the zeroth up to the degree
2063 of the min poly will be recorded here, as we compute them.
2065 Returns
2066 =======
2068 :py:class:`~.Poly`, ``None``
2069 The minimal polynomial for alpha, or ``None`` if no polynomial could be
2070 found over the desired domain.
2072 Raises
2073 ======
2075 MissingUnityError
2076 If the module to which alpha belongs does not start with unity.
2077 ClosureFailure
2078 If the module to which alpha belongs is not closed under
2079 multiplication.
2081 """
2082 R = alpha.module
2083 if not R.starts_with_unity():
2084 raise MissingUnityError("alpha must belong to finitely generated ring with unity.")
2085 if powers is None:
2086 powers = []
2087 one = R(0)
2088 powers.append(one)
2089 powers_matrix = one.column(domain=domain)
2090 ak = alpha
2091 m = None
2092 for k in range(1, R.n + 1):
2093 powers.append(ak)
2094 ak_col = ak.column(domain=domain)
2095 try:
2096 X = powers_matrix._solve(ak_col)[0]
2097 except DMBadInputError:
2098 # This means alpha^k still isn't in the domain-span of the lower powers.
2099 powers_matrix = powers_matrix.hstack(ak_col)
2100 ak *= alpha
2101 else:
2102 # alpha^k is in the domain-span of the lower powers, so we have found a
2103 # minimal-degree poly for alpha.
2104 coeffs = [1] + [-c for c in reversed(X.to_list_flat())]
2105 x = x or Dummy('x')
2106 if domain.is_FF:
2107 m = Poly(coeffs, x, modulus=domain.mod)
2108 else:
2109 m = Poly(coeffs, x, domain=domain)
2110 break
2111 return m