Services & ResourcesWolfram ForumsMathGroup Archive

[Date Index] [Thread Index] [Author Index]

Re: Finding the Number of Pythagorean Triples below a bound

• To: mathgroup at smc.vnet.net
• Subject: [mg68280] 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:49 -0400 (EDT)
• References: <200607290500.BAA04928@smc.vnet.net>
• Sender: owner-wri-mathgroup at wolfram.com

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

• References:
  • Finding the Number of Pythagorean Triples below a bound
    • From: titus_piezas@yahoo.com