MathGroup Archive 2006

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

Search the Archive

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

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

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},
       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, \
       total += add;
       , {c, 1, m - 1}];

With Daniel's code I get on my PowerBook:


{2.23016 Second,1219}

With the new code I get:


{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

  • Prev by Date: Re: Using Map with function that has options...
  • Next by Date: Re: Minimize
  • Previous by thread: Re: Re: Finding the Number of Pythagorean Triples below a bound
  • Next by thread: Two strange problems with a notebook...