# Efficient algorithm for calculating polynomial that has roots of a certain form

I have looked into this specific question on Math.SE concerning a more "mechanical" approach to finding a polynomial $$p \in \mathbb{Q}[x]$$ satisfying $$p(\sqrt{2}+\sqrt{3}) = 0$$. The user MJD sketches an approach that could conceivably be put into algorithmic terms:

I will sketch his approach more generally by assuming that we seek to find $$p$$ as above with $$p(\sqrt{a}+\sqrt{b}) = 0$$ for $$a,b \in \mathbb{N}$$. He argues: Because all powers of $$r_{a,b}:= \sqrt{a}+\sqrt{b}$$ can be expressed as a rational linear combination of $$1$$, $$\sqrt{a}$$, $$\sqrt{b}$$ and $$\sqrt{ab}$$ (so a vector space of dimension $$4$$ over $$\mathbb{Q}$$), we can write $$r_{a,b}^{4} = a_3 r_{a,b}^3 + a_2 r_{a,b}^2 + a_1 r_{a,b}+a_0$$

for adequate $$a_i \in \mathbb{Q}$$. By construction, this yields a polynomial of the wanted form.

The $$a_i$$ are found by comparing the coefficient of the above mentioned base vectors $$1$$, $$\sqrt{a}$$, $$\sqrt{b}$$ and $$\sqrt{ab}$$ for the calculated powers. In this case, we get

$$\begin{pmatrix} 1 & 0 & a+b & 0 \\ 0 & 1 & 0 & a+3b \\ 0 & 1 & 0 & 3a+b \\ 0 & 0 & 2 & 0 \end{pmatrix} \begin{pmatrix}a_0 \\ a_1 \\ a_2 \\ a_3\end{pmatrix} = \begin{pmatrix}a^2 + b^2 + 6ab \\ 0 \\ 0 \\ 4(a+b)\end{pmatrix}$$ (I am sorry for the extremly bad formatting, but the low reputation filter prevents me from writing perfectly normal LaTeX code)

By close examination or by inverting the displayed matrix (not assuming specific $$a$$ and $$b$$), we can identify $$p$$ as $$p(x) := x^4-2(a+b)x^2+(a-b)^2$$ which after inspection turns out to actually have the root $$r_{a,b}$$.

The same can be done for the $$\sqrt{a}+\sqrt[3]{b}$$ case, giving $$p(x) := x^6 - 3ax^4 - 2bx^3+3a^2x^2-6abx+b^2-a^3$$

All of this strongly suggests that with a bit of clever manipulation, one can automate this process. Maybe not efficiently, but I am willing to wait the time out.

Everything comes down to the following problem: It is quite laborious to single out the basis elements for any specific case and the powers of certain expressions can get complicated quite quickly, for example $$\sqrt{\sqrt[3]{a}+\sqrt{b}}+\sqrt[4]{c}$$. This is where I involved Mathematica.

My approach was:

• Calculate the powers of the term in question.

• Group terms with the same basis vector together (Mathematica is hilariously bad at factoring out whole powers from fractional exponents)

• Somehow sort the terms of the expression so that they can be subsequently treated to a matrix representation (filling up with zeroes for missing coefficients if necessary)

• Get the coefficients by left-multiplying the vector of the highest order necessary with the inverse of the just generated matrix

I got a prototype working for the case $$\sqrt[p]{a}+\sqrt[q]{b}$$, but it is very slow and it already fails at nested roots. I'm not really sure if my solution is even consistent, as it relies on pattern recognition of powers and follows an iterative approach. I can post my code if requested, but it is ugly beyond comparison and would just clutter up this already way too large question.

So my question is:

Given an parametric algebraic number $$s$$ that can be expressed as the sum of finitely many $$n$$-th (nested) roots of parameters and assuming we can figure out the dimension of $$\mathbb{Q}(s)$$ over $$\mathbb{Q}$$, can an algorithmic approach to finding a rational polynomial $$p$$ with $$p(s) = 0$$ be formulated in Mathematica?

I'm happy if even a bodged solution should do the job, as I am just interested in finding out more about when these generated polynomials actually are minimal polynomials in the algebraic sense.

