Services & Resources / Wolfram Forums / MathGroup Archive

MathGroup Archive 2005

[Date Index] [Thread Index] [Author Index]

Search the Archive

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

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. 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.7069458370457057382690732028121697907931968951640139830015752491061 
>>> 55514352\
>>> 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
>
  • Prev by Date: Re: four approaches to do a simple sum
  • Next by Date: Re: Re: Re: How to get an answer as a Root object?
  • Previous by thread: Re: Re: How to get an answer as a Root object?
  • Next by thread: Re: Re: Re: How to get an answer as a Root object?