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

*To*: mathgroup at smc.vnet.net*Subject*: [mg57189] Re: [mg57156] Re: [mg57137] How to get an answer as a Root object?*From*: DrBob <drbob at bigfoot.com>*Date*: Fri, 20 May 2005 04:43:23 -0400 (EDT)*References*: <200505170520.BAA25934@smc.vnet.net> <200505190708.DAA13114@smc.vnet.net>*Reply-to*: drbob at bigfoot.com*Sender*: owner-wri-mathgroup at wolfram.com

Two questions, Daniel: 1) Is this the missing definition of sroot? sroot = FindRoot[expr == 0, {s, 2.7}] 2) Doesn't Drop do nothing at all in this statement? fax2 = Drop[Map[First,FactorList[poly, Extension->Automatic]]]; Bobby On Thu, 19 May 2005 03:08:33 -0400 (EDT), Daniel Lichtblau <danl at wolfram.com> wrote: > 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 > > > > > -- DrBob at bigfoot.com

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

**Re: How to get an answer as a Root object?***From:*Daniel Lichtblau <danl@wolfram.com>