Re: Finding the Number of Pythagorean Triples below a bound
- To: mathgroup at smc.vnet.net
- Subject: [mg68232] Re: Finding the Number of Pythagorean Triples below a bound
- From: Peter Pein <petsie at dordos.net>
- Date: Sun, 30 Jul 2006 04:47:56 -0400 (EDT)
- References: <eaeqa3$53v$1@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
titus_piezas at yahoo.com schrieb:
> Hello all,
>
> Given the equation a^2+b^2 = c^2 with (a,b,c) positive integers,
> 0<a<=b. (For argument's sake, a=b will be allowed.) Let S(10^m) be the
> number of solutions with hypotenuse c < bound 10^m. (They need not be
> primitive.)
>
> Problem: Find S(10^m).
>
> What would be an efficient Mathematica code for this? The explicit
> (a,b,c)'s are not needed, just S(10^m). I also need "m" as high as
> possible, maybe m=6,7, with a reasonable run-time, say, less than an
> hour. (A kind person gave me a code but told me it may be practical
> only up to m=3.)
>
> I stumbled upon some code some time ago, but it only gives the explicit
> (a,b,c)'s with c<10^m, not the count:
>
> Block[{a,b}, For[a=1,a<(10^m/Sqrt[2]),a++, For[b=a,b<10^m,b++,
> If[Sqrt[a^2+b^2]==Round[Sqrt[a^2+b^2]]&&Sqrt[a^2+b^2]<10^m,
> Print[{a,b,Sqrt[a^2+b^2]}]]]]]
>
> for some positive integer "m". Can the above program be modified to
> give S(10^m)? Or what program would be fast, for m=6,7?
>
> Any help would be appreciated. Thanks.
>
> -Titus
>
Hi Titus,
using a huge amount of RAM, on can obtain quite short times:
See http://mathworld.wolfram.com/PythagoreanTriple.html, formula (11)
In[1]:= Reduce[{0 < z2^2 - z1^2 && 0 < 2*z1*z2 && z1^2 + z2^2 < 10^m && z1 > 0}, {z1, z2}]
Out[1]= 10^m > 0 && 0 < z1 < Sqrt[10^m]/Sqrt[2] && z1 < z2 < Sqrt[10^m - z1^2]
In[2]:= pythagoCount[m_] :=
Total[Floor[(10^m - 1/2)/Sow[Flatten[Table[(Pick[z^2 + #1^2, (GCD[z, #1] == 1 & ) /@ #1] & )[
Range[z + 1, Floor[Sqrt[10^m - z^2]], 2]], {z, Floor[Sqrt[10^m/2]]}], 1]]]]
In[3]:= {thetime, {allcount, primcount}} =
Timing[Transpose[(Through[{First, Length[hyp = #1[[2,1,1]]] & }[Reap[pythagoCount[#1]]]] & ) /@
Range[7]]]
Out[3]= {11.141*Second,
{{1, 50, 878, 12467, 161431, 1980636, 23471468}, {1, 16, 158, 1593, 15919, 159139, 1591579}}}
In[4]:= N[1/(2*Pi) - primcount/10^Range[7]]
Out[4]= {0.05915494309189534, -0.0008450569081046577, 0.001154943091895344, -0.00014505690810465155,
-0.00003505690810465256, 0.00001594309189534293, -2.9569081046454393*^-6}
For primcount see http://www.research.att.com/~njas/sequences/A101931 and for allcount see http://www.research.att.com/~njas/sequences/A101929.
HTH,
Peter