Integer Solutions to Multi Variable Polynomials

[Eric Martin]

I'm new to Sympy, and I was wondering if there was anyway to find
integer solutions to multivariable polynomials. An example of this
polynomial would be 4x^2 + y^2 - 17 = 0 , or 4*x**2 + y**2 - 17 = 0 if
you prefer sympy notation. The integer solutions to this polynomial
are (2,1), (-2,1), (2,-1), and (-2, -1). Some polynomials do not have
integer solutions, for example if 18 was subtracted instead of 17 in
the example above.

Is there any way I can input these polynomials and sympy can return
the integer solutions?

[Aaron S. Meurer]

Sorry, I don't think this is implemented. Do you know what kind of algorithms would be required to solve these kinds of Diophantine problems? It would be nice to have something like Maple's isolve():

(in Maple)
> isolve(4*x**2 + y**2 - 17);
{x = -2, y = -1}, {x = -2, y = 1}, {x = 2, y = -1}, {x = 2, y = 1}

Aaron Meurer

> --
> You received this message because you are subscribed to the Google Groups "sympy" group.
> To post to this group, send email to sy...@googlegroups.com.
> To unsubscribe from this group, send email to sympy+un...@googlegroups.com.
> For more options, visit this group at http://groups.google.com/group/sympy?hl=en.
>

[Ondrej Certik]

On Wed, Jul 7, 2010 at 9:17 AM, Aaron S. Meurer <asme...@gmail.com> wrote:
> Sorry, I don't think this is implemented.  Do you know what kind of algorithms would be required to solve these kinds of Diophantine problems?  It would be nice to have something like Maple's isolve():
>
> (in Maple)
>> isolve(4*x**2 + y**2 - 17);
>     {x = -2, y = -1}, {x = -2, y = 1}, {x = 2, y = -1}, {x = 2, y = 1}

I wrote some algorithm for that long time (at least 10 years) ago in
Delphi. It was pretty short.

I think I implemented something like this:

http://en.wikipedia.org/wiki/Extended_Euclidean_algorithm

imho it should not be that difficult to implement this again in sympy,
at least for some simple cases, like the one above.

Ondrej

[Aaron S. Meurer]

Well, we already have the Extended Euclidean Algorithm in sympy :)

In [1]: gcdex(12, 32)
Out[1]: (3, -1, 4)

But that is only useful for solving equations of the form a*x + b*x == c. I guess for this specific problem, you could see if the solutions to 4*x + 1*y == 17 are perfect squares and take +/- sqrt() of them. Also look at the diophantine version of gcdex, which I have implemented buried in my integration branch (see http://github.com/asmeurer/sympy/blob/integration/sympy/integrals/risch.py#L50; all it essentially does more than gcdex() is multiply through by c/gcd(a, b)).

We also have the Chinese Remainder Theorem buried in galiostools.py, and it shouldn't be too hard to write a high-level algorithm for it.

Solving Diophantine equations in general is undecidable (Hilbert's tenth problem), but clearly there are heuristics (I don't know what they are, though, other than these).

Aaron Meurer

[Eric Martin]

Thanks, I'll investigate into possibly using gcdex to implement/
develop an algorithm. I'm doing this to implement the Sieve of Atkin
(prime number finder, http://en.wikipedia.org/wiki/Sieve_of_atkin ),
and it requires the number of solutions to equations similar to what I
originally posted. Does anyone have better advice for how I can find
the number of solutions to those Diophantine equations?

On Jul 7, 5:14 pm, "Aaron S. Meurer" <asmeu...@gmail.com> wrote:
> Well, we already have the Extended Euclidean Algorithm in sympy :)
>
> In [1]: gcdex(12, 32)
> Out[1]: (3, -1, 4)
>

> But that is only useful for solving equations of the form a*x + b*x == c.  I guess for this specific problem, you could see if the solutions to 4*x + 1*y == 17 are perfect squares and take +/- sqrt() of them.  Also look at the diophantine version of gcdex, which I have implemented buried in my integration branch (seehttp://github.com/asmeurer/sympy/blob/integration/sympy/integrals/ris...all it essentially does more than gcdex() is multiply through by c/gcd(a, b)).  

[Aaron S. Meurer]

If it's just a*x**2 + b*y**2 == c, with a, b, and c given (it looks like it is), then the algorithm I described should be complete.

See also section 1.3 of this book [0] (page 13 in particular).

Aaron Meurer

P.S. I was thinking that we needed to implement the Sieve of Atkin instead of Eratosthenes just the other day.

[0] - http://books.google.com/books?id=8SAaSd89sSkC&printsec=frontcover&source=gbs_slider_thumb#v=onepage&q&f=false

[smichr]

On Jul 8, 8:59 am, "Aaron S. Meurer" <asmeu...@gmail.com> wrote:
> If it's just a*x**2 + b*y**2 == c, with a, b, and c given (it looks like it is), then the algorithm I described should be complete.
>

See also http://www.alpertron.com.ar/QUAD.HTM for integer solutions
to
A*x**2 + B*x*y + C*y**2 + D*x + E*y = 0. The answers, step-by-step and
code (java) are presented there.

[smichr]

> P.S. I was thinking that we needed to implement the Sieve of Atkin instead of Eratosthenes just the other day.  
>

You might check out http://www.daniweb.com/code/snippet216558.html#post1254825
and look for sieveOfAtkin lower on the page.

Also http://programmingpraxis.com/2010/02/19/sieve-of-atkin-improved
which refers to Daniel Bernstein's very fast primegen at
http://cr.yp.to/primegen.html

[smichr]

On Aug 31, 9:52 am, smichr <smi...@gmail.com> wrote:
> > P.S. I was thinking that we needed to implement the Sieve of Atkin instead of Eratosthenes just the other day.  
>

> You might check outhttp://www.daniweb.com/code/snippet216558.html#post1254825

> and look for sieveOfAtkin lower on the page.

That version is 10X faster at generating primes to 1000000 than our
present primerange function.