Mathematica 9 is now available
Services & Resources / Wolfram Forums
-----
 /
MathGroup Archive
2006
*January
*February
*March
*April
*May
*June
*July
*August
*September
*October
*November
*December
*Archive Index
*Ask about this page
*Print this page
*Give us feedback
*Sign up for the Wolfram Insider

MathGroup Archive 2006

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

Search the Archive

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


  • Prev by Date: Re: Re: Finding the Number of Pythagorean Triples below a bound
  • Next by Date: Re: Re: Re: Re: Finding the Number of Pythagorean Triples below a bound
  • Previous by thread: Re: Re: Finding the Number of Pythagorean Triples below a bound
  • Next by thread: Re: x=2;Composition[f,FindMinimum][x+1,{x,a}]