Problem Statement

A Pythagorean triple is defined as three positive integers

( a , b , c )

{\displaystyle (a,b,c)}

where

a < b < c

{\displaystyle a<b<c}

, and

a

2

b

2

=

c

2

.

{\displaystyle a^{2}+b^{2}=c^{2}.}

They are called primitive triples if

a , b , c

{\displaystyle a,b,c}

are co-prime, that is, if their pairwise greatest common divisors

g c d

( a , b )

g c d

( a , c )

g c d

( b , c )

1

{\displaystyle {\rm {gcd}}(a,b)={\rm {gcd}}(a,c)={\rm {gcd}}(b,c)=1}

. Because of their relationship through the Pythagorean theorem, a, b, and c are co-prime if a and b are co-prime (

g c d

( a , b )

1

{\displaystyle {\rm {gcd}}(a,b)=1}

).   Each triple forms the length of the sides of a right triangle, whose perimeter is

P

a + b + c

{\displaystyle P=a+b+c}

.

The task is to determine how many Pythagorean triples there are with a perimeter no larger than 100 and the number of these that are primitive.

Deal with large values.   Can your program handle a maximum perimeter of 1,000,000?   What about 10,000,000?   100,000,000? Note: the extra credit is not for you to demonstrate how fast your language is compared to others;   you need a proper algorithm to solve them in a timely manner.

Let's start with the solution:

let isqrt n =
   let rec iter t =
      let d = n - t*t in
      if (0 <= d) && (d < t+t+1) (* t*t <= n < (t+1)*(t+1) *)
      then t else iter ((t+(n/t))/2)
   in iter 1

let rec gcd a b =
   let t = a mod b in
   if t = 0 then b else gcd b t

let coprime a b = gcd a b = 1

let num_to ms =
   let ctr = ref 0 in
   let prim_ctr = ref 0 in
   let max_m = isqrt (ms/2) in
   for m = 2 to max_m do
      for j = 0 to (m/2) - 1 do
         let n = m-(2*j+1) in
         if coprime m n then
            let s = 2*m*(m+n) in
            if s <= ms then
               (ctr := !ctr + (ms/s); incr prim_ctr)
      done
   done;
  (!ctr, !prim_ctr)

let show i =
  let s, p = num_to i in
  Printf.printf "For perimeters up to %d there are %d total and %d primitive\n%!" i s p;;

List.iter show [ 100; 1000; 10000; 100000; 1000000; 10000000; 100000000 ]