Services & Resources / Wolfram Forums
-----
 /
MathGroup Archive
2005
*January
*February
*March
*April
*May
*June
*July
*August
*September
*October
*November
*December
*Archive Index
*Ask about this page
*Print this page
*Give us feedback
*Sign up for the Wolfram Insider

MathGroup Archive 2005

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

Search the Archive

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

  • To: mathgroup at smc.vnet.net
  • Subject: [mg57156] Re: [mg57137] How to get an answer as a Root object?
  • From: Daniel Lichtblau <danl at wolfram.com>
  • Date: Thu, 19 May 2005 03:08:33 -0400 (EDT)
  • References: <200505170520.BAA25934@smc.vnet.net>
  • Sender: owner-wri-mathgroup at wolfram.com

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



  • Prev by Date: Re: UpSetDelayed and N
  • Next by Date: Errors from FindFit
  • Previous by thread: How to get an answer as a Root object?
  • Next by thread: Re: Re: How to get an answer as a Root object?