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}

.



Task

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.



Extra credit
 
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.



Related tasks

   Euler's sum of powers conjecture 
   List comprehensions
   Pythagorean quadruples 


 
show :-
    Data = [100, 1_000, 10_000, 100_000, 1_000_000, 10_000_000, 100_000_000],
    forall(
        member(Max, Data),
        (count_triples(Max, Total, Prim),
         format("upto ~D, there are ~D Pythagorean triples (~D primitive.)~n", [Max, Total, Prim]))).
 
div(A, B, C) :- C is A div B.
 
count_triples(Max, Total, Prims) :-
    findall(S, (triple(Max, A, B, C), S is A + B + C), Ps),
    length(Ps, Prims),
    maplist(div(Max), Ps, Counts), sumlist(Counts, Total).
 
% - between_by/4
 
between_by(A, B, N, K) :-
    C is (B - A) div N,
    between(0, C, J),
    K is N*J + A.
 
% - Pythagorean triple generator
 
triple(P, A, B, C) :-
    Max is floor(sqrt(P/2)) - 1,
    between(0, Max, M),
    Start is (M /\ 1) + 1, succ(Pm, M),
    between_by(Start, Pm, 2, N),
    gcd(M, N) =:= 1,
    X is M*M - N*N,
    Y is 2*M*N,
    C is M*M + N*N,
    order2(X, Y, A, B),
    (A + B + C) =< P.
 
order2(A, B, A, B) :- A < B, !.
order2(A, B, B, A).
 