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>