Arithmetic on algebraic numbers

Question

I'd like to perform some elementary operations on algebraic numbers.

p1 = x^5 + 6*x^4 - 42*x^3 - 142*x^2 + 467*x + 422;
p2 = Expand[p1 /. {x -> ((x - 1)^2)}];
r1 = Root[{p2 /. x -> # &, 1 + 2.14850301089246970680618135585609253400250*I}];
r2 = Root[{p2 /. x -> # &, 1 + 2.87513448574410025133873802142153923543448*I}];
Timing[MinimalPolynomial[r1-r2]]

I'm trying to compare Mathematica's behavior to what I'm currently working on with Sage. In particular I'd like to see whether division r1/r2 takes much longer than r1+r2, r1-r2, r1*r2.

But it turns out that Mathematica apparently isn't able to perform elementary arithmetic operations on these algebraic numbers at all. The above example fails with

MinimalPolynomial::nalg: "-Root[…]+Root[…] is not an explicit algebraic number."

Is there some formulation I can use to perform such computations? Or is Mathematica simply not up to this task?

Answer 1

Root[{f, approx}] format was meant for representing roots of non-polynomial functions. Mathematica will never produce algebraic numbers (of degree <= $MaxRootDegree) represented in this way. If you manually enter Root[{f, approx}] with a polynomial f and a correct approximation the result is a valid exact number, but it is not recognized as an algebraic number. I will make RootReduce check for such objects and convert them to the standard algebraic number representation -- Root[poly, root number]. For now you can automate the method suggested by Daniel.

In[1]:= rootFromApprox[poly_, x_, approx_]:=
   Select[RootReduce[x/.Solve[poly==0, x, Cubics->False, Quartics->False]],
      #-approx==0&][[1]]

In[2]:= p1=x^5+6*x^4-42*x^3-142*x^2+467*x+422;
   p2=Expand[p1/.{x->((x-1)^2)}];
   r1=rootFromApprox[p2, x, 1+2.14850301089246970680618135585609253400250*I];
   r2=rootFromApprox[p2, x, 1+2.87513448574410025133873802142153923543448*I];
   Timing[Short[MinimalPolynomial[#]]]&/@{r1+r2, r1-r2, r1 r2, r1/r2}

Out[6]= {{0.124801,-3568379502241854768+<<58>>+#1^40&},
   {0.140401,<<34>>+48 #1^38+#1^40&},  
   {0.140401,66045000696445844586496+82370731205679648866304 #1+<<58>>+#1^40&},
   {0.327602,<<118>>+257988283970491580416 #1^80&}}

Answer 2

I do not know what is going awry and have sent a question to the relevant person here. I can suggest an alternative route. Form the Root objects using say Solve, figure out which are the ones you want, and use those.

p1 = x^5 + 6*x^4 - 42*x^3 - 142*x^2 + 467*x + 422;
p2 = Expand[p1 /. {x -> ((x - 1)^2)}];
rts = RootReduce[x /. Solve[p2 == 0, x]];
N[rts]

(* Out[38]= {1. - 2.87513448574 I, 1. + 2.87513448574 I, 
 1. - 2.14850301089 I, 1. + 2.14850301089 I, 1. - 0.876768035257 I, 
 1. + 0.876768035257 I, -0.823985071296, 2.8239850713, \
-1.07948651024, 3.07948651024} *)

So it is the fourth and second that we want.

{rt1, rt2} = rts[[{4, 2}]]
Timing[MinimalPolynomial[rt1 - rt2]]

(* Out[40]= {Root[
  712 - 172 #1 - 802 #1^2 + 952 #1^3 - 142 #1^4 - 336 #1^5 + 
    336 #1^6 - 168 #1^7 + 51 #1^8 - 10 #1^9 + #1^10 &, 8], 
 Root[712 - 172 #1 - 802 #1^2 + 952 #1^3 - 142 #1^4 - 336 #1^5 + 
    336 #1^6 - 168 #1^7 + 51 #1^8 - 10 #1^9 + #1^10 &, 6]}

Out[41]= {0.050724, 
 515591535280592 - 6604303796727744 #1^2 - 19213629465035648 #1^4 - 
   9238448577096832 #1^6 - 246985961012832 #1^8 - 
   408212448315616 #1^10 - 179542482936184 #1^12 - 
   20639353188832 #1^14 - 6030161483696 #1^16 - 874984589616 #1^18 - 
   187565974272 #1^20 - 33966863800 #1^22 - 4586283639 #1^24 - 
   530873932 #1^26 - 40291328 #1^28 - 2037008 #1^30 + 5146 #1^32 + 
   10684 #1^34 + 924 #1^36 + 48 #1^38 + #1^40 &} *)