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

*To*: mathgroup at smc.vnet.net*Subject*: [mg68431] Re: [mg68382] Re: [mg68345] Re: Finding the Number of Pythagorean Triples below a bound*From*: Andrzej Kozlowski <akoz at mimuw.edu.pl>*Date*: Sat, 5 Aug 2006 03:46:14 -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> <79C36C70-E091-4A82-8EC5-0EDC743D081D@mimuw.edu.pl> <200608031007.GAA15743@smc.vnet.net> <76C59C34-5A55-4318-B0AF-3A572E71B421@mimuw.edu.pl> <C56DFF0B-6F0B-4E64-B72A-D805B3C6C063@mimuw.edu.pl> <7AA8804E-73F2-4D7B-BCAB-633CE33CD1D0@mimuw.edu.pl> <F7520FE0-16DB-4F8F-A39E-DF64749B85F0@mimuw.edu.pl>*Sender*: owner-wri-mathgroup at wolfram.com

I have good news: the code I just posted can be compiled and then it becomes really fast ;-) countT = Compile[{{m, _Integer}, {k, _Integer}}, Module[ {i = 0, c2, diff, sdiff}, Do [ If[Mod[c^2 + k, 4] != 3, c2 = Floor[Sqrt[(c^2 + k)/2]]; Do [ diff = c^2 + k - b^2; sdiff = Sqrt[N[diff]]; If [sdiff >= b && sdiff == Round[sdiff], i++]; , {b, 1, c2}]]; , {c, 1, 10^m - 1}]; i]] countT[3,2]+countT[3,-2]//Timing {1.0032 Second,282} I will try the case m=4 and m=5 and send you the results. I am not promising to do this very soon, just in case ;-) Andrzej Kozlowski On 4 Aug 2006, at 12:09, Andrzej Kozlowski wrote: > Here is the fastest code I have seen so far that works (it is based > on the one that Daniel sent you). I have corrected and enhanced > and enhanced it. It gives approximate answers for reasons that I > explained in earlier postings (use of machine arithmetic to test > for perfect integers). Of course it is easy to replace the code by > exact code by replacing the numerical test for a perfect square by > an exact one. > > > countTriples[m_,k_] := Module[ > {i=0, c2, diff, sdiff}, > Do [ > If[Mod[c^2+k,4]!=3, c2 = Floor[Sqrt[(c^2+k)/2]]; > Do [ > diff = c^2+k-b^2; > sdiff = Sqrt[N[diff]]; > If [sdiff>=b&&sdiff==Round[sdiff],i++]; > , {b,1,c2}]]; > ,{c,1,10^m-1}]; > i] > > > countTriples[3,2]+countTriples[3,-2]//Timing > > > {12.3746 Second,282} > > The correct answer is 283. > > This code should easily deal with the case m=4 (I have not yet > tried it) and I think even m=5 should now be within reach. > > Andrzej Kozlowski > > > > On 4 Aug 2006, at 11:27, Andrzej Kozlowski wrote: > >> The "improvement" below which I sent a little earlier is wrong >> (even though it returned correct answers). Obviously the point is >> that c^2+k can't be of the form 4n + 3, but there is no reason >> why a^2+b^2-k can't be of that form. Since my code does not >> explicitly select c it can't make use of this additional >> improvement. A different code, which uses explicit choices of >> (say) a and c and tests for b being a perfect square could exploit >> this fact and perhaps gain extra speed. It should not be difficult >> to write such a code along the lines I have been using. >> >> Andrzej >> >> >> On 4 Aug 2006, at 10:47, Andrzej Kozlowski wrote: >> >>> >>> On 3 Aug 2006, at 16:42, Andrzej Kozlowski wrote: >>> >>>> >>>> On 3 Aug 2006, at 12:07, Andrzej Kozlowski wrote: >>>> >>>>> On 2 Aug 2006, at 20:01, Andrzej Kozlowski wrote: >>>>> >>>>>> >>>>>> On 2 Aug 2006, at 11:23, titus_piezas at yahoo.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 >>>>>>> >>>>>> >>>>>> >>>>>> Here is a piece code which utilises the ideas I have described in >>>>>> my previous posts: >>>>>> >>>>>> ls = Prime /@ Range[3, 10]; >>>>>> >>>>>> test[n_] := >>>>>> Not[MemberQ[JacobiSymbol[n, ls], -1]] && Element[Sqrt[n], >>>>>> Integers] >>>>>> >>>>>> f[P_, k_] := Sum[If[(w = >>>>>> a^2 + b^2 - k) < P^2 && test[w], 1, 0], {a, 1, P}, {b, a, >>>>>> Floor[Sqrt[P^2 - a^2]]}] >>>>>> >>>>>> g[m_, k_] := f[10^m, k] + f[10^m, -k] >>>>>> >>>>>> We can easily confirm the results of your computations,.e.g. >>>>>> for k=2. >>>>>> >>>>>> >>>>>> Table[g[i,2],{i,3}] >>>>>> >>>>>> >>>>>> {4,30,283} >>>>>> >>>>>> >>>>>> Since you have not revealed the "simple code" you have used it is >>>>>> hard to tell if the above one is any better. It is however, >>>>>> certainly capable of solving the problem for m=4: >>>>>> >>>>>> >>>>>> g[4,2]//Timing >>>>>> >>>>>> >>>>>> {4779.39 Second,2763} >>>>>> >>>>>> The time it took on my 1 Gigahertz PowerBook was over 70 minutes, >>>>>> which is longer than you thought "reasonable", so I am still not >>>>>> sure if this is any improvement on what you already have. The >>>>>> time >>>>>> complexity of this algorithm seems somewhat larger than >>>>>> exponential >>>>>> so I would expect that it will take about 6 hours to deal with >>>>>> n=5 >>>>>> on my computer, and perhaps 2 weeks to deal with n=6. >>>>>> >>>>>> Andrzej Kozlowski >>>>>> >>>>>> >>>>>> >>>>> >>>>> >>>>> I mistakenly copied and pasted a wrong (earlier) definition of f. >>>>> Here is the correct one: >>>>> >>>>> f[P_, k_] := Block[{w}, Sum[If[GCD[a, b, k] == >>>>> 1 && (w = a^2 + >>>>> b^2 - k) < P^2 && test[w], 1, 0], {a, 1, >>>>> P}, {b, a, Floor[Sqrt[P^2 - a^2]]}]] >>>>> >>>>> The definition of g is as before: >>>>> >>>>> g[m_, k_] := f[10^m, k] + f[10^m, -k] >>>>> >>>>> Andrzej Kozlowski >>>>> >>>> >>>> >>>> Below is a faster version of the above code. (It owes a >>>> significant improvement to Daniel Lichtblau, which I borrowed >>>> from him without his knowledge ;-)) >>>> >>>> test[n_] := >>>> JacobiSymbol[n, >>>> Prime[Random[Integer, {2, 20}]]] -1 && Element[Sqrt[n], >>>> Integers] >>>> >>>> >>>> f[P_, k_] := Block[{w}, Sum[If[GCD[a, b, k] == >>>> 1 && (w = a^2 + >>>> b^2 - k) < P^2 && test[w], 1, 0], {a, 1, >>>> Floor[Sqrt[(P^2+k)/2]]}, {b, a, Floor[Sqrt[P^2 - a^2 >>>> +k]]}]] >>>> >>>> g[m_, k_] := f[10^m, k] + f[10^m, -k] >>>> >>>> The improvement is in the upper bound on a in the sum. Since a >>>> is the smaller of the two squares whose sum is equal to P^2+k it >>>> can't be larger than Floor[Sqrt[(P^2+k)/2]]. >>>> >>>> Note that you can improve the performance by loosing some >>>> accuracy if you use a cruder test for a perfect square: >>>> >>>> test1[n_] := With[{w = Sqrt[N[n]]}, w == Round[w] >>>> ] >>>> >>>> f1[P_, k_] := Block[{w}, Sum[If[GCD[a, b, k] == >>>> 1 && (w = a^2 + >>>> b^2 - k) < P^2 && test1[w], 1, 0], {a, 1, >>>> Floor[Sqrt[(P^2+k)/2]]}, {b, a, Floor[Sqrt[P^2 - a^2 >>>> +k]]}]] >>>> >>>> >>>> g1[m_, k_] := f1[10^m, k] + f1[10^m, -k] >>>> >>>> >>>> Let's compare the two cases. >>>> >>>> In[7]:= >>>> g[3,2]//Timing >>>> >>>> Out[7]= >>>> {89.554 Second,283} >>>> >>>> In[8]:= >>>> g1[3,2]//Timing >>>> >>>> Out[8]= >>>> {37.376 Second,283} >>>> >>>> So we see that we get the same answer and the second approach is >>>> considerably faster. However: >>>> >>>> In[9]:= >>>> g[1,6]//Timing >>>> >>>> Out[9]= >>>> {0.008863 Second,3} >>>> >>>> In[10]:= >>>> g1[1,6]//Timing >>>> >>>> Out[10]= >>>> {0.005429 Second,5} >>>> >>>> >>>> The correct answer is 3 (returned by the first method). The >>>> faster method found two false solutions. This should not matter >>>> if you are interested only in approximate answers (as you seem >>>> to be) but it is worth keeping in mind. >>>> >>>> Andrzej Kozlowski >>> >>> >>> I have noticed one more obvious improvement. We can replace test by: >>> >>> test[n_] := >>> Mod[n, 4] =!= 3 && JacobiSymbol[n, >>> Prime[Random[Integer, {2, 100}]]] -1 && Element[Sqrt[n], >>> Integers] >>> >>> and test1 by >>> >>> test1[n_] := With[{w = Sqrt[N[n]]}, Mod[n, 4] =!= 3 && w == Round[w] >>> ] >>> >>> We are simply making use of the easily to prove fact that an >>> integer of the form 4 k + 1 can never be the sum of two squares. >>> There is a noticeable improvement in the performance of g: >>> >>> >>> g[3,2]//Timing >>> >>> {58.0786 Second,283} >>> >>> However, the performance of g1 seems to actually decline slightly: >>> >>> >>> g1[3,2]//Timing >>> >>> {40.8776 Second,283} >>> >>> >>> However, there are fewer cases of "false solutions": >>> >>> In[22]:= >>> g1[1,6]//Timing >>> >>> Out[22]= >>> {0.006694 Second,4} >>> >>> >>> I am still not sure what is the most efficient use of >>> JacobiSymbol in this kind of problems. In the first version of my >>> code I used a test involving the first few odd primes, afterwards >>> I switched to using just one random prime in taken form a certain >>> range. Evaluation of JacobiSympol[m,p] != -1 is certainly much >>> faster than that of Element[Sqrt[m],Integers], but it is not >>> clear what is the most efficient number of primes to use, which >>> are the best primes and whether it is better to choose them at >>> random or just use a fixed selection. >>> The numerical test Sqrt[N[n]]==Round[Sqrt[N[n]] is even faster, >>> but it will sometimes produce "false squares". >>> >>> Andrzej Kozlowski >>> >>> >> >

**References**:**Re: Finding the Number of Pythagorean Triples below a bound***From:*titus_piezas@yahoo.com

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

**Re: Re: Re: Re: "No more memory available" -- a recurring problem**

**Re: Re: Re: Re: "No more memory available" -- a recurring problem**

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

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