Re: Re: Finding the Number of Pythagorean Triples below a bound
- To: mathgroup at smc.vnet.net
- Subject: [mg68463] Re: [mg68345] Re: Finding the Number of Pythagorean Triples below a bound
- From: danl at wolfram.com
- Date: Sun, 6 Aug 2006 02:56:37 -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>
- Sender: owner-wri-mathgroup at wolfram.com
> 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
- References:
- Re: Finding the Number of Pythagorean Triples below a bound
- From: titus_piezas@yahoo.com
- Re: Finding the Number of Pythagorean Triples below a bound