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

*To*: mathgroup at smc.vnet.net*Subject*: [mg57224] Re: [mg57198] Re: [mg57156] Re: [mg57137] How to get an answer as a Root object?*From*: Andrzej Kozlowski <akoz at mimuw.edu.pl>*Date*: Sat, 21 May 2005 02:39:45 -0400 (EDT)*References*: <200505170520.BAA25934@smc.vnet.net> <200505190708.DAA13114@smc.vnet.net> <acbec1a405051916553ad63c12@mail.gmail.com> <200505200843.EAA00645@smc.vnet.net> <8a391a1d99566d06d1861364e4190b82@mimuw.edu.pl>*Sender*: owner-wri-mathgroup at wolfram.com

On 20 May 2005, at 22:21, Andrzej Kozlowski wrote: > I am a little surprised that nobody wrote to point out that some > gremlins also got into method 3. In fact I like this solution best and > it is the only method that seems to me mathematically fully > satisfactory. I would dismiss method 2 as basically a hack relying on > the internals of the Mathematica implementation of Groebner basis: as > far as I know "mathematical" Groebner basis is not "supposed to" work > with expression of this kind. > Mathod 1 is of course classical (discussed for example in 1991 Henri > Cohen's "A Course in Computational Algebraic Number Theory" on page > 100) but I don't like the fact that you have to (or at least it seems > to me you have to) guess the degree of the polynomial used in > Recognize. I forgot to add that actually there can be no guarantee that the solution found with method 1 is the correct one; it could actually be a different algebraic number very close to the correct one. I believe that only the elimination method (whether accomplished by Groebner basis or Resultant) gives a solution that is guaranteed to be exact and does not rely on "accidental" properties of a particular implementation of Groebner basis. Andrzej > On the other hand method 3 is, I think, fully rigorous although it > does indeed require some "by hand working". Anyway, since some strange > errors also got into that code here is a corrected version: > > > 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]); > > 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}]]; > > Now here comes the point where something appears to have gone funny. > It should have been: > > expr2 = expr /. rep1 /. rep2 /. rep3 > > and finally as in the original version: > > > InputForm[Timing[poly4 = > First[GroebnerBasis[Join[repssq,{expr2}], s, > {rad0,rad1,rad2,rad3}, MonomialOrder->EliminationOrder]]]] > > > {1.0300000000000007*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} > > It is indeed clear that one can carry out the elimination process by > using Resultant wiht pairs of polynomials, starting with expr2 and > using the polynomials in repssq to eliminate the unneeded variables > one by one but this looks tedious compared with the Groebner basis > approach. > > Andrzej Kozlowski > > > > > On 20 May 2005, at 17:43, Daniel Lichtblau wrote: > >> *This message was transferred with a trial version of CommuniGate(tm) >> Pro* >> 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.706945837045705738269073202812169790793196895164013983001575249106 >>>> 155514352\ >>>> 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: Re: How to get an answer as a Root object?***From:*Daniel Lichtblau <danl@wolfram.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>

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