Re: How to get an answer as a Root object?

*To*: mathgroup at smc.vnet.net*Subject*: [mg57156] Re: [mg57137] How to get an answer as a Root object?*From*: Daniel Lichtblau <danl at wolfram.com>*Date*: Thu, 19 May 2005 03:08:33 -0400 (EDT)*References*: <200505170520.BAA25934@smc.vnet.net>*Sender*: owner-wri-mathgroup at wolfram.com

David W. Cantrell wrote: > This question must have an obvious answer! > > A certain expression, given below my signature, involves nothing worse than > some nested square roots. It equals 0 for a value of the variable s near > 2.707, easily approximated using FindRoot. However, I would like to get > that value expressed precisely as a Root object. Trying to use Solve or > Reduce to do that, Mathematica seems to take "forever". > > Of course, I could do this by hand; in the days before CASs -- before I got > spoiled by them! -- I wouldn't have hesitated to do so. But surely there > must be an easy way to get Mathematica to give that value as a Root object. > How? > > David Cantrell > > > (-1 + Sqrt[2] - Sqrt[(1 + 2*Sqrt[2] - s)*(-1 + s)] + s)*(-1 - Sqrt[2] + > Sqrt[2]*s - Sqrt[-1 - 2*Sqrt[2] + 2*(2 + Sqrt[2] - s)*s]) - > (-1 - Sqrt[2] + Sqrt[2]*s + Sqrt[-1 - 2*Sqrt[2] + 2*(2 + Sqrt[2] - s)*s])* > (-2 - 2*Sqrt[2] + 2*s + 2*Sqrt[1 - (1/4)*(-1 - Sqrt[2] + Sqrt[(1 + > 2*Sqrt[2] - s)*(-1 + s)] + s)^2]) There are a few reasonable methods. Here are some. expr = (-1 + Sqrt[2] - Sqrt[(1 + 2*Sqrt[2] - s)*(-1 + s)] + s)*(-1 - Sqrt[2] + Sqrt[2]*s - Sqrt[-1 - 2*Sqrt[2] + 2*(2 + Sqrt[2] - s)*s]) - (-1 - Sqrt[2] + Sqrt[2]*s + Sqrt[-1 - 2*Sqrt[2] + 2*(2 + Sqrt[2] - s)*s])* (-2 - 2*Sqrt[2] + 2*s + 2*Sqrt[1 - (1/4)*(-1 - Sqrt[2] + Sqrt[(1 + 2*Sqrt[2] - s)*(-1 + s)] + s)^2]); Method 1. First we use numerical approximation to recognize a root polynomial. Due, I suspect, to nested radicals, we'll get a small spurious imaginary part which I remove using Re. In[5]:= InputForm[srootval = Re[s/.sroot]] Out[5]//InputForm= 2.706945837045705738269073202812169790793196895164013983001575249106155514352\ 266542032260316431664549796229819466076004645196713561849037`120. To find a candidate polynomial we'll use Recognize (uses a method based on lattice reduction), then factor, and find the factor with the right root. <<NumberTheory`Recognize` fax = Drop[Map[First,FactorList[Recognize[srootval, 30, s]]],1]; InputForm[minpoly = First[Select[fax, Chop[#/.sroot,10^(-acc)]===0&]]] Why is the factorization needed? Because the larger product polynomial has smaller vector norm in the lattice than does the actual minimal polynomial. Method 2. We'll compute a Groebner basis with 's' the main variable. The hope, born out in this case, is that internal handling of algebraics will force to the fore a "polynomial" in s. gb = GroebnerBasis[expr, s]; I do not show it here, but the polynomial we want is readily seen to be the first one. Again we'll factor and find the one with the correct root. As there are algebraics in the coefficients of this polynomial we use Extension->Automatic in the factoring. poly = First[gb]; fax2 = Drop[Map[First,FactorList[poly, Extension->Automatic]]]; In[21]:= InputForm[minpoly2 = First[Select[fax2, Chop[#/.sroot,10^(-acc)]===0&]]] Out[21]//InputForm= -136 - 96*Sqrt[2] + (512 + 364*Sqrt[2])*s + (-804 - 578*Sqrt[2])*s^2 + (676 + 486*Sqrt[2])*s^3 + (-329 - 220*Sqrt[2])*s^4 + (96 + 44*Sqrt[2])*s^5 - 16*s^6 How to represent as a polynomial over the integers? Just multiply by an algebraic conjugate. sqrtConjugate[ee_] := ee /. aa_^(1/2)->-aa^(1/2) In[23]:= InputForm[minpoly2 = Expand[minpoly2*sqrtConjugate[minpoly2]]] Out[23]//InputForm= 64 + 512*s - 6112*s^2 + 21024*s^3 - 32136*s^4 + 10832*s^5 + 43568*s^6 - 86152*s^7 + 81425*s^8 - 46080*s^9 + 15872*s^10 - 3072*s^11 + 256*s^12 Method 3. Replace radicals by hand, use defining relations explicitly, formulate Groebner basis, find polynomial in s alone. rep1 = {Sqrt[2]->rad0}; rep2 = {Sqrt[(1+2*Sqrt[2]-s)*(-1+s)] -> rad1} /. rep1; rep3 = ({Sqrt[-1-2*Sqrt[2]+2*(2+Sqrt[2]-s)*s] -> rad2, Sqrt[1-(1/4)*(-1-Sqrt[2]+Sqrt[(1+2*Sqrt[2]-s)*(-1+s)]+s)^2] -> rad3} /. rep1 /. rep2) repssq = Map[#[[1]]^2-#[[2]]^2&, Flatten[{rep1,rep2,rep3}]] expr2 = (expr /. rep1) /. reps; In[64]:= InputForm[Timing[poly4 = First[GroebnerBasis[Join[repssq,{expr2}], s, {rad0,rad1,rad2,rad3}, MonomialOrder->EliminationOrder]]]] Out[64]//InputForm= {0.30000000000000115*Second, 200704 + 3956736*s - 16568320*s^2 - 259776512*s^3 + 2455161856*s^4 - 9217882112*s^5 + 16736476672*s^6 - 4378572288*s^7 - 50113953728*s^8 + 119132812800*s^9 - 109419876224*s^10 - 38011997376*s^11 + 225422464224*s^12 - 259623250704*s^13 + 80225006088*s^14 + 142924054128*s^15 - 196411664327*s^16 + 53713807776*s^17 + 132853129136*s^18 - 220218741248*s^19 + 191796620928*s^20 - 115209149440*s^21 + 51080841728*s^22 - 17021165568*s^23 + 4236362752*s^24 - 766771200*s^25 + 95551488*s^26 - 7340032*s^27 + 262144*s^28} Again one would factor and find the one with the correct root. I'm sure there is a resultant formulation as well that will get a polynomial in s for the root in question. Daniel Lichtblau Wolfram Research

**Follow-Ups**:**Re: Re: How to get an answer as a Root object?***From:*Chris Chiasson <chris.chiasson@gmail.com>

**Re: Re: How to get an answer as a Root object?***From:*DrBob <drbob@bigfoot.com>

**References**:**How to get an answer as a Root object?***From:*"David W. Cantrell" <DWCantrell@sigmaxi.org>