WOLFRAM

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?

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
>>
>
  • Prev by Date: Re: Re: Re: How to get an answer as a Root object?
  • Next by Date: Re: Nestwhile
  • Previous by thread: Re: Re: Re: How to get an answer as a Root object?
  • Next by thread: Re: Re: Re: Re: How to get an answer as a Root object?