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
- Re: Finding the Number of Pythagorean Triples below a bound