Re: Re: Finding the Number of Pythagorean Triples below a bound

*To*: mathgroup at smc.vnet.net*Subject*: [mg68469] Re: [mg68345] Re: Finding the Number of Pythagorean Triples below a bound*From*: Andrzej Kozlowski <akoz at mimuw.edu.pl>*Date*: Sun, 6 Aug 2006 02:57:01 -0400 (EDT)*References*: <eaeqa3$53v$1@smc.vnet.net><200607300848.EAA25171@smc.vnet.net> <eakfgm$rl6$1@smc.vnet.net> <200608020923.FAA28520@smc.vnet.net> <1484.74.136.205.171.1154809948.squirrel@webmail.wolfram.com>*Sender*: owner-wri-mathgroup at wolfram.com

On 5 Aug 2006, at 22:32, danl at wolfram.com wrote: >> Hello all, >> >> My thanks to Peter and Andrzej, as well as those who privately >> emailed >> me. >> >> To recall, the problem was counting the number of solutions to a >> bivariate polynomial equal to a square, >> >> Poly(a,b) = c^2 >> >> One form that interested me was the Pythagorean-like equation: >> >> a^2 + b^2 = c^2 + k >> >> for {a,b} a positive integer, 0<a<=b, and k any small integer. I was >> wondering about the density of solutions to this since I knew in the >> special case of k=0, let S(N) be the number of primitive solutions >> with >> c < N, then S(N)/N = 1/(2pi) as N -> inf. >> >> For k a squarefree integer, it is convenient that any solution is >> also >> primitive. I used a simple code that allowed me to find S(10^m) with >> m=1,2,3 for small values of k (for m=4 took my code more than 30 mins >> so I aborted it). The data is given below: >> >> Note: Values are total S(N) for *both* k & -k: >> >> k = 2 >> S(N) = 4, 30, 283 >> >> k = 3 >> S(N) = 3, 41, 410 >> >> k = 5 >> S(N) = 3, 43, 426 >> >> k = 6 >> S(N) = 3, 36, 351 >> >> Question: Does S(N)/N for these also converge? For example, for the >> particular case of k = -6, we have >> >> S(N) = 2, 20, 202 >> >> which looks suspiciously like the ratio might be converging. >> >> Anybody know of a code for this that can find m=4,5,6 in a reasonable >> amount of time? >> >> >> Yours, >> >> Titus > > Having been involved in some of the private email I have a bit of an > advantage. > > I'd say up to 10^6 or maybe 10^7 the way to go is by looking at the > factorization of c2=c^2+k. This way you only need iterate over c. > Also as > Andrzej Kozlowski pointed out, one can remove from consideration > all c2 > for which, after removing factors of 2, the remaining odd is 3 > modulo 4. I > use a prescan to enforce that if factors of 3,7, or 11 are present > then so > are their squares (note: I did not check whether this helps for > speed but > I suspect it does for large n). The full FactorInteger check > enforces that > any factor equal to 3 mod 4 appears to even degree. After that I > punt to > NumberTheory`NumberTheoryFunctions` to use > OrderedSumOfSquaresRepresentations. Finally, I found it faster to > use this > on c2 with powers of 2 removed. One can show this gives the correct > result. > > countTriples[m_, k_] := Module[ > {c2, c2odd, total = 0, fax, add, g}, > Do[ > c2 = c^2 + k; > If[c2 < 5, Continue[]]; > c2odd = c2; > While[EvenQ[c2odd], c2odd /= 2]; > If[Mod[c2odd, 4] == 3, Continue[]]; > g = GCD[c2odd, 231]; > If[g ª 1 && g^2 ª GCD[c2odd, 53361], Continue[]]; > fax = Mod[FactorInteger[c2odd], 4]; > If[Apply[Or, Map[#[[1]] == 3 && OddQ[#[[2]]] &, fax]], > Continue[]]; > add = OrderedSumOfSquaresRepresentations[2, c2odd]; > total += Length[add]; > , {c, 1, m - 1}]; > total] > > Quick test for agreement with your results: > > In[96]:=countTriples[10^3,-5]+countTriples[10^3,5] > Out[96]=426 > > Here are some timings. > > In[54]:=Timing[countTriples[10^4,3]] > Out[54]={0.591 Second,1154} > > In[93]:=Timing[countTriples[10^5,3]] > Out[93]={8.562 Second,12115} > > In[95]:=Timing[countTriples[10^5,2]] > Out[95]={24.054 Second,9912} > > In[94]:=Timing[countTriples[10^6,3]] > Out[94]={131.91 Second,121054} > > I would think it will handle 10^7 reasonably well (maybe an hour or > two). > But that depends on notgetting bogged in too many relatively expensive > factorizations. > > As I noted in email, I would not be surprised if there are further > improvements, maybe by avoiding the OrderedSumOfSquaresRepresentations > usage. > > Daniel Lichtblau > Wolfram Research > > > > Congratulations, this is rather spectacular, I think. I have to confess, I was very sceptical about the prospects of this method; now I think it was because I forgot how much good mathematics has gone into making FactorInteger as fast as it is! Actually, it is vastly faster not only than my compiled code but also than the C code that can be found here: http://pat7.com/jp/221k.c Which only goes to show, what all of us here have always known anyway ;-) namely, that to solve mathematical problems it is better to use mathematical methods and Mathematica than general purpose programming languages like C. Andrzej Kozlowski

**References**:**Re: Finding the Number of Pythagorean Triples below a bound***From:*titus_piezas@yahoo.com