Services & ResourcesWolfram ForumsMathGroup Archive

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

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

• To: mathgroup at smc.vnet.net
• Subject: [mg68505] Re: [mg68479] Re: [mg68345] Re: Finding the Number of Pythagorean Triples below a bound
• From: Andrzej Kozlowski <akoz at mimuw.edu.pl>
• Date: Tue, 8 Aug 2006 06:28:49 -0400 (EDT)
• References: <20060806080123.15221.qmail@web50704.mail.yahoo.com> <18D1E833-7EBC-45EC-83A3-949C201D34EE@mimuw.edu.pl> <200608070540.BAA24073@smc.vnet.net>
• Sender: owner-wri-mathgroup at wolfram.com

On 7 Aug 2006, at 07:40, Andrzej Kozlowski wrote:

>
> On 6 Aug 2006, at 13:33, Andrzej Kozlowski wrote:
>
>>
>> On 6 Aug 2006, at 10:01, Titus Piezas wrote:
>>
>>> Hello Dan and Andrzej,
>>>
>>>> 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.
>>>
>>> Yes, turns out the problem was equivalent to counting the number
>>> of representations of N as
>>>
>>> a^2+b^2 = N
>>>
>>> where 0<a<=b (I should have phrased it that way, but that is the
>>> benefit of high-sight) which I assume can be given by the function
>>> OrderedSumOfSquaresRepresentations[2,N]?  Hence the code you have
>>> given me, in essence, is, say for k = 5,
>>>
>>> Sum[OrderedSumOfSquaresRepresentations[2,c^2+5],{c,1,10^m-1}]
>>>
>>> roughly speaking? (With finer details of removing from
>>> consideration those c^2+5 that are of from 4m+3, etc.)
>>
>> Well, no, it is rather:
>>
>> Sum[Length[OrderedSumOfSquaresRepresentations[2,c^2+5]],{c,1,10^m-1}]
>>
>> In fact I forgot that there was OrderedSumOfSquaresRepresentations
>> as well as SumOfSquaresRepresentations, which is the main reason
>> why I did not think it would be very useful. Another reason is
>> that, I expected, computing Length of many lists and then adding
>> these Lengths up should be slower than just counting the total
>> number of elements.
>>
>>>
>>>> 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[94]:=Timing[countTriples[10^6,3]]
>>>> Out[94]={131.91 Second,121054}
>>> There is something odd with the first value. James' count (up to
>>> m=5) vs this gives the sequences,
>>>
>>> {1219, 12115}
>>> {1154, 12115, 121054}
>>>
>>> Why the disparity? If the ratio S(N)/N would converge, it should
>>> start as 0.121... (Note however, that James' table has bound 10^m
>>> *not* 10^m-1.}
>>
>> My code also gives:
>>
>>
>> countTriplesP[4,3]
>>
>> 1219
>>
>> Changing 10^m-1 to 10^m should give a value at least as large, and
>> in fact gives the same answer.
>>
>>
>> which suggests that something still needs to be looked carefully at
>> in Daniel's code. I can't do it now but will think about it when I
>> can.
>>
>> Andrzej
>
>
> I have looked carefully at Daniel's code and I found I needn't have
> bothered ;-
> On my PowerBook Daniel's code gives:
>
>
> countTriples[10^4,3]//Timing
>
> {0.94218 Second,1219}
>
> which shows that my computer is about half the speed of his but also
> that somehow he probably pasted the wrong output into his posting?
> Everything about the code seems perfect.
>
> In fact, I am sure we could probably make it even faster if we looked
> inside the code for  OrderedSumOfSquaresRepresentations and changed
> it so that it would count the number of representations rather than
> make a list of them. But I am not sure that the extra speed we would
> gain would be worth the extra effort.
>
> Andrzej Kozlowski
>


As Daniel pointed out, my last comment above was based on forgetting  
how fast Length is in Mathematica - now "I am sure that we  
probably" ;-) could not make the code significantly faster by  
tinkling with OrderedSumOfSquaresRepresentations.

However, there are two additional comments, that I would like to  
make, which I think are of some interest. First, I think I know how  
Daniel got his wrong answer

In[54]:=Timing[countTriples[10^4,3]]
Out[54]={0.591 Second,1154}

If you forget to load the package

<< NumberTheory`NumberTheoryFunctions`

And run Daniel's code above you will get the answer 1154. This is  
what must have happened.

Secondly, I think I have found another speed improvement to Daniel's  
code. This is my new version:

countTriples1[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 = If[PrimeQ[c2odd], 1, Length 
[OrderedSumOfSquaresRepresentations[2, \
c2odd]]];
       total += add;
       , {c, 1, m - 1}];
     total]

With Daniel's code I get on my PowerBook:

In[4]:=
countTriples[10^4,3]//Timing

Out[4]=
{2.23016 Second,1219}

With the new code I get:

In[5]:=
countTriples1[10^4,3]//Timing

Out[5]=
{0.728627 Second,1219}

Which is more than 3 times faster. The improvement is again based on  
a bit of mathematics and a some powerful Mathematica technology. The  
mathematics is that it is known and not hard to prove that a prime  
number  the form 4m+1 there has exactly one representation as a sum  
of two squares. Of course to use this fact we need PrimeQ to be as  
fast as it is; there is no way you could use this fact without having  
a function like that at your disposal.

It seems the code is now pretty close to the optimum.

Andrzej Kozlowski

• References:
  • Re: Re: Finding the Number of Pythagorean Triples below a bound
    • From: Andrzej Kozlowski <akoz@mimuw.edu.pl>