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

*To*: mathgroup at smc.vnet.net*Subject*: [mg68426] 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:02 -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>*Sender*: owner-wri-mathgroup at wolfram.com

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: Finding the Number of Pythagorean Triples below a bound**

**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**