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

*To*: mathgroup at smc.vnet.net*Subject*: [mg68468] Re: [mg68442] Re: [mg68382] Re: [mg68345] Re: Finding the Number of Pythagorean Triples below a bound*From*: Andrzej Kozlowski <akoz at mimuw.edu.pl>*Date*: Sun, 6 Aug 2006 02:56:58 -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> <F7F7DB69-7894-466A-99BC-065835567D5C@mimuw.edu.pl> <200608050746.DAA14183@smc.vnet.net> <71350A02-470E-4AFE-87B5-DCFACFC6483B@mimuw.edu.pl>*Sender*: owner-wri-mathgroup at wolfram.com

A small correction: 5136.24 seconds is about 1 hour and 25 minutes; hence the total time taken was about 1 hour 26.5 minute, not 1 hour 45 minutes as I wrote earlier. Andrzej Kozlowski On 5 Aug 2006, at 15:40, Andrzej Kozlowski wrote: > I have succeeded in computing the number of solutions of a^2 > +b^2==c^2+2 where c<10^5. This was done by means of the follwoign > code: > > test1 = 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, 10.^(m - 1), 10.^m - 1}]; > i]] > > I got > > > test1[5,2]//Timing > > > {5136.24 Second,8929} > > To this one should add the outcome of > > countTriplesP[4,2] > > 983 > > where countTriplesP is defined in my earlier post. So the actual > number of solutions is > > > 8929+983 > > 9912 > This means that compiled Mathematica code can quite successfully > deal with the case m=5 of Titus's problem. This computation took 1 > hour and 43 minutes on a 1 Gigahertz PowerBook G4, to which we have > to add about 1.5 minute for solving the 10^4 case. That gives the > total time of 1.45 minutes, which seems to me pretty impressive, > and shows, I think, that Mathematica is quite usable for hard > numerical computations, provided one one is lucky enough to have > ones code compile properly. I really hope some work at WRI is going > into making Compile work better (in more cases and more predictably). > > > Andrzej Kozlowski > > On 5 Aug 2006, at 09:46, Andrzej Kozlowski wrote: > >> Daniel Lichtblau has already informed me that running the code below >> on a 64 bit computer will not help at all, so there is no point >> trying :-( >> However, he also suggested a way around the Compile problem with >> integers larger than machine integers :-) >> I am trying it out right now. >> >> Andrzej Kozlowski >> >> >> On 4 Aug 2006, at 21:09, Andrzej Kozlowski wrote: >> >>> 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>

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

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

**NDSolve with boundary condition(s) at infinity?**

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

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