Mathematica 9 is now available
Services & Resources / Wolfram Forums / MathGroup Archive
-----

MathGroup Archive 2009

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

Search the Archive

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


  • Prev by Date: Re: Is it possible to make NIntegrate faster?
  • Next by Date: Re: Help with algorithm to find rational roots of a
  • Previous by thread: Re: Help with algorithm to find rational roots of a bivariate equation?
  • Next by thread: Re: Help with algorithm to find rational roots of a