galois_resolvents.py 24 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673
  1. r"""
  2. Galois resolvents
  3. Each of the functions in ``sympy.polys.numberfields.galoisgroups`` that
  4. computes Galois groups for a particular degree $n$ uses resolvents. Given the
  5. polynomial $T$ whose Galois group is to be computed, a resolvent is a
  6. polynomial $R$ whose roots are defined as functions of the roots of $T$.
  7. One way to compute the coefficients of $R$ is by approximating the roots of $T$
  8. to sufficient precision. This module defines a :py:class:`~.Resolvent` class
  9. that handles this job, determining the necessary precision, and computing $R$.
  10. In some cases, the coefficients of $R$ are symmetric in the roots of $T$,
  11. meaning they are equal to fixed functions of the coefficients of $T$. Therefore
  12. another approach is to compute these functions once and for all, and record
  13. them in a lookup table. This module defines code that can compute such tables.
  14. The tables for polynomials $T$ of degrees 4 through 6, produced by this code,
  15. are recorded in the resolvent_lookup.py module.
  16. """
  17. from sympy.core.evalf import (
  18. evalf, fastlog, _evalf_with_bounded_error, quad_to_mpmath,
  19. )
  20. from sympy.core.symbol import symbols, Dummy
  21. from sympy.polys.densetools import dup_eval
  22. from sympy.polys.domains import ZZ
  23. from sympy.polys.numberfields.resolvent_lookup import resolvent_coeff_lambdas
  24. from sympy.polys.orderings import lex
  25. from sympy.polys.polyroots import preprocess_roots
  26. from sympy.polys.polytools import Poly
  27. from sympy.polys.rings import xring
  28. from sympy.polys.specialpolys import symmetric_poly
  29. from sympy.utilities.lambdify import lambdify
  30. from mpmath import MPContext
  31. from mpmath.libmp.libmpf import prec_to_dps
  32. class GaloisGroupException(Exception):
  33. ...
  34. class ResolventException(GaloisGroupException):
  35. ...
  36. class Resolvent:
  37. r"""
  38. If $G$ is a subgroup of the symmetric group $S_n$,
  39. $F$ a multivariate polynomial in $\mathbb{Z}[X_1, \ldots, X_n]$,
  40. $H$ the stabilizer of $F$ in $G$ (i.e. the permutations $\sigma$ such that
  41. $F(X_{\sigma(1)}, \ldots, X_{\sigma(n)}) = F(X_1, \ldots, X_n)$), and $s$
  42. a set of left coset representatives of $H$ in $G$, then the resolvent
  43. polynomial $R(Y)$ is the product over $\sigma \in s$ of
  44. $Y - F(X_{\sigma(1)}, \ldots, X_{\sigma(n)})$.
  45. For example, consider the resolvent for the form
  46. $$F = X_0 X_2 + X_1 X_3$$
  47. and the group $G = S_4$. In this case, the stabilizer $H$ is the dihedral
  48. group $D4 = < (0123), (02) >$, and a set of representatives of $G/H$ is
  49. $\{I, (01), (03)\}$. The resolvent can be constructed as follows:
  50. >>> from sympy.combinatorics.permutations import Permutation
  51. >>> from sympy.core.symbol import symbols
  52. >>> from sympy.polys.numberfields.galoisgroups import Resolvent
  53. >>> X = symbols('X0 X1 X2 X3')
  54. >>> F = X[0]*X[2] + X[1]*X[3]
  55. >>> s = [Permutation([0, 1, 2, 3]), Permutation([1, 0, 2, 3]),
  56. ... Permutation([3, 1, 2, 0])]
  57. >>> R = Resolvent(F, X, s)
  58. This resolvent has three roots, which are the conjugates of ``F`` under the
  59. three permutations in ``s``:
  60. >>> R.root_lambdas[0](*X)
  61. X0*X2 + X1*X3
  62. >>> R.root_lambdas[1](*X)
  63. X0*X3 + X1*X2
  64. >>> R.root_lambdas[2](*X)
  65. X0*X1 + X2*X3
  66. Resolvents are useful for computing Galois groups. Given a polynomial $T$
  67. of degree $n$, we will use a resolvent $R$ where $Gal(T) \leq G \leq S_n$.
  68. We will then want to substitute the roots of $T$ for the variables $X_i$
  69. in $R$, and study things like the discriminant of $R$, and the way $R$
  70. factors over $\mathbb{Q}$.
  71. From the symmetry in $R$'s construction, and since $Gal(T) \leq G$, we know
  72. from Galois theory that the coefficients of $R$ must lie in $\mathbb{Z}$.
  73. This allows us to compute the coefficients of $R$ by approximating the
  74. roots of $T$ to sufficient precision, plugging these values in for the
  75. variables $X_i$ in the coefficient expressions of $R$, and then simply
  76. rounding to the nearest integer.
  77. In order to determine a sufficient precision for the roots of $T$, this
  78. ``Resolvent`` class imposes certain requirements on the form ``F``. It
  79. could be possible to design a different ``Resolvent`` class, that made
  80. different precision estimates, and different assumptions about ``F``.
  81. ``F`` must be homogeneous, and all terms must have unit coefficient.
  82. Furthermore, if $r$ is the number of terms in ``F``, and $t$ the total
  83. degree, and if $m$ is the number of conjugates of ``F``, i.e. the number
  84. of permutations in ``s``, then we require that $m < r 2^t$. Again, it is
  85. not impossible to work with forms ``F`` that violate these assumptions, but
  86. this ``Resolvent`` class requires them.
  87. Since determining the integer coefficients of the resolvent for a given
  88. polynomial $T$ is one of the main problems this class solves, we take some
  89. time to explain the precision bounds it uses.
  90. The general problem is:
  91. Given a multivariate polynomial $P \in \mathbb{Z}[X_1, \ldots, X_n]$, and a
  92. bound $M \in \mathbb{R}_+$, compute an $\varepsilon > 0$ such that for any
  93. complex numbers $a_1, \ldots, a_n$ with $|a_i| < M$, if the $a_i$ are
  94. approximated to within an accuracy of $\varepsilon$ by $b_i$, that is,
  95. $|a_i - b_i| < \varepsilon$ for $i = 1, \ldots, n$, then
  96. $|P(a_1, \ldots, a_n) - P(b_1, \ldots, b_n)| < 1/2$. In other words, if it
  97. is known that $P(a_1, \ldots, a_n) = c$ for some $c \in \mathbb{Z}$, then
  98. $P(b_1, \ldots, b_n)$ can be rounded to the nearest integer in order to
  99. determine $c$.
  100. To derive our error bound, consider the monomial $xyz$. Defining
  101. $d_i = b_i - a_i$, our error is
  102. $|(a_1 + d_1)(a_2 + d_2)(a_3 + d_3) - a_1 a_2 a_3|$, which is bounded
  103. above by $|(M + \varepsilon)^3 - M^3|$. Passing to a general monomial of
  104. total degree $t$, this expression is bounded by
  105. $M^{t-1}\varepsilon(t + 2^t\varepsilon/M)$ provided $\varepsilon < M$,
  106. and by $(t+1)M^{t-1}\varepsilon$ provided $\varepsilon < M/2^t$.
  107. But since our goal is to make the error less than $1/2$, we will choose
  108. $\varepsilon < 1/(2(t+1)M^{t-1})$, which implies the condition that
  109. $\varepsilon < M/2^t$, as long as $M \geq 2$.
  110. Passing from the general monomial to the general polynomial is easy, by
  111. scaling and summing error bounds.
  112. In our specific case, we are given a homogeneous polynomial $F$ of
  113. $r$ terms and total degree $t$, all of whose coefficients are $\pm 1$. We
  114. are given the $m$ permutations that make the conjugates of $F$, and
  115. we want to bound the error in the coefficients of the monic polynomial
  116. $R(Y)$ having $F$ and its conjugates as roots (i.e. the resolvent).
  117. For $j$ from $1$ to $m$, the coefficient of $Y^{m-j}$ in $R(Y)$ is the
  118. $j$th elementary symmetric polynomial in the conjugates of $F$. This sums
  119. the products of these conjugates, taken $j$ at a time, in all possible
  120. combinations. There are $\binom{m}{j}$ such combinations, and each product
  121. of $j$ conjugates of $F$ expands to a sum of $r^j$ terms, each of unit
  122. coefficient, and total degree $jt$. An error bound for the $j$th coeff of
  123. $R$ is therefore
  124. $$\binom{m}{j} r^j (jt + 1) M^{jt - 1} \varepsilon$$
  125. When our goal is to evaluate all the coefficients of $R$, we will want to
  126. use the maximum of these error bounds. It is clear that this bound is
  127. strictly increasing for $j$ up to the ceiling of $m/2$. After that point,
  128. the first factor $\binom{m}{j}$ begins to decrease, while the others
  129. continue to increase. However, the binomial coefficient never falls by more
  130. than a factor of $1/m$ at a time, so our assumptions that $M \geq 2$ and
  131. $m < r 2^t$ are enough to tell us that the constant coefficient of $R$,
  132. i.e. that where $j = m$, has the largest error bound. Therefore we can use
  133. $$r^m (mt + 1) M^{mt - 1} \varepsilon$$
  134. as our error bound for all the coefficients.
  135. Note that this bound is also (more than) adequate to determine whether any
  136. of the roots of $R$ is an integer. Each of these roots is a single
  137. conjugate of $F$, which contains less error than the trace, i.e. the
  138. coefficient of $Y^{m - 1}$. By rounding the roots of $R$ to the nearest
  139. integers, we therefore get all the candidates for integer roots of $R$. By
  140. plugging these candidates into $R$, we can check whether any of them
  141. actually is a root.
  142. Note: We take the definition of resolvent from Cohen, but the error bound
  143. is ours.
  144. References
  145. ==========
  146. .. [1] Cohen, H. *A Course in Computational Algebraic Number Theory*.
  147. (Def 6.3.2)
  148. """
  149. def __init__(self, F, X, s):
  150. r"""
  151. Parameters
  152. ==========
  153. F : :py:class:`~.Expr`
  154. polynomial in the symbols in *X*
  155. X : list of :py:class:`~.Symbol`
  156. s : list of :py:class:`~.Permutation`
  157. representing the cosets of the stabilizer of *F* in
  158. some subgroup $G$ of $S_n$, where $n$ is the length of *X*.
  159. """
  160. self.F = F
  161. self.X = X
  162. self.s = s
  163. # Number of conjugates:
  164. self.m = len(s)
  165. # Total degree of F (computed below):
  166. self.t = None
  167. # Number of terms in F (computed below):
  168. self.r = 0
  169. for monom, coeff in Poly(F).terms():
  170. if abs(coeff) != 1:
  171. raise ResolventException('Resolvent class expects forms with unit coeffs')
  172. t = sum(monom)
  173. if t != self.t and self.t is not None:
  174. raise ResolventException('Resolvent class expects homogeneous forms')
  175. self.t = t
  176. self.r += 1
  177. m, t, r = self.m, self.t, self.r
  178. if not m < r * 2**t:
  179. raise ResolventException('Resolvent class expects m < r*2^t')
  180. M = symbols('M')
  181. # Precision sufficient for computing the coeffs of the resolvent:
  182. self.coeff_prec_func = Poly(r**m*(m*t + 1)*M**(m*t - 1))
  183. # Precision sufficient for checking whether any of the roots of the
  184. # resolvent are integers:
  185. self.root_prec_func = Poly(r*(t + 1)*M**(t - 1))
  186. # The conjugates of F are the roots of the resolvent.
  187. # For evaluating these to required numerical precisions, we need
  188. # lambdified versions.
  189. # Note: for a given permutation sigma, the conjugate (sigma F) is
  190. # equivalent to lambda [sigma^(-1) X]: F.
  191. self.root_lambdas = [
  192. lambdify((~s[j])(X), F)
  193. for j in range(self.m)
  194. ]
  195. # For evaluating the coeffs, we'll also need lambdified versions of
  196. # the elementary symmetric functions for degree m.
  197. Y = symbols('Y')
  198. R = symbols(' '.join(f'R{i}' for i in range(m)))
  199. f = 1
  200. for r in R:
  201. f *= (Y - r)
  202. C = Poly(f, Y).coeffs()
  203. self.esf_lambdas = [lambdify(R, c) for c in C]
  204. def get_prec(self, M, target='coeffs'):
  205. r"""
  206. For a given upper bound *M* on the magnitude of the complex numbers to
  207. be plugged in for this resolvent's symbols, compute a sufficient
  208. precision for evaluating those complex numbers, such that the
  209. coefficients, or the integer roots, of the resolvent can be determined.
  210. Parameters
  211. ==========
  212. M : real number
  213. Upper bound on magnitude of the complex numbers to be plugged in.
  214. target : str, 'coeffs' or 'roots', default='coeffs'
  215. Name the task for which a sufficient precision is desired.
  216. This is either determining the coefficients of the resolvent
  217. ('coeffs') or determining its possible integer roots ('roots').
  218. The latter may require significantly lower precision.
  219. Returns
  220. =======
  221. int $m$
  222. such that $2^{-m}$ is a sufficient upper bound on the
  223. error in approximating the complex numbers to be plugged in.
  224. """
  225. # As explained in the docstring for this class, our precision estimates
  226. # require that M be at least 2.
  227. M = max(M, 2)
  228. f = self.coeff_prec_func if target == 'coeffs' else self.root_prec_func
  229. r, _, _, _ = evalf(2*f(M), 1, {})
  230. return fastlog(r) + 1
  231. def approximate_roots_of_poly(self, T, target='coeffs'):
  232. """
  233. Approximate the roots of a given polynomial *T* to sufficient precision
  234. in order to evaluate this resolvent's coefficients, or determine
  235. whether the resolvent has an integer root.
  236. Parameters
  237. ==========
  238. T : :py:class:`~.Poly`
  239. target : str, 'coeffs' or 'roots', default='coeffs'
  240. Set the approximation precision to be sufficient for the desired
  241. task, which is either determining the coefficients of the resolvent
  242. ('coeffs') or determining its possible integer roots ('roots').
  243. The latter may require significantly lower precision.
  244. Returns
  245. =======
  246. list of elements of :ref:`CC`
  247. """
  248. ctx = MPContext()
  249. # Because sympy.polys.polyroots._integer_basis() is called when a CRootOf
  250. # is formed, we proactively extract the integer basis now. This means that
  251. # when we call T.all_roots(), every root will be a CRootOf, not a Mul
  252. # of Integer*CRootOf.
  253. coeff, T = preprocess_roots(T)
  254. coeff = ctx.mpf(str(coeff))
  255. scaled_roots = T.all_roots(radicals=False)
  256. # Since we're going to be approximating the roots of T anyway, we can
  257. # get a good upper bound on the magnitude of the roots by starting with
  258. # a very low precision approx.
  259. approx0 = [coeff * quad_to_mpmath(_evalf_with_bounded_error(r, m=0)) for r in scaled_roots]
  260. # Here we add 1 to account for the possible error in our initial approximation.
  261. M = max(abs(b) for b in approx0) + 1
  262. m = self.get_prec(M, target=target)
  263. n = fastlog(M._mpf_) + 1
  264. p = m + n + 1
  265. ctx.prec = p
  266. d = prec_to_dps(p)
  267. approx1 = [r.eval_approx(d, return_mpmath=True) for r in scaled_roots]
  268. approx1 = [coeff*ctx.mpc(r) for r in approx1]
  269. return approx1
  270. @staticmethod
  271. def round_mpf(a):
  272. if isinstance(a, int):
  273. return a
  274. # If we use python's built-in `round()`, we lose precision.
  275. # If we use `ZZ` directly, we may add or subtract 1.
  276. return ZZ(a.context.nint(a))
  277. def round_roots_to_integers_for_poly(self, T):
  278. """
  279. For a given polynomial *T*, round the roots of this resolvent to the
  280. nearest integers.
  281. Explanation
  282. ===========
  283. None of the integers returned by this method is guaranteed to be a
  284. root of the resolvent; however, if the resolvent has any integer roots
  285. (for the given polynomial *T*), then they must be among these.
  286. If the coefficients of the resolvent are also desired, then this method
  287. should not be used. Instead, use the ``eval_for_poly`` method. This
  288. method may be significantly faster than ``eval_for_poly``.
  289. Parameters
  290. ==========
  291. T : :py:class:`~.Poly`
  292. Returns
  293. =======
  294. dict
  295. Keys are the indices of those permutations in ``self.s`` such that
  296. the corresponding root did round to a rational integer.
  297. Values are :ref:`ZZ`.
  298. """
  299. approx_roots_of_T = self.approximate_roots_of_poly(T, target='roots')
  300. approx_roots_of_self = [r(*approx_roots_of_T) for r in self.root_lambdas]
  301. return {
  302. i: self.round_mpf(r.real)
  303. for i, r in enumerate(approx_roots_of_self)
  304. if self.round_mpf(r.imag) == 0
  305. }
  306. def eval_for_poly(self, T, find_integer_root=False):
  307. r"""
  308. Compute the integer values of the coefficients of this resolvent, when
  309. plugging in the roots of a given polynomial.
  310. Parameters
  311. ==========
  312. T : :py:class:`~.Poly`
  313. find_integer_root : ``bool``, default ``False``
  314. If ``True``, then also determine whether the resolvent has an
  315. integer root, and return the first one found, along with its
  316. index, i.e. the index of the permutation ``self.s[i]`` it
  317. corresponds to.
  318. Returns
  319. =======
  320. Tuple ``(R, a, i)``
  321. ``R`` is this resolvent as a dense univariate polynomial over
  322. :ref:`ZZ`, i.e. a list of :ref:`ZZ`.
  323. If *find_integer_root* was ``True``, then ``a`` and ``i`` are the
  324. first integer root found, and its index, if one exists.
  325. Otherwise ``a`` and ``i`` are both ``None``.
  326. """
  327. approx_roots_of_T = self.approximate_roots_of_poly(T, target='coeffs')
  328. approx_roots_of_self = [r(*approx_roots_of_T) for r in self.root_lambdas]
  329. approx_coeffs_of_self = [c(*approx_roots_of_self) for c in self.esf_lambdas]
  330. R = []
  331. for c in approx_coeffs_of_self:
  332. if self.round_mpf(c.imag) != 0:
  333. # If precision was enough, this should never happen.
  334. raise ResolventException(f"Got non-integer coeff for resolvent: {c}")
  335. R.append(self.round_mpf(c.real))
  336. a0, i0 = None, None
  337. if find_integer_root:
  338. for i, r in enumerate(approx_roots_of_self):
  339. if self.round_mpf(r.imag) != 0:
  340. continue
  341. if not dup_eval(R, (a := self.round_mpf(r.real)), ZZ):
  342. a0, i0 = a, i
  343. break
  344. return R, a0, i0
  345. def wrap(text, width=80):
  346. """Line wrap a polynomial expression. """
  347. out = ''
  348. col = 0
  349. for c in text:
  350. if c == ' ' and col > width:
  351. c, col = '\n', 0
  352. else:
  353. col += 1
  354. out += c
  355. return out
  356. def s_vars(n):
  357. """Form the symbols s1, s2, ..., sn to stand for elem. symm. polys. """
  358. return symbols([f's{i + 1}' for i in range(n)])
  359. def sparse_symmetrize_resolvent_coeffs(F, X, s, verbose=False):
  360. """
  361. Compute the coefficients of a resolvent as functions of the coefficients of
  362. the associated polynomial.
  363. F must be a sparse polynomial.
  364. """
  365. import time, sys
  366. # Roots of resolvent as multivariate forms over vars X:
  367. root_forms = [
  368. F.compose(list(zip(X, sigma(X))))
  369. for sigma in s
  370. ]
  371. # Coeffs of resolvent (besides lead coeff of 1) as symmetric forms over vars X:
  372. Y = [Dummy(f'Y{i}') for i in range(len(s))]
  373. coeff_forms = []
  374. for i in range(1, len(s) + 1):
  375. if verbose:
  376. print('----')
  377. print(f'Computing symmetric poly of degree {i}...')
  378. sys.stdout.flush()
  379. t0 = time.time()
  380. G = symmetric_poly(i, *Y)
  381. t1 = time.time()
  382. if verbose:
  383. print(f'took {t1 - t0} seconds')
  384. print('lambdifying...')
  385. sys.stdout.flush()
  386. t0 = time.time()
  387. C = lambdify(Y, (-1)**i*G)
  388. t1 = time.time()
  389. if verbose:
  390. print(f'took {t1 - t0} seconds')
  391. sys.stdout.flush()
  392. coeff_forms.append(C)
  393. coeffs = []
  394. for i, f in enumerate(coeff_forms):
  395. if verbose:
  396. print('----')
  397. print(f'Plugging root forms into elem symm poly {i+1}...')
  398. sys.stdout.flush()
  399. t0 = time.time()
  400. g = f(*root_forms)
  401. t1 = time.time()
  402. coeffs.append(g)
  403. if verbose:
  404. print(f'took {t1 - t0} seconds')
  405. sys.stdout.flush()
  406. # Now symmetrize these coeffs. This means recasting them as polynomials in
  407. # the elementary symmetric polys over X.
  408. symmetrized = []
  409. symmetrization_times = []
  410. ss = s_vars(len(X))
  411. for i, A in list(enumerate(coeffs)):
  412. if verbose:
  413. print('-----')
  414. print(f'Coeff {i+1}...')
  415. sys.stdout.flush()
  416. t0 = time.time()
  417. B, rem, _ = A.symmetrize()
  418. t1 = time.time()
  419. if rem != 0:
  420. msg = f"Got nonzero remainder {rem} for resolvent (F, X, s) = ({F}, {X}, {s})"
  421. raise ResolventException(msg)
  422. B_str = str(B.as_expr(*ss))
  423. symmetrized.append(B_str)
  424. symmetrization_times.append(t1 - t0)
  425. if verbose:
  426. print(wrap(B_str))
  427. print(f'took {t1 - t0} seconds')
  428. sys.stdout.flush()
  429. return symmetrized, symmetrization_times
  430. def define_resolvents():
  431. """Define all the resolvents for polys T of degree 4 through 6. """
  432. from sympy.combinatorics.galois import PGL2F5
  433. from sympy.combinatorics.permutations import Permutation
  434. R4, X4 = xring("X0,X1,X2,X3", ZZ, lex)
  435. X = X4
  436. # The one resolvent used in `_galois_group_degree_4_lookup()`:
  437. F40 = X[0]*X[1]**2 + X[1]*X[2]**2 + X[2]*X[3]**2 + X[3]*X[0]**2
  438. s40 = [
  439. Permutation(3),
  440. Permutation(3)(0, 1),
  441. Permutation(3)(0, 2),
  442. Permutation(3)(0, 3),
  443. Permutation(3)(1, 2),
  444. Permutation(3)(2, 3),
  445. ]
  446. # First resolvent used in `_galois_group_degree_4_root_approx()`:
  447. F41 = X[0]*X[2] + X[1]*X[3]
  448. s41 = [
  449. Permutation(3),
  450. Permutation(3)(0, 1),
  451. Permutation(3)(0, 3)
  452. ]
  453. R5, X5 = xring("X0,X1,X2,X3,X4", ZZ, lex)
  454. X = X5
  455. # First resolvent used in `_galois_group_degree_5_hybrid()`,
  456. # and only one used in `_galois_group_degree_5_lookup_ext_factor()`:
  457. F51 = ( X[0]**2*(X[1]*X[4] + X[2]*X[3])
  458. + X[1]**2*(X[2]*X[0] + X[3]*X[4])
  459. + X[2]**2*(X[3]*X[1] + X[4]*X[0])
  460. + X[3]**2*(X[4]*X[2] + X[0]*X[1])
  461. + X[4]**2*(X[0]*X[3] + X[1]*X[2]))
  462. s51 = [
  463. Permutation(4),
  464. Permutation(4)(0, 1),
  465. Permutation(4)(0, 2),
  466. Permutation(4)(0, 3),
  467. Permutation(4)(0, 4),
  468. Permutation(4)(1, 4)
  469. ]
  470. R6, X6 = xring("X0,X1,X2,X3,X4,X5", ZZ, lex)
  471. X = X6
  472. # First resolvent used in `_galois_group_degree_6_lookup()`:
  473. H = PGL2F5()
  474. term0 = X[0]**2*X[5]**2*(X[1]*X[4] + X[2]*X[3])
  475. terms = {term0.compose(list(zip(X, s(X)))) for s in H.elements}
  476. F61 = sum(terms)
  477. s61 = [Permutation(5)] + [Permutation(5)(0, n) for n in range(1, 6)]
  478. # Second resolvent used in `_galois_group_degree_6_lookup()`:
  479. F62 = X[0]*X[1]*X[2] + X[3]*X[4]*X[5]
  480. s62 = [Permutation(5)] + [
  481. Permutation(5)(i, j + 3) for i in range(3) for j in range(3)
  482. ]
  483. return {
  484. (4, 0): (F40, X4, s40),
  485. (4, 1): (F41, X4, s41),
  486. (5, 1): (F51, X5, s51),
  487. (6, 1): (F61, X6, s61),
  488. (6, 2): (F62, X6, s62),
  489. }
  490. def generate_lambda_lookup(verbose=False, trial_run=False):
  491. """
  492. Generate the whole lookup table of coeff lambdas, for all resolvents.
  493. """
  494. jobs = define_resolvents()
  495. lambda_lists = {}
  496. total_time = 0
  497. time_for_61 = 0
  498. time_for_61_last = 0
  499. for k, (F, X, s) in jobs.items():
  500. symmetrized, times = sparse_symmetrize_resolvent_coeffs(F, X, s, verbose=verbose)
  501. total_time += sum(times)
  502. if k == (6, 1):
  503. time_for_61 = sum(times)
  504. time_for_61_last = times[-1]
  505. sv = s_vars(len(X))
  506. head = f'lambda {", ".join(str(v) for v in sv)}:'
  507. lambda_lists[k] = ',\n '.join([
  508. f'{head} ({wrap(f)})'
  509. for f in symmetrized
  510. ])
  511. if trial_run:
  512. break
  513. table = (
  514. "# This table was generated by a call to\n"
  515. "# `sympy.polys.numberfields.galois_resolvents.generate_lambda_lookup()`.\n"
  516. f"# The entire job took {total_time:.2f}s.\n"
  517. f"# Of this, Case (6, 1) took {time_for_61:.2f}s.\n"
  518. f"# The final polynomial of Case (6, 1) alone took {time_for_61_last:.2f}s.\n"
  519. "resolvent_coeff_lambdas = {\n")
  520. for k, L in lambda_lists.items():
  521. table += f" {k}: [\n"
  522. table += " " + L + '\n'
  523. table += " ],\n"
  524. table += "}\n"
  525. return table
  526. def get_resolvent_by_lookup(T, number):
  527. """
  528. Use the lookup table, to return a resolvent (as dup) for a given
  529. polynomial *T*.
  530. Parameters
  531. ==========
  532. T : Poly
  533. The polynomial whose resolvent is needed
  534. number : int
  535. For some degrees, there are multiple resolvents.
  536. Use this to indicate which one you want.
  537. Returns
  538. =======
  539. dup
  540. """
  541. degree = T.degree()
  542. L = resolvent_coeff_lambdas[(degree, number)]
  543. T_coeffs = T.rep.rep[1:]
  544. return [ZZ(1)] + [c(*T_coeffs) for c in L]
  545. # Use
  546. # (.venv) $ python -m sympy.polys.numberfields.galois_resolvents
  547. # to reproduce the table found in resolvent_lookup.py
  548. if __name__ == "__main__":
  549. import sys
  550. verbose = '-v' in sys.argv[1:]
  551. trial_run = '-t' in sys.argv[1:]
  552. table = generate_lambda_lookup(verbose=verbose, trial_run=trial_run)
  553. print(table)