If this question misses the scope of this forum, feel free to migrate it to wherever you deem it fit.

Thank you and have a good day!

## 2 Answers

In Mathematica this can be done using GroebnerBasis or the resource function DixonResultant. Remarkably similar to something I just responded to a few minutes ago.

First[
GroebnerBasis[{y - (sqrt2 + sqrt3), sqrt2^2 - 2, sqrt3^2 - 3},
y, {sqrt2, sqrt3}]]

(* Out[611]= 1 - 10 y^2 + y^4 *)

ResourceFunction["DixonResultant"][{y - (sqrta + cbrtb), sqrta^2 - a,
cbrtb^3 - b}, {sqrta, cbrtb}]

(* Out[614]= -a^3 + b^2 - 6 a b y + 3 a^2 y^2 - 2 b y^3 - 3 a y^4 + y^6 *)


Be warned that this can become intractable fairly easily. The nested example takes considerable time, even with some undocumented option settings.

Timing[
First[GroebnerBasis[polys, expr, {sqrtbig, rt4c, cbrta, sqrtb},
MonomialOrder -> EliminationOrder,
CoefficientDomain -> RationalFunctions,
Method -> {"GroebnerWalk", "EarlyEliminate" -> True}]]]

(* Out[621]= {4628.73, -a^8 + 4 a^6 b^3 - 6 a^4 b^6 + 4 a^2 b^9 - b^12 +
60 a^6 b^2 c - 108 a^4 b^5 c + 36 a^2 b^8 c + 12 b^11 c +
60 a^6 b c^2 - 1050 a^4 b^4 c^2 - 240 a^2 b^7 c^2 - 66 b^10 c^2 +
4 a^6 c^3 - 1768 a^4 b^3 c^3 + 464 a^2 b^6 c^3 + 220 b^9 c^3 -
1050 a^4 b^2 c^4 - 264 a^2 b^5 c^4 - 495 b^8 c^4 - 108 a^4 b c^5 -
264 a^2 b^4 c^5 + 792 b^7 c^5 - 6 a^4 c^6 + 464 a^2 b^3 c^6 -
924 b^6 c^6 - 240 a^2 b^2 c^7 + 792 b^5 c^7 + 36 a^2 b c^8 -
495 b^4 c^8 + 4 a^2 c^9 + 220 b^3 c^9 - 66 b^2 c^10 + 12 b c^11 -
c^12 + (24 a^7 b - 72 a^5 b^4 + 72 a^3 b^7 - 24 a b^10 +
120 a^7 c + 1440 a^5 b^3 c - 648 a^3 b^6 c - 912 a b^9 c +
6480 a^5 b^2 c^2 - 9144 a^3 b^5 c^2 + 5256 a b^8 c^2 +
1440 a^5 b c^3 - 38472 a^3 b^4 c^3 - 9408 a b^7 c^3 - 72 a^5 c^4
+ 16152 a^3 b^3 c^4 + 336 a b^6 c^4 + 23976 a^3 b^2 c^5 +
22176 a b^5 c^5 + 8280 a^3 b c^6 - 34608 a b^4 c^6 -
216 a^3 c^7 + 24384 a b^3 c^7 - 7992 a b^2 c^8 + 624 a b c^9 +
168 a c^10) expr^2 + (-228 a^6 b^2 + 468 a^4 b^5 - 252 a^2 b^8 +
12 b^11 + 1320 a^6 b c - 5820 a^4 b^4 c + 6816 a^2 b^7 c +
276 b^10 c - 5220 a^6 c^2 + 104040 a^4 b^3 c^2 -
98064 a^2 b^6 c^2 - 2268 b^9 c^2 + 12648 a^4 b^2 c^3 -
16800 a^2 b^5 c^3 + 6780 b^8 c^3 - 110652 a^4 b c^4 -
71400 a^2 b^4 c^4 - 9672 b^7 c^4 - 6828 a^4 c^5 +
519264 a^2 b^3 c^5 + 4872 b^6 c^5 - 314640 a^2 b^2 c^6 +
4872 b^5 c^6 - 16224 a^2 b c^7 - 9672 b^4 c^7 - 8700 a^2 c^8 +
6780 b^3 c^8 - 2268 b^2 c^9 + 276 b c^10 +
12 c^11) expr^4 + (8 a^7 + 1056 a^5 b^3 - 1272 a^3 b^6 +
208 a b^9 - 3360 a^5 b^2 c + 11472 a^3 b^5 c - 13296 a b^8 c -
90912 a^5 b c^2 + 263736 a^3 b^4 c^2 - 79296 a b^7 c^2 +
99872 a^5 c^3 + 427616 a^3 b^3 c^3 - 258496 a b^6 c^3 -
1670664 a^3 b^2 c^4 + 899424 a b^5 c^4 + 855504 a^3 b c^5 -
282144 a b^4 c^5 + 129992 a^3 c^6 - 504000 a b^3 c^6 +
118080 a b^2 c^7 + 117840 a b c^8 +
1680 a c^9) expr^6 + (-132 a^6 b - 2394 a^4 b^4 + 1296 a^2 b^7 -
66 b^10 + 1020 a^6 c + 25320 a^4 b^3 c - 17808 a^2 b^6 c -
2268 b^9 c - 733212 a^4 b^2 c^2 + 776016 a^2 b^5 c^2 -
40986 b^8 c^2 + 1429224 a^4 b c^3 - 1055952 a^2 b^4 c^3 +
97968 b^7 c^3 - 821082 a^4 c^4 - 494544 a^2 b^3 c^4 +
14940 b^6 c^4 + 4261968 a^2 b^2 c^5 - 139176 b^5 c^5 -
1290384 a^2 b c^6 + 14940 b^4 c^6 - 214512 a^2 c^7 +
97968 b^3 c^7 - 40986 b^2 c^8 - 2268 b c^9 -
66 c^10) expr^8 + (720 a^5 b^2 + 2664 a^3 b^5 - 792 a b^8 +
30816 a^5 b c - 104088 a^3 b^4 c + 25536 a b^7 c -
79152 a^5 c^2 - 435504 a^3 b^3 c^2 - 75168 a b^6 c^2 +
1244688 a^3 b^2 c^3 + 85056 a b^5 c^3 - 7544376 a^3 b c^4 +
7046256 a b^4 c^4 + 2044296 a^3 c^5 + 2338368 a b^3 c^5 +
2374752 a b^2 c^6 + 798144 a b c^7 -
9240 a c^8) expr^10 + (-28 a^6 - 1128 a^4 b^3 - 2736 a^2 b^6 +
220 b^9 - 147864 a^4 b^2 c + 136800 a^2 b^5 c + 6780 b^8 c +
720744 a^4 b c^2 - 1188432 a^2 b^4 c^2 + 97968 b^7 c^2 +
942232 a^4 c^3 - 2924736 a^2 b^3 c^3 + 2546800 b^6 c^3 -
16836432 a^2 b^2 c^4 + 5736840 b^5 c^4 - 2889120 a^2 b c^5 +
5736840 b^4 c^5 - 2206640 a^2 c^6 + 2546800 b^3 c^6 +
97968 b^2 c^7 + 6780 b c^8 + 220 c^9) expr^12 + (288 a^5 b -
1368 a^3 b^4 + 1728 a b^7 - 4320 a^5 c + 231648 a^3 b^3 c +
79680 a b^6 c + 866160 a^3 b^2 c^2 - 1394496 a b^5 c^2 +
6007776 a^3 b c^3 - 6522048 a b^4 c^3 + 1518888 a^3 c^4 -
8438208 a b^3 c^4 - 254016 a b^2 c^5 + 787008 a b c^6 +
11712 a c^7) expr^14 + (-666 a^4 b^2 + 2808 a^2 b^5 - 495 b^8 +
19332 a^4 b c + 80856 a^2 b^4 c - 9672 b^7 c - 647322 a^4 c^2 +
1337904 a^2 b^3 c^2 + 14940 b^6 c^2 + 5581872 a^2 b^2 c^3 +
5736840 b^5 c^3 - 11392296 a^2 b c^4 + 10536870 b^4 c^4 -
5896200 a^2 c^5 + 5736840 b^3 c^5 + 14940 b^2 c^6 - 9672 b c^7 -
495 c^8) expr^16 + (56 a^5 - 296 a^3 b^3 - 2352 a b^6 +
174792 a^3 b^2 c - 253152 a b^5 c + 1673832 a^3 b c^2 -
1313616 a b^4 c^2 - 2480648 a^3 c^3 + 4522432 a b^3 c^3 -
6233808 a b^2 c^4 - 2478816 a b c^5 +
8528 a c^6) expr^18 + (-300 a^4 b - 1224 a^2 b^4 + 792 b^7 +
4692 a^4 c - 494496 a^2 b^3 c + 4872 b^6 c -
4144176 a^2 b^2 c^2 - 139176 b^5 c^2 - 7300512 a^2 b c^3 +
5736840 b^4 c^3 - 6115272 a^2 c^4 + 5736840 b^3 c^4 -
139176 b^2 c^5 + 4872 b c^6 + 792 c^7) expr^20 + (24 a^3 b^2 +
2016 a b^5 - 116400 a^3 b c + 222432 a b^4 c - 687912 a^3 c^2 +
1556928 a b^3 c^2 + 3132096 a b^2 c^3 - 775584 a b c^4 -
33696 a c^5) expr^22 + (-70 a^4 - 48 a^2 b^3 - 924 b^6 +
244848 a^2 b^2 c + 4872 b^5 c + 823152 a^2 b c^2 +
14940 b^4 c^2 - 2074928 a^2 c^3 + 2546800 b^3 c^3 +
14940 b^2 c^4 + 4872 b c^5 - 924 c^6) expr^24 + (120 a^3 b -
1008 a b^4 - 168 a^3 c - 24384 a b^3 c + 927072 a b^2 c^2 +
1156800 a b c^3 + 30480 a c^4) expr^26 + (144 a^2 b^2 + 792 b^5 +
45216 a^2 b c - 9672 b^4 c - 258672 a^2 c^2 + 97968 b^3 c^2 +
97968 b^2 c^3 - 9672 b c^4 + 792 c^5) expr^28 + (56 a^3 +
192 a b^3 - 56256 a b^2 c + 374592 a b c^2 -
9280 a c^3) expr^30 + (36 a^2 b - 495 b^4 - 2268 a^2 c +
6780 b^3 c - 40986 b^2 c^2 + 6780 b c^3 -
495 c^4) expr^32 + (72 a b^2 + 19440 a b c -
1272 a c^2) expr^34 + (-28 a^2 + 220 b^3 - 2268 b^2 c -
2268 b c^2 + 220 c^3) expr^36 + (-48 a b +
912 a c) expr^38 + (-66 b^2 + 276 b c - 66 c^2) expr^40 +
8 a expr^42 + (12 b + 12 c) expr^44 - expr^48} *)

