Re: Finding the Number of Pythagorean Triples below a bound
- To: mathgroup at smc.vnet.net
- Subject: [mg68283] Re: [mg68189] Finding the Number of Pythagorean Triples below a bound
- From: Andrzej Kozlowski <akoz at mimuw.edu.pl>
- Date: Mon, 31 Jul 2006 03:45:52 -0400 (EDT)
- References: <200607290500.BAA04928@smc.vnet.net> <52E622AC-6781-4ABD-8BFD-ED66D11D2FB7@mimuw.edu.pl>
- Sender: owner-wri-mathgroup at wolfram.com
On 30 Jul 2006, at 11:51, Andrzej Kozlowski wrote:
> On 29 Jul 2006, at 07:00, titus_piezas at yahoo.com wrote:
>
>> 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
>>
>
> Let me start by describing two possible way of attempting to solve
> your problem with the help of Mathematica. Afterwords I make a few
> remarks about the feasibility of solving your problem in this way
> for large values of m.
> I will describe two methods: a slow one and a much faster one. The
> reason for giving two methods is that in this way we can be more
> confident that the answers we are getting are correct ones. Even
> the slow method will be vastly faster than the one included in your
> post.
>
> The first (slow) method makes use of the function
> SumOfSquaresRepresentations , which is defined in the package
> NumberTheory`NumberTheoryFunctions`. First we load the package:
>
> << NumberTheory`NumberTheoryFunctions`
>
>
> ?SumOfSquaresRepresentations
>
> SumOfSquaresRepresentations[d, n] gives the list of all \
> representations of the positive integer n as a sum of d squares of \
> integers.
>
> For example:
>
> SumOfSquaresRepresentations[2, 50^2]
>
>
> {{-50, 0}, {-48, -14}, {-48, 14},
> {-40, -30}, {-40, 30},
> {-30, -40}, {-30, 40},
> {-14, -48}, {-14, 48},
> {0, -50}, {0, 50}, {14, -48},
> {14, 48}, {30, -40}, {30, 40},
> {40, -30}, {40, 30}, {48, -14},
> {48, 14}, {50, 0}}
>
> As you see, this is not quite right, as it includes non-positive
> integers as well as the same representation more then once. With
> Mathematica it is easy to convert this to the form you want:
>
> pythagoreans1[m_] := Module[{g}, g[n_] := Cases
> [SumOfSquaresRepresentations[2,
> n], {x_?Positive, y_?Positive} /; x > y]; Flatten[DeleteCases
> [g /@
> Table[i^2, {i, 1, m}], {}], 1]]
>
>
> Actually this is not quite the form you wanted because pythagoreans
> [m] finds all pythagorean pairs (x,y) with x^2+y^2<=m^2. However, I
> find this form most convenient so I will continue using it here. It
> can easily be converted to the form you wanted. So we have:
>
>
> pythagoreans1[50]
>
>
> {{4, 3}, {8, 6}, {12, 5}, {12, 9}, {15, 8}, {16, 12},
> {20, 15}, {24, 7}, {24, 10}, {21, 20}, {24, 18},
> {30, 16}, {28, 21}, {35, 12}, {36, 15}, {32, 24},
> {40, 9}, {36, 27}, {40, 30}, {48, 14}}
>
> If we do not wish to see the output but only compute the number of
> representations we can use the function NumberOFPythagoreans1
> defined as follows:
>
> NumberOfPythagoreans1[m_] := Length[pythagoreans1[m]]
>
> So that:
>
>
> NumberOfPythagoreans1[50]
>
> 20
>
>
> The method is inefficient and to make it much faster we shall use
> some elementary number theory. We know that (e.g., see Hardy and
> Wright, "An Introduction to The Theory Of Numbers) that all the
> Pythagorean triples (a,b,c) with coprime a and b have the form:
> {2*n*m, m^2 - n^2, m^2 + n^2} where n and m are coprime positive
> integers and m>n. So we first implement a function
> CoprimePythagoreans as follows:
>
>
> CoprimePythagoreans[p_] := DeleteCases[
> With[{q = Sqrt[p]}, Flatten[
> Table[If[GCD[n, m] == 1, {2*n*m, m^2 - n^2,
> m^2 + n^2}], {n, 1, q}, {m, n + 1,
> Sqrt[q^2 - n^2], 2}], 1]], Null]
>
> this computes coprime puthagorean triples in which the "hypotenuse"
> is smaller than p, e.g.
>
>
> CoprimePythagoreans[50]
>
>
> {{4, 3, 5}, {8, 15, 17}, {12, 35, 37}, {12, 5, 13},
> {20, 21, 29}, {24, 7, 25}, {40, 9, 41}}
>
>
> We can now compute the total number of all pythagorean triples with
> hypotenuse less than 50 without actually having to list them:
>
> NumberOfPythagoreans2[k_] := Total[Floor[k/#] & /@
> CoprimePythagoreans[k][[All, 3]]]
>
>
> NumberOfPythagoreans2[50]
>
> 20
>
>
> Happily this agrees with the previous answer. However, the
> difference in performance is pretty huge:
>
>
>
> NumberOfPythagoreans1[100]//Timing
>
>
> {0.427973 Second,52}
>
> v.s.
>
> NumberOfPythagoreans2[100]//Timing
>
>
> {0.005594 Second,52}
>
>
> In spite of its impressive speed I doubt that NumberOfPythagoreans2
> can cope with your problem for values of m as large as you want (I
> have not tried testing it). The problem is not actually speed but
> memory. The number of triples that have to be stored is a huge one,
> and probably will stretch or exceed the capabilities of even a
> powerful computer. This is, of course, something you can test
> yourself. The alternative approach would be to try to work out the
> number of pythagorean triples without actually storing them as a
> list but this is, of course, an entirely different type of problem.
>
> Andrzej Kozlowski
When I wrote "an entirely different type of problem" I had, of
course, in mind, an entirely mathematical solution. But almost
immediately after posting my message I realized that there is a very
simple computer solution based on the idea above, that avoids the
storage difficulty altogether. Here it is:
NumberOfPythagoreans[p_] := With[{q = Sqrt[p]}, Sum[If[GCD[n, m] ==
1, Floor[p/(n^2 + m^2)], 0], {n, 1, q},
{m, n + 1, Sqrt[q^2 - n^2], 2}]]
Now,
NumberOfPythagoreans[1000]//Timing
{0.010012 Second,881}
and
NumberOfPythagoreans[10^6]//Timing
{4.70295 Second,1980642}
Which seems to me to solve your problem completely.
Andrzej Kozlowski
- References:
- Finding the Number of Pythagorean Triples below a bound
- From: titus_piezas@yahoo.com
- Finding the Number of Pythagorean Triples below a bound