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

*To*: mathgroup at smc.vnet.net*Subject*: [mg68441] 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:42 -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> <B45555BF-3093-4487-A3EA-0981783CA75A@mimuw.edu.pl> <B935EC43-291F-4C01-B646-534789B6E908@mimuw.edu.pl>*Sender*: owner-wri-mathgroup at wolfram.com

I need to add some corrections. There were some some small mistakes in the code I posted earlier, which caused problems with running the compiled code for negative values of k, and which also probably accounted for the slightly incorrect answers which the "approximate" code returned. I attributed the inaccuracy to the use of numerical precision in the test for a number being a perfect square but it seems that the cause was (probably) elsewhere. I don't want to try to make this message too long so I won't bother explaining what I think the mistakes were; but I will give what I think is the correct code. I have decided to separate the code for negative and positive values of k. The code for negative k works also for positive k's but is slightly slower, due to the extra test that needs to be performed. For this reason, and for the sake of greater clarity I have decided to separate the two codes. The code for positive k: countTriplesP = 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]] The code for negative k: countTriplesN = Compile[{{m, _Integer}, {k, _Integer}}, Module[ {i = 1, c2, diff, sdiff}, Do [ If[Mod[c^2 + k, 4] != 3 && c^2 + k ³ 0, 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]] Now we get: countTriplesP[1,6]+countTriplesN[1,-6]//Timing {0.000221 Second,3} countTriplesP[3,2]+countTriplesN[3,-2]//Timing {0.95186 Second,282} countTriplesP[4,2]+countTriplesN[4,-2]//Timing {95.2177 Second,2762} Note that these values are consistently less by one form the values obtained by Titus and also by me using my earlier "exact" code, but actually I believe that this disparity was due to mistakes in the earlier code. In any case, if we replace the "numerical" test for a perfect square by the much slower "exact" test, the answers will be the same, so the difference of 1 is certainly not due to the use of numerical precision. Anyway, everything works fast and looks perfectly satisfactory but then there is a snag. I decided to run the code for m=5 and k=2, started it and went out for several hours. When I came back I was rather disappointed to see that it was still running and then I saw the message: countTriplesP[5,2]//Timing CompiledFunction::"cfn" Numerical error encountered at instruction 10; proceeding with uncompiled evaluation. I assume the cause of this is that on 32 bit computers Compile cannot deal with integers larger than 2^32 but we have: c=10^5; Log[2,c^2]//N 33.2193 I can't at the moment see any way around this problem except to run uncompiled code (far too long) or try a 64 bit computer. Unfortunately I do not have one and I don't think Titus has one, so if any body has a 64-bit computer with a 64 bit version of Mathematica installed, and a little spare time, I think both of us would like to know how the above code performs in this setting. If it works it will provide a nice bit of advertising for 64 bit computing. Andrzej Kozlowski On 4 Aug 2006, at 12:48, Andrzej Kozlowski wrote: > My latest code has just managed: > > countT[4,2]+countT[4,-2]//Timing > > > {93.9638 Second,2762} > > > That is less than two minutes on a 1 gigahertz computer. The > correct answer is actually 2763 (by my earlier computation using > the exact test) so we have lost one solution but gained more than > 50 fold improvement in performance! > The case m=5 is now certainly feasible, although I am not sure if I > wish my not very powerful PowerBook to be occupied for so long, as > I need to use Mathematica for other tasks. Perhaps I can now leave > this to others. > > Andrzej Kozlowski > > On 4 Aug 2006, at 12:38, Andrzej Kozlowski wrote: > >> 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>