• I may have forgot to thank you for your solution at the time, sorry! I may not fully understand the inner workings of how to calculate a Groebner basis or Dixon resultant yet, but your solution works just as advertised (and the problem is probably NP-hard anyway, so no hard feelings about speed). Thank you for your insights! Apr 24, 2022 at 11:33
• Glad it helped. The "naive" method by @MichaelE2 might well be the fast approach though. Apr 24, 2022 at 15:43

Somewhat naive approach to powers: Form conjugates and multiply out. Intermediate simplification helps with speed and memory use. Map distinct powers to inert forms to prevent autosimplification.

nn = 0;
p // ClearAll;
p[n_][x_List, pow_] := p[n][#, pow] & /@ x;
expr = Sqrt[a^(1/3) + Sqrt@b] + c^(1/4) /.
Power[b_, r_Rational] :> p[++nn][b, r]
pwrs = Cases[
expr,
p[n_][a_, b_], Infinity];
pwrs[[All, 1]] = Blank[];
pwrs
(* expr and pwrs:
p[1][p[2][a, 1/3] + p[3][b, 1/2], 1/2] + p[4][c, 1/4]
{p[2][_, 1/3], p[3][_, 1/2], p[1][_, 1/2], p[4][_, 1/4]}
*)

(poly = Fold[
Collect[
Times @@ (# /.
p : #2 :> Power @@ p*Exp[2 Pi I Range[1/Last@p] Last@p]),
x,
Simplify] &,
x - expr,
pwrs
] // Collect[#, x, Simplify] &) // AbsoluteTiming
(* = the negative of @Daniel's
{1.32343, (a^4 + (b - c)^6 -
2 a^2 (b^3 + 15 b^2 c + 15 b c^2 + c^3))^2 -
24 (a^7 (b + 5 c) -
a (b - c)^7 (b^3 + 45 b^2 c + 75 b c^2 + 7 c^3) -
3 a^5 (b^4 - 20 b^3 c - 90 b^2 c^2 - 20 b c^3 + c^4) +
a^3 (3 b^7 - 27 b^6 c - 381 b^5 c^2 - 1603 b^4 c^3 +
673 b^3 c^4 + 999 b^2 c^5 + 345 b c^6 - 9 c^7)) x^2 +
12 (a^6 (19 b^2 - 110 b c + 435 c^2) - (b - c)^8 (b^3 + 31 b^2 c +
31 b c^2 + c^3) +
a^4 (-39 b^5 + 485 b^4 c - 8670 b^3 c^2 - 1054 b^2 c^3 +
9221 b c^4 + 569 c^5) +
a^2 (b - c)^2 (21 b^6 - 526 b^5 c + 7099 b^4 c^2 +
16124 b^3 c^3 + 31099 b^2 c^4 + 2802 b c^5 + 725 c^6)) x^4 -
8 (a^7 + 4 a^5 (33 b^3 - 105 b^2 c - 2841 b c^2 + 3121 c^3) +
2 a (b - c)^3 (13 b^6 - 792 b^5 c - 7371 b^4 c^2 -
35880 b^3 c^3 - 30105 b^2 c^4 - 7680 b c^5 - 105 c^6) +
a^3 (-159 b^6 + 1434 b^5 c + 32967 b^4 c^2 + 53452 b^3 c^3 -
208833 b^2 c^4 + 106938 b c^5 + 16249 c^6)) x^6 +
6 (2 a^6 (11 b - 85 c) +
a^4 (399 b^4 - 4220 b^3 c + 122202 b^2 c^2 - 238204 b c^3 +
136847 c^4) + (b - c)^4 (11 b^6 + 422 b^5 c + 8453 b^4 c^2 +
14996 b^3 c^3 + 8453 b^2 c^4 + 422 b c^5 + 11 c^6) -
8 a^2 (27 b^7 - 371 b^6 c + 16167 b^5 c^2 - 21999 b^4 c^3 -
10303 b^3 c^4 + 88791 b^2 c^5 - 26883 b c^6 -
4469 c^7)) x^8 -
24 a (-33 b^8 + 1064 b^7 c - 3132 b^6 c^2 + 3544 b^5 c^3 +
293594 b^4 c^4 + 97432 b^3 c^5 + 98948 b^2 c^6 + 33256 b c^7 -
385 c^8 + 2 a^4 (15 b^2 + 642 b c - 1649 c^2) +
a^2 (111 b^5 - 4337 b^4 c - 18146 b^3 c^2 + 51862 b^2 c^3 -
314349 b c^4 + 85179 c^5)) x^10 +
4 (7 a^6 - 55 b^9 - 1695 b^8 c - 24492 b^7 c^2 - 636700 b^6 c^3 -
1434210 b^5 c^4 - 1434210 b^4 c^5 - 636700 b^3 c^6 -
24492 b^2 c^7 - 1695 b c^8 - 55 c^9 +
2 a^4 (141 b^3 + 18483 b^2 c - 90093 b c^2 - 117779 c^3) +
4 a^2 (171 b^6 - 8550 b^5 c + 74277 b^4 c^2 + 182796 b^3 c^3 +
1052277 b^2 c^4 + 180570 b c^5 + 137915 c^6)) x^12 -
24 (12 a^5 (b - 15 c) +
a^3 (-57 b^4 + 9652 b^3 c + 36090 b^2 c^2 + 250324 b c^3 +
63287 c^4) +
8 a (9 b^7 + 415 b^6 c - 7263 b^5 c^2 - 33969 b^4 c^3 -
43949 b^3 c^4 - 1323 b^2 c^5 + 4099 b c^6 + 61 c^7)) x^14 +
3 (165 b^8 + 3224 b^7 c - 4980 b^6 c^2 - 1912280 b^5 c^3 -
3512290 b^4 c^4 - 1912280 b^3 c^5 - 4980 b^2 c^6 + 3224 b c^7 +
165 c^8 + a^4 (222 b^2 - 6444 b c + 215774 c^2) -
8 a^2 (117 b^5 + 3369 b^4 c + 55746 b^3 c^2 + 232578 b^2 c^3 -
474679 b c^4 - 245675 c^5)) x^16 -
8 (7 a^5 +
a^3 (-37 b^3 + 21849 b^2 c + 209229 b c^2 - 310081 c^3) -
2 a (147 b^6 + 15822 b^5 c + 82101 b^4 c^2 - 282652 b^3 c^3 +
389613 b^2 c^4 + 154926 b c^5 - 533 c^6)) x^18 +
12 (a^4 (25 b - 391 c) +
2 a^2 (51 b^4 + 20604 b^3 c + 172674 b^2 c^2 + 304188 b c^3 +
254803 c^4) -
2 (33 b^7 + 203 b^6 c - 5799 b^5 c^2 + 239035 b^4 c^3 +
239035 b^3 c^4 - 5799 b^2 c^5 + 203 b c^6 + 33 c^7)) x^20 -
24 (a^3 (b^2 - 4850 b c - 28663 c^2) +
4 a (21 b^5 + 2317 b^4 c + 16218 b^3 c^2 + 32626 b^2 c^3 -
8079 b c^4 - 351 c^5)) x^22 +
2 (35 a^4 + 462 b^6 - 2436 b^5 c - 7470 b^4 c^2 - 1273400 b^3 c^3 -
7470 b^2 c^4 - 2436 b c^5 + 462 c^6 +
8 a^2 (3 b^3 - 15303 b^2 c - 51447 b c^2 + 129683 c^3)) x^24 -
24 a (-42 b^4 + a^2 (5 b - 7 c) - 1016 b^3 c + 38628 b^2 c^2 +
48200 b c^3 + 1270 c^4) x^26 -
24 (33 b^5 - 403 b^4 c + 4082 b^3 c^2 + 4082 b^2 c^3 - 403 b c^4 +
33 c^5 + 2 a^2 (3 b^2 + 942 b c - 5389 c^2)) x^28 -
8 (7 a^3 +
8 a (3 b^3 - 879 b^2 c + 5853 b c^2 -
145 c^3)) x^30 + (495 b^4 - 36 a^2 (b - 63 c) - 6780 b^3 c +
40986 b^2 c^2 - 6780 b c^3 + 495 c^4) x^32 -
24 a (3 b^2 + 810 b c - 53 c^2) x^34 +
4 (7 a^2 - 55 b^3 + 567 b^2 c + 567 b c^2 - 55 c^3) x^36 +
48 a (b - 19 c) x^38 + (66 b^2 - 276 b c + 66 c^2) x^40 -
8 a x^42 - 12 (b + c) x^44 + x^48}
*)


A sophisticated way of getting expr and pwrs:

nn = 0;
p // ClearAll;
p[n_][x_List, pow_] := p[n][#, pow] & /@ x;
expr = FixedPoint[
Fold[ReplaceAll,
#,
Cases[#,
pwr : Inactive[Power][b_, r_Rational] :> (pwr -> p[++nn][b, r]),
Infinity]
] &,
Sqrt[a^(1/3) + Sqrt@b] + c^(1/4) //.
Power[b_, r_Rational] :>
Inactive[Power][b, 1/Denominator@r]^Numerator@r,
5]
pwrs = Reverse@Cases[
expr,
p[n_][a_, b_], Infinity];
pwrs[[All, 1]] = Blank[];


Note the order matters (Reverse@Cases[...]): Have to simplify from the outside in. Not sure how this will hold up with really complicated examples.

• This was the way I was originally going to try, but failed (probably) mostly because of my lack of knowledge of how to tackle problems with Mathematica. Interesting how it can still be done this way despite the very dynamic nature of the problem, thank you! Apr 24, 2022 at 11:36