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

*To*: mathgroup at smc.vnet.net*Subject*: [mg57198] Re: [mg57156] Re: [mg57137] How to get an answer as a Root object?*From*: Daniel Lichtblau <danl at wolfram.com>*Date*: Fri, 20 May 2005 04:43:51 -0400 (EDT)*References*: <200505170520.BAA25934@smc.vnet.net> <200505190708.DAA13114@smc.vnet.net> <acbec1a405051916553ad63c12@mail.gmail.com>*Sender*: owner-wri-mathgroup at wolfram.com

Chris Chiasson wrote: > Daniel, > > Thank you for taking the time to illustrate how to solve for > complicated roots using those exact methods. > > In Method 2, how do you know which part of gb is the correct one? > > In Method 1, there were some undefined variables, namely act and > sroot. I fixed the example code below: > > srootval=s/.(sroot= > MapAt[Re, > FindRoot[expr,{s,2.7},WorkingPrecision->240, > PrecisionGoal->120],{1,2}]); > <<NumberTheory`Recognize` > fax=Drop[Map[First,FactorList[Recognize[srootval,30,s]]],1] > minpoly=First[Select[fax,Chop[#/.sroot]===0&]] > sexactroot=Solve[minpoly==0][[2]] > sexactroot//N > > Thank you for your time, See bottom. For better or worse I decided to keep the whole post so it would at least be self contained. > On 5/19/05, 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 I inadvertantly left out the actual FindRoot. It was something along the lines you have. Here are the first few lines I managed to omit from my last post. acc = 100; sroot = FindRoot[expr==0, {s,2.7}, AccuracyGoal->acc, PrecisionGoal->acc, WorkingPrecision->(1.2*acc)]; As for figuring out which polynomial to use in the Groebner basis from method 2, I knew because I looked. I should have stated this outright so it did not resemble a magic trick. In brief, the expectation should be that numeric algebraics e.g Sqrt[2] will be ordered lowest, so the stated variable 's' will be interms of them. But algebraics containing s will be ordered higher so that we do not get a circularity expressing s in terms of an algebraic in s (which, in essence, is what the given equation is doing). In any case I would neither expect most users to realize this, not would I necessarily trust it to work that way myself without checking. Which is what I did. To finish off by extracting a Root object, one can do as below. As ever, we will select based on numerical value being equal to that we already found, srootval. In[21]:= InputForm[root = First[ Select[s/.Solve[minpoly2==0,s], #==srootval&]]] Out[21]//InputForm= 1 + Sqrt[1 - Root[64 + 256*#1 - 1464*#1^2 + 1896*#1^3 - 495*#1^4 - 512*#1^5 + 256*#1^6 & , 1, 0]] If we prefer a single root form, we can do In[24]:= InputForm[root = RootReduce[root]] Out[24]//InputForm= Root[64 + 512*#1 - 6112*#1^2 + 21024*#1^3 - 32136*#1^4 + 10832*#1^5 + 43568*#1^6 - 86152*#1^7 + 81425*#1^8 - 46080*#1^9 + 15872*#1^10 - 3072*#1^11 + 256*#1^12 & , 8, 0] My apologies for the omissions (which I hope this note will fix but I'm not going to guarantee it). Daniel Lichtblau Wolfram Research

**Follow-Ups**:**Re: Re: Re: How to get an answer as a Root object?***From:*Andrzej Kozlowski <akoz@mimuw.edu.pl>

**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>