Re: Solving Diophantine Equations
- To: mathgroup at smc.vnet.net
- Subject: [mg61413] Re: Solving Diophantine Equations
- From: "Ray Koopman" <koopman at sfu.ca>
- Date: Tue, 18 Oct 2005 02:45:14 -0400 (EDT)
- References: <32592565.1129236555309.JavaMail.root@elwamui-darkeyed.atl.sa.earthlink.net> <9F2DCCEC-34BF-4F53-9475-F686A911F260@akikoz.net> <E39C9EDB-A88A-42F3-9AD9-7DA09C7D790A@yhc.att.ne.jp> <dio1s4$suo$1@smc.vnet.net> <diprnq$hn2$1@smc.vnet.net> <divhhk$gb0$1@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
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)}]; y++], {n, c, b, 2}]][[2,1]]] gg2[10^5,30]//Timing {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.
- Follow-Ups:
- Re: Re: Solving Diophantine Equations
- From: Andrzej Kozlowski <akoz@mimuw.edu.pl>
- Re: Re: Solving Diophantine Equations