MathGroup Archive 2005

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

Search the Archive

Re: Solving Diophantine Equations

  • To: mathgroup at
  • Subject: [mg61413] Re: Solving Diophantine Equations
  • From: "Ray Koopman" <koopman at>
  • Date: Tue, 18 Oct 2005 02:45:14 -0400 (EDT)
  • References: <> <> <> <dio1s4$suo$> <diprnq$hn2$> <divhhk$gb0$>
  • Sender: owner-wri-mathgroup at

Diana wrote:
> Ray,
> I am a beginning - intermediate Mathematica user.
> The algorithm you and Andrzej have coded works great. Could I trouble you to
> explain briefly what the code is doing? To the un-initiated, the
> sophisticated code is hard to interpret.
> As was mentioned in this thread, I am trying to corroborate the findings of
> Pingzhi Yuan in 2004.
> I hope to use your improvements towards the equations of twenty other
> Diophantine equation articles, as well, for a thesis.
> Thanks,
> Diana

Here's an improved gg, with some comments.

gg2[a_, b_, c_:5] := Block[{y, r1}, Reap[
Do[y = 2;
   While[y^(n-1) < (y - 1)a + y,
         r1 = (y^n - 1)/(y - 1) - 1;
         Do[If[(x+1)x == r1, Sow[{x, y, n}]],
            {x, y^((n-1)/2) + Boole[n==3], (y^(n-1) - y)/(y - 1)}];
   {n, c, b, 2}]][[2,1]]]

{6.91 Second,{{5,2,5},{90,2,13}}}

The arguments are
a = upper bound for x
b = upper bound for n
c = lower bound for n; defaults to 5

Only y & r1 are declared explicitly as local variables, because n & x
are created implicitly as local variables by the Do statements.

There are three nested loops. The outermost is Do[...,{n,c,b,2}].
Note that there is no check that c is odd.

The middle loop is y=2;While[y^(n-1)<(y-1)a+y,...;y++]. The final y
is the largest value for which the corresponding upper bound for x,
(y^(n-1) - y)/(y - 1), is less than a.

The innermost loop is Do[...,{x,xmin,xmax}],
with xmin = y^((n-1)/2) + Boole[n==3]
and xmax = (y^(n-1) - y)/(y - 1).
Note that xmin is adjusted so that the test x!=y can be omitted.

Inside the x-loop, Sow[{x,y,n}] saves {x,y,n} whenever
(x^3-1)/(x-1) == (y^n-1)/(y-1). Note that the much of the arithmetic
involved in the comparison is done outside the x-loop.
Note also that (x+1)x is faster than x*x+x.

Finally, Reap[...][[2,1]] returns what was Sown.
If no solutions were found then there will be an error message
"Part::partw: Part 1 of {} does not exist."
and the result will be {Null,{}}[[2,1]].
If that bothers you, change [[2,1]] to simply [[2]],
which will give the solutions with an extra level of nesting
and thus avoid the error message when there are no solutions.

  • Prev by Date: Re: How smooth graphs?
  • Next by Date: Re: Re: big integer for Range
  • Previous by thread: Re: Solving Diophantine Equations
  • Next by thread: Re: Re: Solving Diophantine Equations