Coverage for /usr/lib/python3/dist-packages/sympy/polys/numberfields/utilities.py: 22%
128 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1"""Utilities for algebraic number theory. """
3from sympy.core.sympify import sympify
4from sympy.ntheory.factor_ import factorint
5from sympy.polys.domains.rationalfield import QQ
6from sympy.polys.domains.integerring import ZZ
7from sympy.polys.matrices.exceptions import DMRankError
8from sympy.polys.numberfields.minpoly import minpoly
9from sympy.printing.lambdarepr import IntervalPrinter
10from sympy.utilities.decorator import public
11from sympy.utilities.lambdify import lambdify
13from mpmath import mp
16def is_rat(c):
17 r"""
18 Test whether an argument is of an acceptable type to be used as a rational
19 number.
21 Explanation
22 ===========
24 Returns ``True`` on any argument of type ``int``, :ref:`ZZ`, or :ref:`QQ`.
26 See Also
27 ========
29 is_int
31 """
32 # ``c in QQ`` is too accepting (e.g. ``3.14 in QQ`` is ``True``),
33 # ``QQ.of_type(c)`` is too demanding (e.g. ``QQ.of_type(3)`` is ``False``).
34 #
35 # Meanwhile, if gmpy2 is installed then ``ZZ.of_type()`` accepts only
36 # ``mpz``, not ``int``, so we need another clause to ensure ``int`` is
37 # accepted.
38 return isinstance(c, int) or ZZ.of_type(c) or QQ.of_type(c)
41def is_int(c):
42 r"""
43 Test whether an argument is of an acceptable type to be used as an integer.
45 Explanation
46 ===========
48 Returns ``True`` on any argument of type ``int`` or :ref:`ZZ`.
50 See Also
51 ========
53 is_rat
55 """
56 # If gmpy2 is installed then ``ZZ.of_type()`` accepts only
57 # ``mpz``, not ``int``, so we need another clause to ensure ``int`` is
58 # accepted.
59 return isinstance(c, int) or ZZ.of_type(c)
62def get_num_denom(c):
63 r"""
64 Given any argument on which :py:func:`~.is_rat` is ``True``, return the
65 numerator and denominator of this number.
67 See Also
68 ========
70 is_rat
72 """
73 r = QQ(c)
74 return r.numerator, r.denominator
77@public
78def extract_fundamental_discriminant(a):
79 r"""
80 Extract a fundamental discriminant from an integer *a*.
82 Explanation
83 ===========
85 Given any rational integer *a* that is 0 or 1 mod 4, write $a = d f^2$,
86 where $d$ is either 1 or a fundamental discriminant, and return a pair
87 of dictionaries ``(D, F)`` giving the prime factorizations of $d$ and $f$
88 respectively, in the same format returned by :py:func:`~.factorint`.
90 A fundamental discriminant $d$ is different from unity, and is either
91 1 mod 4 and squarefree, or is 0 mod 4 and such that $d/4$ is squarefree
92 and 2 or 3 mod 4. This is the same as being the discriminant of some
93 quadratic field.
95 Examples
96 ========
98 >>> from sympy.polys.numberfields.utilities import extract_fundamental_discriminant
99 >>> print(extract_fundamental_discriminant(-432))
100 ({3: 1, -1: 1}, {2: 2, 3: 1})
102 For comparison:
104 >>> from sympy import factorint
105 >>> print(factorint(-432))
106 {2: 4, 3: 3, -1: 1}
108 Parameters
109 ==========
111 a: int, must be 0 or 1 mod 4
113 Returns
114 =======
116 Pair ``(D, F)`` of dictionaries.
118 Raises
119 ======
121 ValueError
122 If *a* is not 0 or 1 mod 4.
124 References
125 ==========
127 .. [1] Cohen, H. *A Course in Computational Algebraic Number Theory.*
128 (See Prop. 5.1.3)
130 """
131 if a % 4 not in [0, 1]:
132 raise ValueError('To extract fundamental discriminant, number must be 0 or 1 mod 4.')
133 if a == 0:
134 return {}, {0: 1}
135 if a == 1:
136 return {}, {}
137 a_factors = factorint(a)
138 D = {}
139 F = {}
140 # First pass: just make d squarefree, and a/d a perfect square.
141 # We'll count primes (and units! i.e. -1) that are 3 mod 4 and present in d.
142 num_3_mod_4 = 0
143 for p, e in a_factors.items():
144 if e % 2 == 1:
145 D[p] = 1
146 if p % 4 == 3:
147 num_3_mod_4 += 1
148 if e >= 3:
149 F[p] = (e - 1) // 2
150 else:
151 F[p] = e // 2
152 # Second pass: if d is cong. to 2 or 3 mod 4, then we must steal away
153 # another factor of 4 from f**2 and give it to d.
154 even = 2 in D
155 if even or num_3_mod_4 % 2 == 1:
156 e2 = F[2]
157 assert e2 > 0
158 if e2 == 1:
159 del F[2]
160 else:
161 F[2] = e2 - 1
162 D[2] = 3 if even else 2
163 return D, F
166@public
167class AlgIntPowers:
168 r"""
169 Compute the powers of an algebraic integer.
171 Explanation
172 ===========
174 Given an algebraic integer $\theta$ by its monic irreducible polynomial
175 ``T`` over :ref:`ZZ`, this class computes representations of arbitrarily
176 high powers of $\theta$, as :ref:`ZZ`-linear combinations over
177 $\{1, \theta, \ldots, \theta^{n-1}\}$, where $n = \deg(T)$.
179 The representations are computed using the linear recurrence relations for
180 powers of $\theta$, derived from the polynomial ``T``. See [1], Sec. 4.2.2.
182 Optionally, the representations may be reduced with respect to a modulus.
184 Examples
185 ========
187 >>> from sympy import Poly, cyclotomic_poly
188 >>> from sympy.polys.numberfields.utilities import AlgIntPowers
189 >>> T = Poly(cyclotomic_poly(5))
190 >>> zeta_pow = AlgIntPowers(T)
191 >>> print(zeta_pow[0])
192 [1, 0, 0, 0]
193 >>> print(zeta_pow[1])
194 [0, 1, 0, 0]
195 >>> print(zeta_pow[4]) # doctest: +SKIP
196 [-1, -1, -1, -1]
197 >>> print(zeta_pow[24]) # doctest: +SKIP
198 [-1, -1, -1, -1]
200 References
201 ==========
203 .. [1] Cohen, H. *A Course in Computational Algebraic Number Theory.*
205 """
207 def __init__(self, T, modulus=None):
208 """
209 Parameters
210 ==========
212 T : :py:class:`~.Poly`
213 The monic irreducible polynomial over :ref:`ZZ` defining the
214 algebraic integer.
216 modulus : int, None, optional
217 If not ``None``, all representations will be reduced w.r.t. this.
219 """
220 self.T = T
221 self.modulus = modulus
222 self.n = T.degree()
223 self.powers_n_and_up = [[-c % self for c in reversed(T.rep.rep)][:-1]]
224 self.max_so_far = self.n
226 def red(self, exp):
227 return exp if self.modulus is None else exp % self.modulus
229 def __rmod__(self, other):
230 return self.red(other)
232 def compute_up_through(self, e):
233 m = self.max_so_far
234 if e <= m: return
235 n = self.n
236 r = self.powers_n_and_up
237 c = r[0]
238 for k in range(m+1, e+1):
239 b = r[k-1-n][n-1]
240 r.append(
241 [c[0]*b % self] + [
242 (r[k-1-n][i-1] + c[i]*b) % self for i in range(1, n)
243 ]
244 )
245 self.max_so_far = e
247 def get(self, e):
248 n = self.n
249 if e < 0:
250 raise ValueError('Exponent must be non-negative.')
251 elif e < n:
252 return [1 if i == e else 0 for i in range(n)]
253 else:
254 self.compute_up_through(e)
255 return self.powers_n_and_up[e - n]
257 def __getitem__(self, item):
258 return self.get(item)
261@public
262def coeff_search(m, R):
263 r"""
264 Generate coefficients for searching through polynomials.
266 Explanation
267 ===========
269 Lead coeff is always non-negative. Explore all combinations with coeffs
270 bounded in absolute value before increasing the bound. Skip the all-zero
271 list, and skip any repeats. See examples.
273 Examples
274 ========
276 >>> from sympy.polys.numberfields.utilities import coeff_search
277 >>> cs = coeff_search(2, 1)
278 >>> C = [next(cs) for i in range(13)]
279 >>> print(C)
280 [[1, 1], [1, 0], [1, -1], [0, 1], [2, 2], [2, 1], [2, 0], [2, -1], [2, -2],
281 [1, 2], [1, -2], [0, 2], [3, 3]]
283 Parameters
284 ==========
286 m : int
287 Length of coeff list.
288 R : int
289 Initial max abs val for coeffs (will increase as search proceeds).
291 Returns
292 =======
294 generator
295 Infinite generator of lists of coefficients.
297 """
298 R0 = R
299 c = [R] * m
300 while True:
301 if R == R0 or R in c or -R in c:
302 yield c[:]
303 j = m - 1
304 while c[j] == -R:
305 j -= 1
306 c[j] -= 1
307 for i in range(j + 1, m):
308 c[i] = R
309 for j in range(m):
310 if c[j] != 0:
311 break
312 else:
313 R += 1
314 c = [R] * m
317def supplement_a_subspace(M):
318 r"""
319 Extend a basis for a subspace to a basis for the whole space.
321 Explanation
322 ===========
324 Given an $n \times r$ matrix *M* of rank $r$ (so $r \leq n$), this function
325 computes an invertible $n \times n$ matrix $B$ such that the first $r$
326 columns of $B$ equal *M*.
328 This operation can be interpreted as a way of extending a basis for a
329 subspace, to give a basis for the whole space.
331 To be precise, suppose you have an $n$-dimensional vector space $V$, with
332 basis $\{v_1, v_2, \ldots, v_n\}$, and an $r$-dimensional subspace $W$ of
333 $V$, spanned by a basis $\{w_1, w_2, \ldots, w_r\}$, where the $w_j$ are
334 given as linear combinations of the $v_i$. If the columns of *M* represent
335 the $w_j$ as such linear combinations, then the columns of the matrix $B$
336 computed by this function give a new basis $\{u_1, u_2, \ldots, u_n\}$ for
337 $V$, again relative to the $\{v_i\}$ basis, and such that $u_j = w_j$
338 for $1 \leq j \leq r$.
340 Examples
341 ========
343 Note: The function works in terms of columns, so in these examples we
344 print matrix transposes in order to make the columns easier to inspect.
346 >>> from sympy.polys.matrices import DM
347 >>> from sympy import QQ, FF
348 >>> from sympy.polys.numberfields.utilities import supplement_a_subspace
349 >>> M = DM([[1, 7, 0], [2, 3, 4]], QQ).transpose()
350 >>> print(supplement_a_subspace(M).to_Matrix().transpose())
351 Matrix([[1, 7, 0], [2, 3, 4], [1, 0, 0]])
353 >>> M2 = M.convert_to(FF(7))
354 >>> print(M2.to_Matrix().transpose())
355 Matrix([[1, 0, 0], [2, 3, -3]])
356 >>> print(supplement_a_subspace(M2).to_Matrix().transpose())
357 Matrix([[1, 0, 0], [2, 3, -3], [0, 1, 0]])
359 Parameters
360 ==========
362 M : :py:class:`~.DomainMatrix`
363 The columns give the basis for the subspace.
365 Returns
366 =======
368 :py:class:`~.DomainMatrix`
369 This matrix is invertible and its first $r$ columns equal *M*.
371 Raises
372 ======
374 DMRankError
375 If *M* was not of maximal rank.
377 References
378 ==========
380 .. [1] Cohen, H. *A Course in Computational Algebraic Number Theory*
381 (See Sec. 2.3.2.)
383 """
384 n, r = M.shape
385 # Let In be the n x n identity matrix.
386 # Form the augmented matrix [M | In] and compute RREF.
387 Maug = M.hstack(M.eye(n, M.domain))
388 R, pivots = Maug.rref()
389 if pivots[:r] != tuple(range(r)):
390 raise DMRankError('M was not of maximal rank')
391 # Let J be the n x r matrix equal to the first r columns of In.
392 # Since M is of rank r, RREF reduces [M | In] to [J | A], where A is the product of
393 # elementary matrices Ei corresp. to the row ops performed by RREF. Since the Ei are
394 # invertible, so is A. Let B = A^(-1).
395 A = R[:, r:]
396 B = A.inv()
397 # Then B is the desired matrix. It is invertible, since B^(-1) == A.
398 # And A * [M | In] == [J | A]
399 # => A * M == J
400 # => M == B * J == the first r columns of B.
401 return B
404@public
405def isolate(alg, eps=None, fast=False):
406 """
407 Find a rational isolating interval for a real algebraic number.
409 Examples
410 ========
412 >>> from sympy import isolate, sqrt, Rational
413 >>> print(isolate(sqrt(2))) # doctest: +SKIP
414 (1, 2)
415 >>> print(isolate(sqrt(2), eps=Rational(1, 100)))
416 (24/17, 17/12)
418 Parameters
419 ==========
421 alg : str, int, :py:class:`~.Expr`
422 The algebraic number to be isolated. Must be a real number, to use this
423 particular function. However, see also :py:meth:`.Poly.intervals`,
424 which isolates complex roots when you pass ``all=True``.
425 eps : positive element of :ref:`QQ`, None, optional (default=None)
426 Precision to be passed to :py:meth:`.Poly.refine_root`
427 fast : boolean, optional (default=False)
428 Say whether fast refinement procedure should be used.
429 (Will be passed to :py:meth:`.Poly.refine_root`.)
431 Returns
432 =======
434 Pair of rational numbers defining an isolating interval for the given
435 algebraic number.
437 See Also
438 ========
440 .Poly.intervals
442 """
443 alg = sympify(alg)
445 if alg.is_Rational:
446 return (alg, alg)
447 elif not alg.is_real:
448 raise NotImplementedError(
449 "complex algebraic numbers are not supported")
451 func = lambdify((), alg, modules="mpmath", printer=IntervalPrinter())
453 poly = minpoly(alg, polys=True)
454 intervals = poly.intervals(sqf=True)
456 dps, done = mp.dps, False
458 try:
459 while not done:
460 alg = func()
462 for a, b in intervals:
463 if a <= alg.a and alg.b <= b:
464 done = True
465 break
466 else:
467 mp.dps *= 2
468 finally:
469 mp.dps = dps
471 if eps is not None:
472 a, b = poly.refine_root(a, b, eps=eps, fast=fast)
474 return (a, b)