Re: Help with algorithm to find rational roots of a bivariate

*To*: mathgroup at smc.vnet.net*Subject*: [mg105187] Re: [mg105080] Help with algorithm to find rational roots of a bivariate*From*: Daniel Lichtblau <danl at wolfram.com>*Date*: Tue, 24 Nov 2009 05:45:55 -0500 (EST)*References*: <200911201139.GAA03433@smc.vnet.net>

TPiezas wrote: > Hello all, > > Does anyone know an efficient algorithm using Mathematica that can > find _rational_ roots of a non-homogenous eqn in two variables with > deg > 4? For ex, you want to find rational {x,y} such that, > > F(x,y) = x^n + (P_1)x^(n-1) + (P_2)x^(n-2) + .... = 0 > > where the P_i are polynomials in y. > > I _always_ come across this situation in the course of experimental > mathematics, and it would be great if Mathematica had a built-in > feature that solves _bivariate_ eqns in the rationals. Right now, I > have a 22-deg eqn in x with coefficients in y. I know three non- > trivial rational pairs {x,y} such that F(x,y) = 0, but I want to know > if there are others. If there are, a certain family of identities > would have more members. > > Any help will be appreciated. > > - Titus THis might give you cause to question that last remark. I can outline an approach that might get you rational solutions up to some specified size (in numerator and denominator), but it is not terribly efficient. The idea is to use homomorphic imaging, that is solve modulo a prime, and lift to prime powers. Say size is the larger in absolute value of numerator and denominator we will consider. In outline form, it goes like this. (1) First pick a modulus p of modest size, say 11. (2) For all 0<=x<=p-1, find solutions mod p for y. (3) Use rational recovery to go from solutions mod p to candidate solutions over the rationals. (4) Select from the candidates all actual solutions. (5) Replace p by p^2, and use Use Hensel lifting to take solutions mod the old p into solutions mod the new one. At this step an interesting thing happens: individual solutions can split into multiple solutions. This is bad in that it gives rise to a combinatorial explosion. But it is good because it allows us to find more rational solution candidates. (6) If the new p is larger than twice size^2, we are done. Else go to step (3). There are numerous caveats to this approach. First is that it is likely to be slow in actual applications of interest. Second is that there are ways in which a prime can be "unlucky". Obviously if a solution denominator is divisible by that prime, we'll not recover that rational. But I believe other things can go wrong, just as they can when using homomorphic images to factor a multivariate polynomial. A third is that I am not entirely certain the termination condition guarantees we'll have all solutions up to the specified size (assuming none were lost to the second caveat). I suspect it can be proven, in some manner similar to how one shows rational recovery gives guaranteed results for linear algebra problems after a certain amount of lifting (and assuming the prime is not one of finitely many that are unlucky...). A fourth, perhaps part of the first, is that I've no idea how to cull out lifted "solutions" that will not eventually give rise to actual solutions. I could remove the ones that gave solutions at an earlier lifting step, but in something like the example below this is but a small percentage. I'll illustrate with the usual Pythagorian triples problem, solving x^2+y^2-1 == 0 over the rationals. Note that I am using Mathematica's development version of Solve, so some amount of work would be needed to translate the current version (though this can be done). rationalRecover[x_,pk_]:= ((#[[2,2]]/#[[1,2,2]])&)[Internal`HGCD[pk,x]] bpoly = x^2+y^2-1; (* First we solve mod 11. *) firstsols = Flatten[Table[Solve[{bpoly,y-j}==0, {x,y}, Modulus->11], {j,2,10}] /. {}:>Sequence[], 1]; firstpairs = {x,y}/.firstsols; (* Now get candidates. *) firstrats = Map[rationalRecover[#,11^2]&, firstpairs, {2}]; (* Now lift and iterate the process, We'll harvest triples later. *) nextgbs = Table[GroebnerBasis[{11^2, bpoly, (x-firstpairs[[j,1]])^2, (y-firstpairs[[j,2]])^2}, {x,y}, CoefficientDomain->Integers], {j,Length[firstpairs]}]; nextsols = Flatten[Map[Solve[Rest[#]==0, {x,y}, Modulus->First[#]]&, nextgbs], 1]; nextpairs = {x,y}/.nextsols; nextrats = Map[rationalRecover[#,11^2]&, nextpairs, {2}]; (* Lift and iterate again. It starts to get slow at this point. I think that's the sound of a combinatorial explosion. *) nextgbs2 = Table[GroebnerBasis[{11^4, bpoly, (x-nextpairs[[j,1]])^2, (y-nextpairs[[j,2]])^2}, {x,y}, CoefficientDomain->Integers], {j,Length[nextpairs]}]; nextsols2 = Flatten[Map[Solve[Rest[#]==0, {x,y}, Modulus->First[#]]&, nextgbs2], 1]; nextpairs2 = {x,y}/.nextsols2; nextrats2 = Map[rationalRecover[#,11^4]&, nextpairs2, {2}]; (* Let's see what solutions we obtained. *) In[39]:= InputForm[firsttriples = Select[firstrats, (bpoly/.Thread[{x,y}->#])==0&]] Out[39]//InputForm= {} (* Okay nothing there. *) In[40]:= InputForm[nexttriples = Select[nextrats, (bpoly/.Thread[{x,y}->#])==0&]] Out[40]//InputForm= {{3/5, 4/5}, {-3/5, 4/5}, {4/5, 3/5}, {-4/5, 3/5}, {4/5, -3/5}, {-4/5, -3/5}, {3/5, -4/5}, {-3/5, -4/5}, {0, -1}} (* Something, though not terribly exciting. *) In[41]:= InputForm[nexttriples2 = Select[nextrats2, (bpoly/.Thread[{x,y}->#])==0&]] Out[41]//InputForm= {{-20/29, 21/29}, {-12/13, -5/13}, {3/5, 4/5}, {3/5, 4/5}, {3/5, 4/5}, {3/5, 4/5}, {3/5, 4/5}, {3/5, 4/5}, {3/5, 4/5}, {3/5, 4/5}, {3/5, 4/5}, {3/5, 4/5}, {3/5, 4/5}, {40/41, -9/41}, {-35/37, 12/37}, {-7/25, -24/25}, {8/17, -15/17}, {-16/65, 63/65}, {45/53, -28/53}, {-45/53, -28/53}, {16/65, 63/65}, {-8/17, -15/17}, {7/25, -24/25}, {35/37, 12/37}, {-40/41, -9/41}, {-3/5, 4/5}, {-3/5, 4/5}, {-3/5, 4/5}, {-3/5, 4/5}, {-3/5, 4/5}, {-3/5, 4/5}, {-3/5, 4/5}, {-3/5, 4/5}, {-3/5, 4/5}, {-3/5, 4/5}, {-3/5, 4/5}, {12/13, -5/13}, {20/29, 21/29}, {63/65, -16/65}, {4/5, 3/5}, {4/5, 3/5}, {4/5, 3/5}, {4/5, 3/5}, {4/5, 3/5}, {4/5, 3/5}, {4/5, 3/5}, {4/5, 3/5}, {4/5, 3/5}, {4/5, 3/5}, {4/5, 3/5}, {-28/53, 45/53}, {-9/41, 40/41}, {12/37, -35/37}, {21/29, -20/29}, {-24/25, -7/25}, {-5/13, -12/13}, {-15/17, 8/17}, {15/17, 8/17}, {5/13, -12/13}, {24/25, -7/25}, {-21/29, -20/29}, {-12/37, -35/37}, {9/41, 40/41}, {28/53, 45/53}, {-4/5, 3/5}, {-4/5, 3/5}, {-4/5, 3/5}, {-4/5, 3/5}, {-4/5, 3/5}, {-4/5, 3/5}, {-4/5, 3/5}, {-4/5, 3/5}, {-4/5, 3/5}, {-4/5, 3/5}, {-4/5, 3/5}, {-63/65, -16/65}, {63/65, 16/65}, {4/5, -3/5}, {4/5, -3/5}, {4/5, -3/5}, {4/5, -3/5}, {4/5, -3/5}, {4/5, -3/5}, {4/5, -3/5}, {4/5, -3/5}, {4/5, -3/5}, {4/5, -3/5}, {4/5, -3/5}, {-28/53, -45/53}, {-9/41, -40/41}, {12/37, 35/37}, {21/29, 20/29}, {-24/25, 7/25}, {-5/13, 12/13}, {-15/17, -8/17}, {15/17, -8/17}, {5/13, 12/13}, {24/25, 7/25}, {-21/29, 20/29}, {-12/37, 35/37}, {9/41, -40/41}, {28/53, -45/53}, {-4/5, -3/5}, {-4/5, -3/5}, {-4/5, -3/5}, {-4/5, -3/5}, {-4/5, -3/5}, {-4/5, -3/5}, {-4/5, -3/5}, {-4/5, -3/5}, {-4/5, -3/5}, {-4/5, -3/5}, {-4/5, -3/5}, {-63/65, 16/65}, {-20/29, -21/29}, {-12/13, 5/13}, {3/5, -4/5}, {3/5, -4/5}, {3/5, -4/5}, {3/5, -4/5}, {3/5, -4/5}, {3/5, -4/5}, {3/5, -4/5}, {3/5, -4/5}, {3/5, -4/5}, {3/5, -4/5}, {3/5, -4/5}, {40/41, 9/41}, {-35/37, -12/37}, {-7/25, 24/25}, {8/17, 15/17}, {-16/65, -63/65}, {45/53, 28/53}, {-45/53, 28/53}, {16/65, -63/65}, {-8/17, 15/17}, {7/25, 24/25}, {35/37, -12/37}, {-40/41, 9/41}, {-3/5, -4/5}, {-3/5, -4/5}, {-3/5, -4/5}, {-3/5, -4/5}, {-3/5, -4/5}, {-3/5, -4/5}, {-3/5, -4/5}, {-3/5, -4/5}, {-3/5, -4/5}, {-3/5, -4/5}, {-3/5, -4/5}, {12/13, 5/13}, {20/29, -21/29}, {0, -1}, {0, -1}, {0, -1}, {0, -1}, {0, -1}, {0, -1}, {0, -1}, {0, -1}, {0, -1}, {0, -1}, {0, -1}, {-99/101, 20/101}, {11/61, 60/61}, {-33/65, 56/65}, {-55/73, 48/73}, {77/85, 36/85}, {-77/85, 36/85}, {55/73, 48/73}, {33/65, 56/65}, {-11/61, 60/61}, {99/101, 20/101}} So that was a more substantial harvest. Clearly I skimped on a many details of why this might work at all. In part this is because I don't know them all. But here are a few places one can get background information on some of the tactics indicated above. For details on the upcoming Solve functionality: http://library.wolfram.com/infocenter/Conferences/7552/ For some idea of how the lifting works, at least in the univariate case: http://library.wolfram.com/infocenter/Conferences/7511/ http://library.wolfram.com/infocenter/MathSource/7522/ These show how one can use Groebner bases over the integers to lift p-adic gcds in the univariate case. I do not claim that the same proof gives lifted solutions in the bivariate case. But we do indeed get them by that method. For explanation of Rational Recovery, see the ARRA bill signed into law last February. No, wait, that was a different recovery (though perhaps rational all the same). Instead see: http://library.wolfram.com/infocenter/Conferences/7534/ Daniel Lichtblau Wolfram Research

**References**:**Help with algorithm to find rational roots of a bivariate equation?***From:*TPiezas <tpiezas@gmail.com>

**Re: Is it possible to make NIntegrate faster?**

**Re: Help with algorithm to find rational roots of a**

**Re: Help with algorithm to find rational roots of a bivariate equation?**

**Re: Help with algorithm to find rational roots of a**