Prev: how to find a function
Next: problem about calling Mathematica in Fortran
From: Daniel Lichtblau on 24 Nov 2009 06:10
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

 | 
Pages: 1
Prev: how to find a function
Next: problem about calling Mathematica in Fortran