# Pythagorean triples

Pythagorean triples
You are encouraged to solve this task according to the task description, using any language you may know.

A Pythagorean triple is defined as three positive integers ${\displaystyle (a,b,c)}$ where ${\displaystyle a, and ${\displaystyle a^{2}+b^{2}=c^{2}.}$

They are called primitive triples if ${\displaystyle a,b,c}$ are co-prime, that is, if their pairwise greatest common divisors ${\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 (${\displaystyle {\rm {gcd}}(a,b)=1}$).

Each triple forms the length of the sides of a right triangle, whose perimeter is ${\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.

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.

## 360 Assembly

Translation of: Perl
*        Pythagorean triples -     12/06/2018PYTHTRI  CSECT         USING  PYTHTRI,R13        base register         B      72(R15)            skip savearea         DC     17F'0'             savearea         SAVE   (14,12)            save previous context         ST     R13,4(R15)         link backward         ST     R15,8(R13)         link forward         LR     R13,R15            set addressability         MVC    PMAX,=F'1'         pmax=1         LA     R6,1               i=1       DO WHILE=(C,R6,LE,=F'6')    do i=1 to 6         L      R5,PMAX              pmax         MH     R5,=H'10'            *10         ST     R5,PMAX              pmax=pmax*10         MVC    PRIM,=F'0'           prim=0         MVC    COUNT,=F'0'          count=0              L      R1,PMAX              pmax         BAL    R14,ISQRT            isqrt(pmax)         SRA    R0,1                 /2         ST     R0,NMAX              nmax=isqrt(pmax)/2         LA     R7,1                 n=1       DO WHILE=(C,R7,LE,NMAX)       do n=1 to nmax         LA     R9,1(R7)               m=n+1         LR     R5,R9                  m         AR     R5,R7                  +n         MR     R4,R9                  *m         SLA    R5,1                   *2         LR     R8,R5                  p=2*m*(m+n)       DO WHILE=(C,R8,LE,PMAX)         do while p<=pmax          LR     R1,R9                    m         LR     R2,R7                    n         BAL    R14,GCD                  gcd(m,n)       IF C,R0,EQ,=F'1' THEN             if gcd(m,n)=1 then         L      R2,PRIM                    prim         LA     R2,1(R2)                   +1         ST     R2,PRIM                    prim=prim+1         L      R4,PMAX                    pmax         SRDA   R4,32                      ~         DR     R4,R8                      /p         A      R5,COUNT                   +count         ST     R5,COUNT                   count=count+pmax/p       ENDIF    ,                        endif         LA     R9,2(R9)                 m=m+2         LR     R5,R9                    m         AR     R5,R7                    +n         MR     R4,R9                    *m         SLA    R5,1                     *2         LR     R8,R5                    p=2*m*(m+n)       ENDDO    ,                      enddo n         LA     R7,1(R7)               n++       ENDDO    ,                    enddo n         L      R1,PMAX              pmax         XDECO  R1,XDEC              edit pmax         MVC    PG+15(9),XDEC+3      output pmax         L      R1,COUNT             count         XDECO  R1,XDEC              edit count         MVC    PG+33(9),XDEC+3      output count         L      R1,PRIM              prim         XDECO  R1,XDEC              edit prim         MVC    PG+55(9),XDEC+3      output prim         XPRNT  PG,L'PG              print         LA     R6,1(R6)             i++       ENDDO    ,                  enddo i         L      R13,4(0,R13)       restore previous savearea pointer         RETURN (14,12),RC=0       restore registers from calling savNMAX     DS     F                  nmaxPMAX     DS     F                  pmaxCOUNT    DS     F                  countPRIM     DS     F                  primPG     DC CL80'Max Perimeter: ........., Total: ........., Primitive:'XDEC     DS     CL12GCD      EQU    *  --------------- function gcd(a,b)         STM    R2,R7,GCDSA        save context         LR     R3,R1              c=a         LR     R4,R2              d=b      GCDLOOP  LR     R6,R3              c         SRDA   R6,32              ~         DR     R6,R4              /d         LTR    R6,R6              if c mod d=0         BZ     GCDELOOP           then leave loop         LR     R5,R6              e=c mod d         LR     R3,R4              c=d         LR     R4,R5              d=e         B      GCDLOOP            loopGCDELOOP LR     R0,R4              return(d)         LM     R2,R7,GCDSA        restore context         BR     R14                returnGCDSA    DS     6A                 context storeISQRT    EQU    *  --------------- function isqrt(n)                            STM    R3,R10,ISQRTSA     save context         LR     R6,R1              n=r1         LR     R10,R6             sqrtn=n         SRA    R10,1              sqrtn=n/2       IF LTR,R10,Z,R10 THEN       if sqrtn=0 then          LA     R10,1                sqrtn=1       ELSE     ,                  else          LA     R9,0                 snm2=0         LA     R8,0                 snm1=0         LA     R7,0                 sn=0         LA     R3,0                 okexit=0       DO UNTIL=(C,R3,EQ,=A(1))      do until okexit=1          AR     R10,R7                 sqrtn=sqrtn+sn         LR     R9,R8                  snm2=snm1         LR     R8,R7                  snm1=sn         LR     R4,R6                  n         SRDA   R4,32                  ~         DR     R4,R10                 /sqrtn         SR     R5,R10                 -sqrtn         SRA    R5,1                   /2         LR     R7,R5                  sn=(n/sqrtn-sqrtn)/2       IF C,R7,EQ,=F'0',OR,CR,R7,EQ,R9 THEN  if sn=0 or sn=snm2 then         LA     R3,1                     okexit=1       ENDIF    ,                      endif       ENDDO    ,                    enddo until       ENDIF    ,                  endif         LR     R5,R10             sqrtn         MR     R4,R10             *sqrtn       IF CR,R5,GT,R6 THEN         if sqrtn*sqrtn>n then         BCTR   R10,0                sqrtn=sqrtn-1       ENDIF    ,                  endif         LR     R0,R10             return(sqrtn)         LM     R3,R10,ISQRTSA     restore context         BR     R14                returnISQRTSA  DS     8A                 context store         YREGS         END    PYTHTRI
Output:
Max Perimeter:        10, Total:         0, Primitive:         0
Max Perimeter:       100, Total:        17, Primitive:         7
Max Perimeter:      1000, Total:       325, Primitive:        70
Max Perimeter:     10000, Total:      4858, Primitive:       703
Max Perimeter:    100000, Total:     64741, Primitive:      7026
Max Perimeter:   1000000, Total:    808950, Primitive:     70229


Translation of efficient method from C, see the WP article. Compiles on gnat/gcc.

with Ada.Text_IO; procedure Pythagorean_Triples is    type Large_Natural is range 0 .. 2**63-1;     -- this is the maximum for gnat    procedure New_Triangle(A, B, C: Large_Natural;                          Max_Perimeter: Large_Natural;                          Total_Cnt, Primitive_Cnt: in out Large_Natural) is      Perimeter: constant Large_Natural := A + B + C;   begin      if Perimeter <= Max_Perimeter then         Primitive_Cnt := Primitive_Cnt + 1;         Total_Cnt     := Total_Cnt + Max_Perimeter / Perimeter;         New_Triangle(A-2*B+2*C,     2*A-B+2*C,    2*A-2*B+3*C,   Max_Perimeter, Total_Cnt, Primitive_Cnt);         New_Triangle(A+2*B+2*C,     2*A+B+2*C,    2*A+2*B+3*C,   Max_Perimeter, Total_Cnt, Primitive_Cnt);         New_Triangle(2*B+2*C-A,     B+2*C-2*A,    2*B+3*C-2*A,   Max_Perimeter, Total_Cnt, Primitive_Cnt);      end if;   end New_Triangle;    T_Cnt, P_Cnt: Large_Natural; begin   for I in 1 .. 9 loop      T_Cnt := 0;      P_Cnt := 0;      New_Triangle(3,4,5, 10**I, Total_Cnt => T_Cnt, Primitive_Cnt => P_Cnt);      Ada.Text_IO.Put_Line("Up to 10 **" & Integer'Image(I) & " :" &                             Large_Natural'Image(T_Cnt) & " Triples," &                             Large_Natural'Image(P_Cnt) & " Primitives");   end loop;end Pythagorean_Triples;

Output:

Up to 10 ** 1 : 0 Triples, 0 Primitives
Up to 10 ** 2 : 17 Triples, 7 Primitives
Up to 10 ** 3 : 325 Triples, 70 Primitives
Up to 10 ** 4 : 4858 Triples, 703 Primitives
Up to 10 ** 5 : 64741 Triples, 7026 Primitives
Up to 10 ** 6 : 808950 Triples, 70229 Primitives
Up to 10 ** 7 : 9706567 Triples, 702309 Primitives
Up to 10 ** 8 : 113236940 Triples, 7023027 Primitives
Up to 10 ** 9 : 1294080089 Triples, 70230484 Primitives

## ANSI Standard BASIC

Translation of BBC BASIC Program

100 DECLARE EXTERNAL SUB tri110 !120 PUBLIC NUMERIC U0(3,3), U1(3,3), U2(3,3), all, prim130 DIM seed(3)140 MAT READ U0, U1, U2150 DATA 1, -2, 2, 2, -1, 2, 2, -2, 3160 DATA 1, 2, 2, 2, 1, 2, 2, 2, 3170 DATA -1, 2, 2, -2, 1, 2, -2, 2, 3180 !190 MAT READ seed200 DATA 3, 4, 5210 FOR power  = 1 TO 7220    LET all  = 0230    LET prim  = 0240    CALL tri(seed, 10^power , all , prim)250    PRINT "Up to 10^";power,260    PRINT USING "######### triples ######### primitives":all,prim270 NEXT power280 END290 !300 EXTERNAL SUB tri(i(), mp, all, prim)310 DECLARE EXTERNAL FUNCTION SUM320 DECLARE NUMERIC t(3)330 !340 IF SUM(i) > mp THEN EXIT SUB350 LET prim = prim + 1360 LET all  = all + INT(mp  / SUM(i))370 !380 MAT t = U0 * i390 CALL tri(t, mp , all , prim)400 MAT t = U1 * i410 CALL tri(t, mp , all , prim)420 MAT t = U2 * i430 CALL tri(t, mp , all , prim)440 END SUB450 !460 EXTERNAL FUNCTION SUM(a())470 LET temp = 0480 FOR i=LBOUND(a) TO UBOUND(a)490    LET temp = temp + a(i)500 NEXT i510 LET SUM = temp520 END FUNCTION

## AutoHotkey

#NoEnvSetBatchLines, -1#SingleInstance, Force ; Greatest common divisor, from http://rosettacode.org/wiki/Greatest_common_divisor#AutoHotkeygcd(a,b) {	Return b=0 ? Abs(a) : Gcd(b,mod(a,b))} count_triples(max) {	primitives := 0, triples := 0, m := 2	while m <= (max / 2)**0.5	{		n := mod(m, 2) + 1		,p := 2*m*(m + n)		, delta := 4*m		while n < m and p <= max			gcd(m, n) = 1				? (primitives++				, triples += max // p)				: ""			, n += 2			, p += delta		m++	}	Return primitives " primitives out of " triples " triples"} Loop, 8	Msgbox % 10**A_Index ": " count_triples(10**A_Index)
Output:
10: 0 primitives out of 0 triples
100: 7 primitives out of 17 triples
1000: 70 primitives out of 325 triples
10000: 703 primitives out of 4858 triples
100000: 7026 primitives out of 64741 triples
1000000: 70229 primitives out of 808950 triples
10000000: 702309 primitives out of 9706567 triples
100000000: 7023027 primitives out of 113236940 triples

## AWK

 # syntax: GAWK -f PYTHAGOREAN_TRIPLES.AWK# converted from GoBEGIN {    printf("%5s %11s %11s %11s %s\n","limit","limit","triples","primitives","seconds")    for (max_peri=10; max_peri<=1E9; max_peri*=10) {      t = systime()      prim = 0      total = 0      new_tri(3,4,5)      printf("10^%-2d %11d %11d %11d %d\n",++n,max_peri,total,prim,systime()-t)    }    exit(0)}function new_tri(s0,s1,s2,  p) {    p = s0 + s1 + s2    if (p <= max_peri) {      prim++      total += int(max_peri / p)      new_tri(+1*s0-2*s1+2*s2,+2*s0-1*s1+2*s2,+2*s0-2*s1+3*s2)      new_tri(+1*s0+2*s1+2*s2,+2*s0+1*s1+2*s2,+2*s0+2*s1+3*s2)      new_tri(-1*s0+2*s1+2*s2,-2*s0+1*s1+2*s2,-2*s0+2*s1+3*s2)    }}
Output:
limit       limit     triples  primitives seconds
10^1           10           0           0 0
10^2          100          17           7 0
10^3         1000         325          70 0
10^4        10000        4858         703 0
10^5       100000       64741        7026 0
10^6      1000000      808950       70229 0
10^7     10000000     9706567      702309 2
10^8    100000000   113236940     7023027 12
10^9   1000000000  1294080089    70230484 116


## BBC BASIC

The built-in array arithmetic is very well suited to this task!

      DIM U0%(2,2), U1%(2,2), U2%(2,2), seed%(2)      U0%() =  1, -2, 2,  2, -1, 2,  2, -2, 3      U1%() =  1,  2, 2,  2,  1, 2,  2,  2, 3      U2%() = -1,  2, 2, -2,  1, 2, -2,  2, 3       seed%() = 3, 4, 5      FOR power% = 1 TO 7        all% = 0 : prim% = 0        PROCtri(seed%(), 10^power%, all%, prim%)        PRINT "Up to 10^"; power%, ": " all% " triples" prim% " primitives"      NEXT      END       DEF PROCtri(i%(), mp%, RETURN all%, RETURN prim%)      LOCAL t%() : DIM t%(2)       IF SUM(i%()) > mp% ENDPROC      prim% += 1      all% += mp% DIV SUM(i%())       t%() = U0%() . i%()      PROCtri(t%(), mp%, all%, prim%)      t%() = U1%() . i%()      PROCtri(t%(), mp%, all%, prim%)      t%() = U2%() . i%()      PROCtri(t%(), mp%, all%, prim%)      ENDPROC

Output:

Up to 10^1:          0 triples         0 primitives
Up to 10^2:         17 triples         7 primitives
Up to 10^3:        325 triples        70 primitives
Up to 10^4:       4858 triples       703 primitives
Up to 10^5:      64741 triples      7026 primitives
Up to 10^6:     808950 triples     70229 primitives
Up to 10^7:    9706567 triples    702309 primitives
Up to 10^8:  113236940 triples   7023027 primitives


## Bracmat

Translation of: C
(pythagoreanTriples=  total prim max-peri U.       (.(1,-2,2) (2,-1,2) (2,-2,3))        (.(1,2,2) (2,1,2) (2,2,3))        (.(-1,2,2) (-2,1,2) (-2,2,3))    : ?U  & ( new-tri    =     i t p Urows Urow Ucols        , a b c loop A B C      .     !arg:(,?a,?b,?c)          & !a+!b+!c:~>!max-peri:?p          & 1+!prim:?prim          & div$(!max-peri.!p)+!total:?total & !U:?Urows & ( loop = !Urows:(.?Urow) ?Urows & !Urow:?Ucols & :?t & whl ' ( !Ucols:(?A,?B,?C) ?Ucols & (!t,!a*!A+!b*!B+!c*!C):?t ) & new-tri$!t              & !loop            )          & !loop        |    )  & ( Main    =   seed      .   (,3,4,5):?seed        & 10:?max-peri        &   whl          ' ( 0:?total:?prim            & new-tri$!seed & out$ ( str                $( "Up to " !max-peri ": " !total " triples, " !prim " primitives." ) ) & !max-peri*10:~>10000000:?max-peri ) ) & Main$); pythagoreanTriples$;  Output (under Linux): Up to 10: 0 triples, 0 primitives. Up to 100: 17 triples, 7 primitives. Up to 1000: 325 triples, 70 primitives. Up to 10000: 4858 triples, 703 primitives. Up to 100000: 64741 triples, 7026 primitives. Up to 1000000: 808950 triples, 70229 primitives. Up to 10000000: 9706567 triples, 702309 primitives. Up to 100000000: 113236940 triples, 7023027 primitives. Under Windows XP Command prompt the last result is unattainable due to stack overflow. With very few changes we can get rid of the stack exhausting recursion. Instead of calling new-tri recursively, be push the triples to test onto a stack and return to the Main function. In the innermost loop we pop a triple from the stack and call new-tri. The memory overhead is only a few megabytes for a max perimeter of 100,000,000. On my Windows XP box the whole computation takes at least 15 minutes! Given enough time (and memory), the program can compute results for larger perimeters. (pythagoreanTriples= total prim max-peri U stack. (.(1,-2,2) (2,-1,2) (2,-2,3)) (.(1,2,2) (2,1,2) (2,2,3)) (.(-1,2,2) (-2,1,2) (-2,2,3)) : ?U & ( new-tri = i t p Urows Urow Ucols Ucol , a b c loop A B C . !arg:(,?a,?b,?c) & !a+!b+!c:~>!max-peri:?p & 1+!prim:?prim & div$(!max-peri.!p)+!total:?total          & !U:?Urows          & ( loop            =   !Urows:(.?Urow) ?Urows              & !Urow:?Ucols              & :?t              &   whl                ' ( !Ucols:(?A,?B,?C) ?Ucols                  & (!t,!a*!A+!b*!B+!c*!C):?t                  )              & !t !stack:?stack              & !loop            )          & !loop        |    )  & ( Main    =   seed      .   10:?max-peri        &   whl          ' ( 0:?total:?prim            & (,3,4,5):?stack            &   whl              ' (!stack:%?seed ?stack&new-tri$!seed) & out$ ( str                $( "Up to " !max-peri ": " !total " triples, " !prim " primitives." ) ) & !max-peri*10:~>100000000:?max-peri ) ) & Main$); pythagoreanTriples$; ## C Sample implemention; naive method, patentedly won't scale to larger numbers, despite the attempt to optimize it. Calculating up to 10000 is already a test of patience. #include <stdio.h>#include <stdlib.h> typedef unsigned long long xint;typedef unsigned long ulong; inline ulong gcd(ulong m, ulong n){ ulong t; while (n) { t = n; n = m % n; m = t; } return m;} int main(){ ulong a, b, c, pytha = 0, prim = 0, max_p = 100; xint aa, bb, cc; for (a = 1; a <= max_p / 3; a++) { aa = (xint)a * a; printf("a = %lu\r", a); /* show that we are working */ fflush(stdout); /* max_p/2: valid limit, because one side of triangle * must be less than the sum of the other two */ for (b = a + 1; b < max_p/2; b++) { bb = (xint)b * b; for (c = b + 1; c < max_p/2; c++) { cc = (xint)c * c; if (aa + bb < cc) break; if (a + b + c > max_p) break; if (aa + bb == cc) { pytha++; if (gcd(a, b) == 1) prim++; } } } } printf("Up to %lu, there are %lu triples, of which %lu are primitive\n", max_p, pytha, prim); return 0;} output: Up to 100, there are 17 triples, of which 7 are primitive Efficient method, generating primitive triples only as described in the same WP article: #include <stdio.h>#include <stdlib.h>#include <stdint.h> /* should be 64-bit integers if going over 1 billion */typedef unsigned long xint;#define FMT "%lu" xint total, prim, max_peri;xint U[][9] = {{ 1, -2, 2, 2, -1, 2, 2, -2, 3}, { 1, 2, 2, 2, 1, 2, 2, 2, 3}, {-1, 2, 2, -2, 1, 2, -2, 2, 3}}; void new_tri(xint in[]){ int i; xint t[3], p = in[0] + in[1] + in[2]; if (p > max_peri) return; prim ++; /* for every primitive triangle, its multiples would be right-angled too; * count them up to the max perimeter */ total += max_peri / p; /* recursively produce next tier by multiplying the matrices */ for (i = 0; i < 3; i++) { t[0] = U[i][0] * in[0] + U[i][1] * in[1] + U[i][2] * in[2]; t[1] = U[i][3] * in[0] + U[i][4] * in[1] + U[i][5] * in[2]; t[2] = U[i][6] * in[0] + U[i][7] * in[1] + U[i][8] * in[2]; new_tri(t); }} int main(){ xint seed[3] = {3, 4, 5}; for (max_peri = 10; max_peri <= 100000000; max_peri *= 10) { total = prim = 0; new_tri(seed); printf( "Up to "FMT": "FMT" triples, "FMT" primitives.\n", max_peri, total, prim); } return 0;} Output Up to 10: 0 triples, 0 primitives.Up to 100: 17 triples, 7 primitives.Up to 1000: 325 triples, 70 primitives.Up to 10000: 4858 triples, 703 primitives.Up to 100000: 64741 triples, 7026 primitives.Up to 1000000: 808950 triples, 70229 primitives.Up to 10000000: 9706567 triples, 702309 primitives.Up to 100000000: 113236940 triples, 7023027 primitives. Same as above, but with loop unwound and third recursion eliminated: #include <stdio.h>#include <stdlib.h>#include <stdint.h> /* should be 64-bit integers if going over 1 billion */typedef unsigned long xint;#define FMT "%lu" xint total, prim, max_peri; void new_tri(xint in[]){ int i; xint t[3], p; xint x = in[0], y = in[1], z = in[2]; recur: p = x + y + z; if (p > max_peri) return; prim ++; total += max_peri / p; t[0] = x - 2 * y + 2 * z; t[1] = 2 * x - y + 2 * z; t[2] = t[1] - y + z; new_tri(t); t[0] += 4 * y; t[1] += 2 * y; t[2] += 4 * y; new_tri(t); z = t[2] - 4 * x; y = t[1] - 4 * x; x = t[0] - 2 * x; goto recur;} int main(){ xint seed[3] = {3, 4, 5}; for (max_peri = 10; max_peri <= 100000000; max_peri *= 10) { total = prim = 0; new_tri(seed); printf( "Up to "FMT": "FMT" triples, "FMT" primitives.\n", max_peri, total, prim); } return 0;} ## C# Based on Ada example, which is a translation of efficient method from C, see the WP article. using System; namespace RosettaCode.CSharp{ class Program { static void Count_New_Triangle(ulong A, ulong B, ulong C, ulong Max_Perimeter, ref ulong Total_Cnt, ref ulong Primitive_Cnt) { ulong Perimeter = A + B + C; if (Perimeter <= Max_Perimeter) { Primitive_Cnt = Primitive_Cnt + 1; Total_Cnt = Total_Cnt + Max_Perimeter / Perimeter; Count_New_Triangle(A + 2 * C - 2 * B, 2 * A + 2 * C - B, 2 * A + 3 * C - 2 * B, Max_Perimeter, ref Total_Cnt, ref Primitive_Cnt); Count_New_Triangle(A + 2 * B + 2 * C, 2 * A + B + 2 * C, 2 * A + 2 * B + 3 * C, Max_Perimeter, ref Total_Cnt, ref Primitive_Cnt); Count_New_Triangle(2 * B + 2 * C - A, B + 2 * C - 2 * A, 2 * B + 3 * C - 2 * A, Max_Perimeter, ref Total_Cnt, ref Primitive_Cnt); } } static void Count_Pythagorean_Triples() { ulong T_Cnt, P_Cnt; for (int I = 1; I <= 8; I++) { T_Cnt = 0; P_Cnt = 0; ulong ExponentNumberValue = (ulong)Math.Pow(10, I); Count_New_Triangle(3, 4, 5, ExponentNumberValue, ref T_Cnt, ref P_Cnt); Console.WriteLine("Perimeter up to 10E" + I + " : " + T_Cnt + " Triples, " + P_Cnt + " Primitives"); } } static void Main(string[] args) { Count_Pythagorean_Triples(); } }} Output: Perimeter up to 10E1 : 0 Triples, 0 Primitives Perimeter up to 10E2 : 17 Triples, 7 Primitives Perimeter up to 10E3 : 325 Triples, 70 Primitives Perimeter up to 10E4 : 4858 Triples, 703 Primitives Perimeter up to 10E5 : 64741 Triples, 7026 Primitives Perimeter up to 10E6 : 808950 Triples, 70229 Primitives Perimeter up to 10E7 : 9706567 Triples, 702309 Primitives Perimeter up to 10E8 : 113236940 Triples, 7023027 Primitives ## Clojure This version is based on Euclid's formula: for each pair (m,n) such that m>n>0, m and n coprime and of opposite polarity (even/odd), there is a primitive Pythagorean triple. It can be proven that the converse is true as well. (defn gcd [a b] (if (zero? b) a (recur b (mod a b)))) (defn pyth [peri] (for [m (range 2 (Math/sqrt (/ peri 2))) n (range (inc (mod m 2)) m 2) ; n<m, opposite polarity :let [p (* 2 m (+ m n))] ; = a+b+c for this (m,n) :while (<= p peri) :when (= 1 (gcd m n)) :let [m2 (* m m), n2 (* n n), [a b] (sort [(- m2 n2) (* 2 m n)]), c (+ m2 n2)] k (range 1 (inc (quot peri p)))] [(= k 1) (* k a) (* k b) (* k c)])) (defn rcount [ts] ; (->> peri pyth rcount) produces [total, primitive] counts (reduce (fn [[total prims] t] [(inc total), (if (first t) (inc prims) prims)]) [0 0] ts)) To handle really large perimeters, we can dispense with actually generating the triples and just calculate the counts: (defn pyth-count [peri] (reduce (fn [[total prims] k] [(+ total k), (inc prims)]) [0 0] (for [m (range 2 (Math/sqrt (/ peri 2))) n (range (inc (mod m 2)) m 2) ; n<m, opposite polarity :let [p (* 2 m (+ m n))] ; = a+b+c for this (m,n) :while (<= p peri) :when (= 1 (gcd m n))] (quot peri p)))) ## CoffeeScript This algorithm scales linearly with the max perimeter. It uses two loops that are capped by the square root of the half-perimeter to examine/count provisional values of m and n, where m and n generate a, b, c, and p using simple number theory.  gcd = (x, y) -> return x if y == 0 gcd(y, x % y) # m,n generate primitive Pythag triples## preconditions:# m, n are integers of different parity# m > n# gcd(m,n) == 1 (coprime)## m, n generate: [m*m - n*n, 2*m*n, m*m + n*n]# perimeter is 2*m*m + 2*m*n = 2 * m * (m+n)count_triples = (max_perim) -> num_primitives = 0 num_triples = 0 m = 2 upper_limit = Math.sqrt max_perim / 2 while m <= upper_limit n = m % 2 + 1 p = 2*m*m + 2*m*n delta = 4*m while n < m and p <= max_perim if gcd(m, n) == 1 num_primitives += 1 num_triples += Math.floor max_perim / p n += 2 p += delta m += 1 console.log num_primitives, num_triples max_perim = Math.pow 10, 9 # takes under a minutecount_triples(max_perim)  output time coffee pythag_triples.coffee 70230484 1294080089 real 0m45.989s  ## Common Lisp (defun mmul (a b) (loop for x in a collect (loop for y in x for z in b sum (* y z)))) (defun count-tri (lim &aux (prim 0) (cnt 0)) (labels ((count1 (tr &aux (peri (reduce #'+ tr))) (when (<= peri lim) (incf prim) (incf cnt (truncate lim peri)) (count1 (mmul '(( 1 -2 2) ( 2 -1 2) ( 2 -2 3)) tr)) (count1 (mmul '(( 1 2 2) ( 2 1 2) ( 2 2 3)) tr)) (count1 (mmul '((-1 2 2) (-2 1 2) (-2 2 3)) tr))))) (count1 '(3 4 5)) (format t "~a: ~a prim, ~a all~%" lim prim cnt))) (loop for p from 2 do (count-tri (expt 10 p))) output 100: 7 prim, 17 all1000: 70 prim, 325 all10000: 703 prim, 4858 all100000: 7026 prim, 64741 all1000000: 70229 prim, 808950 all10000000: 702309 prim, 9706567 all... ## Crystal Translation of: Ruby class PythagoranTriplesCounter def initialize(limit = 0) @limit = limit @total = 0 @primitives = 0 generate_triples(3, 4, 5) end def total; @total end def primitives; @primitives end private def generate_triples(a, b, c) perim = a + b + c return if perim > @limit @primitives += 1 @total += @limit / perim generate_triples( a-2*b+2*c, 2*a-b+2*c, 2*a-2*b+3*c ) generate_triples( a+2*b+2*c, 2*a+b+2*c, 2*a+2*b+3*c ) generate_triples(-a+2*b+2*c,-2*a+b+2*c,-2*a+2*b+3*c ) endend perim = 10while perim <= 100_000_000 c = PythagoranTriplesCounter.new perim p [perim, c.total, c.primitives] perim *= 10end output [10, 0, 0] [100, 17, 7] [1000, 325, 70] [10000, 4858, 703] [100000, 64741, 7026] [1000000, 808950, 70229] [10000000, 9706567, 702309] [100000000, 113236940, 7023027] ## D ### Lazy Functional Version With hints from the Haskell solution. void main() @safe { import std.stdio, std.range, std.algorithm, std.typecons, std.numeric; enum triples = (in uint n) pure nothrow @safe /*@nogc*/ => iota(1, n + 1) .map!(z => iota(1, z + 1) .map!(x => iota(x, z + 1).map!(y => tuple(x, y, z)))) .joiner.joiner .filter!(t => t[0] ^^ 2 + t[1] ^^ 2 == t[2] ^^ 2 && t[].only.sum <= n) .map!(t => tuple(t[0 .. 2].gcd == 1, t[])); auto xs = triples(100); writeln("Up to 100 there are ", xs.count, " triples, ", xs.filter!q{ a[0] }.count, " are primitive.");} Output: Up to 100 there are 17 triples, 7 are primitive. ### Shorter Version ulong[2] tri(ulong lim, ulong a=3, ulong b=4, ulong c=5)pure nothrow @safe @nogc { immutable l = a + b + c; if (l > lim) return [0, 0]; typeof(return) r = [1, lim / l]; r[] += tri(lim, a - 2*b + 2*c, 2*a - b + 2*c, 2*a - 2*b + 3*c)[]; r[] += tri(lim, a + 2*b + 2*c, 2*a + b + 2*c, 2*a + 2*b + 3*c)[]; r[] += tri(lim, -a + 2*b + 2*c, -2*a + b + 2*c, -2*a + 2*b + 3*c)[]; return r;} void main() /*@safe*/ { import std.stdio; foreach (immutable p; 1 .. 9) writeln(10 ^^ p, ' ', tri(10 ^^ p));} Output: 10 [0, 0] 100 [7, 17] 1000 [70, 325] 10000 [703, 4858] 100000 [7026, 64741] 1000000 [70229, 808950] 10000000 [702309, 9706567] 100000000 [7023027, 113236940] Run-time (32 bit system): about 0.80 seconds with ldc2. ### Short SIMD Version With LDC compiler this is a little faster than the precedent version (remove @nogc to compile it with the current version of LDC compiler). import std.stdio, core.simd; ulong2 tri(in ulong lim, in ulong a=3, in ulong b=4, in ulong c=5)pure nothrow @safe @nogc { immutable l = a + b + c; if (l > lim) return [0, 0]; typeof(return) r = [1, lim / l]; r += tri(lim, a - 2*b + 2*c, 2*a - b + 2*c, 2*a - 2*b + 3*c); r += tri(lim, a + 2*b + 2*c, 2*a + b + 2*c, 2*a + 2*b + 3*c); r += tri(lim, -a + 2*b + 2*c, -2*a + b + 2*c, -2*a + 2*b + 3*c); return r;} void main() /*@safe*/ { foreach (immutable p; 1 .. 9) writeln(10 ^^ p, ' ', tri(10 ^^ p).array);} The output is the same. Run-time (32 bit system): about 0.67 seconds with ldc2. ### Faster Version Translation of: C import std.stdio; alias Xuint = uint; // ulong if going over 1 billion. __gshared Xuint nTriples, nPrimitives, limit; void countTriples(Xuint x, Xuint y, Xuint z) nothrow @nogc { while (true) { immutable p = x + y + z; if (p > limit) return; nPrimitives++; nTriples += limit / p; auto t0 = x - 2 * y + 2 * z; auto t1 = 2 * x - y + 2 * z; auto t2 = t1 - y + z; countTriples(t0, t1, t2); t0 += 4 * y; t1 += 2 * y; t2 += 4 * y; countTriples(t0, t1, t2); z = t2 - 4 * x; y = t1 - 4 * x; x = t0 - 2 * x; }} void main() { foreach (immutable p; 1 .. 9) { limit = Xuint(10) ^^ p; nTriples = nPrimitives = 0; countTriples(3, 4, 5); writefln("Up to %11d: %11d triples, %9d primitives.", limit, nTriples, nPrimitives); }} Output: Up to 10: 0 triples, 0 primitives. Up to 100: 17 triples, 7 primitives. Up to 1000: 325 triples, 70 primitives. Up to 10000: 4858 triples, 703 primitives. Up to 100000: 64741 triples, 7026 primitives. Up to 1000000: 808950 triples, 70229 primitives. Up to 10000000: 9706567 triples, 702309 primitives. Up to 100000000: 113236940 triples, 7023027 primitives. Run-time: about 0.27 seconds with ldc2. Using the power p up to 11, using ulong for xuint, and compiling with the dmd -L/STACK:10000000 switch to increase the stack size to about 10MB: Up to 10: 0 triples, 0 primitives. Up to 100: 17 triples, 7 primitives. Up to 1000: 325 triples, 70 primitives. Up to 10000: 4858 triples, 703 primitives. Up to 100000: 64741 triples, 7026 primitives. Up to 1000000: 808950 triples, 70229 primitives. Up to 10000000: 9706567 triples, 702309 primitives. Up to 100000000: 113236940 triples, 7023027 primitives. Up to 1000000000: 1294080089 triples, 70230484 primitives. Up to 10000000000: 14557915466 triples, 702304875 primitives. Total run-time up to 10_000_000_000: about 63 seconds. Waiting less than half an hour: Up to 10: 0 triples, 0 primitives. Up to 100: 17 triples, 7 primitives. Up to 1000: 325 triples, 70 primitives. Up to 10000: 4858 triples, 703 primitives. Up to 100000: 64741 triples, 7026 primitives. Up to 1000000: 808950 triples, 70229 primitives. Up to 10000000: 9706567 triples, 702309 primitives. Up to 100000000: 113236940 triples, 7023027 primitives. Up to 1000000000: 1294080089 triples, 70230484 primitives. Up to 10000000000: 14557915466 triples, 702304875 primitives. Up to 100000000000: 161750315680 triples, 7023049293 primitives. ## EDSAC order code Not much optimization is done in this code, which is best run on a fast simulator if the maximum perimeter is more than 10^5 or so. A maximum perimeter of 10^7 takes 340 million EDSAC orders (about 6 days on the original EDSAC). The number of primitive triples divided by the maximum perimeter seems to tend to a limit, which looks very much like (ln 2)/pi^2 = 0.07023049277 (see especially the FreeBASIC output below).  [Pythagorean triples for Rosetta code.Counts (1) all Pythagorean triples (2) primitive Pythagorean triples,with perimeter not greater than a given value. Library subroutine M3, Prints header and is then overwritten.Here, the last character sets the teleprinter to figures.] ..PZ [simulate blank tape] [email protected]@E8FEZPF @&*[email protected]&#. ..PZ [Library subroutine P7, prints long strictly positive integer;10 characters, right justified, padded left with spaces.Closed, even; 35 storage locations; working position 4D.] T 56 K [email protected]#@[email protected]@[email protected]@SFLDUFOFFFSFL4F [email protected]@XFT28#[email protected][email protected]@ [Subroutine for positive integer division. Input: 4D = dividend, 6D = divisor. Output: 4D = remainder, 6D = quotient. 37 locations; working locations 0D, 8D.] T 100 K [email protected]@[email protected]@[email protected]@ [email protected]@[email protected]@EFPF [Subroutine to return GCD of two non-negative 35-bit integers. Input: Integers at 4D, 6D. Output: GCD at 4D; changes 6D. 41 locations; working location 0D.] T 200 K [email protected]@[email protected]@[email protected]@[email protected]@T6D [email protected]@[email protected]@[email protected]@A6DT4DEFPF [************************ ROSETTA CODE TASK ************************* Subroutine to count Pythagorean triples with given maximum perimeter. Input: 0D = maximum perimeter. Output: 4D = number of triples, 6D = number of primitive. 0D is changed. Must be loaded at an even address. Uses the well-known fact that a primitive Pythagorean triple is of the form (m^2 - n^2, 2*m*n, m^2 + n^2) where m, n are coprime and of opposite parity.] T 300 K G K A 3 F [make link] E 16 @ [jump over variables and constants] [Double values are put here to ensure even address] [Variables] [2] P F P F [maximum perimeter] [4] P F P F [total number of Pythagorean triples] [6] P F P F [number of primitive Pythagorean triples] [8] P F P F [m] [10] P F P F [n] [Constants] [12] P D P F [double-value 1] [14] P1F P F [double-value 2] [Continue with code] [16] T 69 @ [plant link for return] A D [load maximum perimeter] T 2#@ [store locally] T 4#@ [initialize counts of triangles to 0] T 6#@ A 12#@ [load 1] T 8#@ [m := 1] [Next m, inc by 1] [23] T F [clear acc] A 8#@ [load m] A 12#@ [add 1] T 8#@ [update m] H 8#@ [mult reg := m] C 12#@ [acc := m AND 1] A 12#@ [add 1] T 10#@ [n := 1 if m even, 2 if m odd] [Here to count triangles arising from m, n. It's assumed m and n are known coprime.] [31] A 31 @ [call the count subroutine,] G 70 @ [result is in 6D] S 6 D [load negative count] G 40 @ [jump if count > 0] [No triangles found for this n. If n = 1 or 2 then whole thing is finished. Else move on to next m.] T F [clear acc] A 14#@ [load 2] S 10#@ [2 - n] G 23 @ [if n > 2, go to next m] E 64 @ [if n <= 2, exit] [Found triangles, count is in 6D] [40] T F [clear acc] A 4#@ [load total count] A 6 D [add count just found] T 4#@ [update total count] A 6#@ [load primitive count] A 12#@ [add 1] T 6#@ [update primitive count] [47] T F [clear acc] A 10#@ [load n] A 14#@ [add 2] U 10#@ [update n] S 8#@ [is n > m?] E 23 @ [if so, loop back for next m] [Test whether m and n are coprime.] T F [clear acc] A 8#@ [load m] T 4 D [to 4D for GCD routine] A 10#@ [load n] T 6 D [to 6D for GCD routine] A 58 @ [call GCD routine,] G 200 F [GCD is returned in 4D] A 4 D [load GCD] S 14#@ [is GCD = 1? (test by subtracting 2)] E 47 @ [no, go straight to next n] G 31 @ [yes, count triangles, then next n] [64] T F [exit, clear acc] A 4#@ [load total number of triples] T 4 D [return in 4D] A 6#@ [load number of primitive triples] T 6 D [return in 6D] [69] E F [2nd-level subroutine to count triangles arising from m, n. Assumes m, n are coprime and of opposite parity, and m is in the multiplier register. Result is returned in 6D.] [70] A 3 F [make and plant link for return] T 91 @ A 2#@ [acc := maximum perimeter] T 4 D [to 4D for division routine] A 8#@ [load m] A 10#@ [add n] T D [m + n to 0D] V D [acc := m*(m + n)] [Need to shift product 34 left to restore integer scaling. Since we want 2*m*(m+n), shift 35 left.] L F [13 left (maximum possible)] L F [13 more] L 128 F [9 more] T 6 D [perimeter to 6D for division routine] A 4 D [load maximum perimeter] S 6 D [is perimeter > maximum?] G 89 @ [quick exit if so] T F [clear acc] A 86 @ [call division routine,] G 100 F [leaves count in 6D] E 91 @ [jump to exit] [89] T F [acc := 0] T 6 D [return count = 0] [91] E F [Main routine. Load at an even address.] T 500 K G K [The initial maximum perimeter is repeatedly multiplied by 10] [0] P50F PF [initial maximum perimeter <---------- EDIT HERE] [2] P 3 F [number of values to calculate <---------- EDIT HERE] [3] P D [1] [4] P F P F [maximum perimeter] [6] P F P F [total number of triples] [8] P F P F [number of primitive triples] [10] P F [negative count of values] [11] # F [figures shift] [12] @ F [carriage return] [13] & F [line feed] [14] K 4096 F [null char] [Enter with acc = 0] [15] S 2 @ [initialize a negative counter] T 10 @ [(standard EDSAC practice)] A #@ [initialize maximum perimeter] T 4#@ [19] T F [clear acc] A 4#@ [load maximum perimeter] T D [to 0D for subroutine] A 22 @ [call subroutine to count triples] G 300 F A 4 D [returns total number in 4D] T 6#@ [save locally] A 6 D [returns number of primitive in 6D] T 8#@ [save locally] [Print the result] A 4#@ [load maximum perimeter] T D [to 0D for print subroutine] A 30 @ [call print subroutine] G 56 F A 6#@ [repeat for total number of triples] T D A 34 @ G 56 F A 8#@ [repeat for number of primitive triples] T D A 38 @ G 56 F O 12 @ O 13 @ A 10 @ [load negative count] A 3 @ [add 1] E 53 @ [out if reached 0] T 10 @ [else update count] A 4#@ [load max perimeter] U D [temp store] L 1 F [times 4] A D [times 5] L D [times 10] T 4#@ [update] E 19 @ [loop back] [53] O 14 @ [done; print null to flush printer buffer] Z F [stop] E 15 Z [define entry point] P F [acc = 0 on entry]  Output:  MAX PERIM TOTAL PRIM 100 17 7 1000 325 70 10000 4858 703 100000 64741 7026 1000000 808950 70229 10000000 9706567 702309  ## Eiffel  class APPLICATION create make feature make local perimeter: INTEGER do perimeter := 100 from until perimeter > 1000000 loop total := 0 primitive_triples := 0 count_pythagorean_triples (3, 4, 5, perimeter) io.put_string ("There are " + total.out + " triples, below " + perimeter.out + ". Of which " + primitive_triples.out + " are primitives.%N") perimeter := perimeter * 10 end end count_pythagorean_triples (a, b, c, perimeter: INTEGER) -- Total count of pythagorean triples and total count of primitve triples below perimeter. local p: INTEGER do p := a + b + c if p <= perimeter then primitive_triples := primitive_triples + 1 total := total + perimeter // p count_pythagorean_triples (a + 2 * (- b + c), 2 * (a + c) - b, 2 * (a - b + c) + c, perimeter) count_pythagorean_triples (a + 2 * (b + c), 2 * (a + c) + b, 2 * (a + b + c) + c, perimeter) count_pythagorean_triples (- a + 2 * (b + c), 2 * (- a + c) + b, 2 * (- a + b + c) + c, perimeter) end end feature {NONE} primitive_triples: INTEGER total: INTEGER end  Output: There are 17 triples, below 100. Of which 7 are primitives. There are 325 triples, below 1000. Of which 70 are primitives. There are 4858 triples, below 10000. Of which 703 are primitives. There are 64741 triples, below 100000. Of which 7026 are primitives. There are 808950 triples, below 1000000. Of which 70229 are primitives.  ## Elixir Translation of: Ruby defmodule RC do def count_triples(limit), do: count_triples(limit,3,4,5) defp count_triples(limit, a, b, c) when limit<(a+b+c), do: {0,0} defp count_triples(limit, a, b, c) do {p1, t1} = count_triples(limit, a-2*b+2*c, 2*a-b+2*c, 2*a-2*b+3*c) {p2, t2} = count_triples(limit, a+2*b+2*c, 2*a+b+2*c, 2*a+2*b+3*c) {p3, t3} = count_triples(limit,-a+2*b+2*c,-2*a+b+2*c,-2*a+2*b+3*c) {1+p1+p2+p3, div(limit, a+b+c)+t1+t2+t3} endend list = for n <- 1..8, do: Enum.reduce(1..n, 1, fn(_,acc)->10*acc end)Enum.each(list, fn n -> IO.inspect {n, RC.count_triples(n)} end) Output: {10, {0, 0}} {100, {7, 17}} {1000, {70, 325}} {10000, {703, 4858}} {100000, {7026, 64741}} {1000000, {70229, 808950}} {10000000, {702309, 9706567}} {100000000, {7023027, 113236940}}  ## Erlang %%%% Pythagorian triples in Erlang, J.W. Luiten%%-module(triples).-export([main/1]). %% Transformations t1, t2 and t3 to generate new triples t1(A, B, C) -> {A-2*B+2*C, 2*A-B+2*C, 2*A-2*B+3*C}.t2(A, B, C) -> {A+2*B+2*C, 2*A+B+2*C, 2*A+2*B+3*C}.t3(A, B, C) -> {2*B+2*C-A, B+2*C-2*A, 2*B+3*C-2*A}. %% Generation of triplescount_triples(A, B, C, Tot_acc, Cnt_acc, Max_perimeter) when (A+B+C) =< Max_perimeter -> Tot1 = Tot_acc + Max_perimeter div (A+B+C), {A1, B1, C1} = t1(A, B, C), {Tot2, Cnt2} = count_triples(A1, B1, C1, Tot1, Cnt_acc+1, Max_perimeter), {A2, B2, C2} = t2(A, B, C), {Tot3, Cnt3} = count_triples(A2, B2, C2, Tot2, Cnt2, Max_perimeter), {A3, B3, C3} = t3(A, B, C), {Tot4, Cnt4} = count_triples(A3, B3, C3, Tot3, Cnt3, Max_perimeter), {Tot4, Cnt4};count_triples(_A, _B, _C, Tot_acc, Cnt_acc, _Max_perimeter) -> {Tot_acc, Cnt_acc}. count_triples(A, B, C, Pow) -> Max = trunc(math:pow(10, Pow)), {Tot, Prim} = count_triples(A, B, C, 0, 0, Max), {Pow, Tot, Prim}. count_triples(Pow) -> count_triples(3, 4, 5, Pow). %% Display a single result.display_result({Pow, Tot, Prim}) -> io:format("Up to 10 ** ~w : ~w triples, ~w primitives~n", [Pow, Tot, Prim]). main(Max) -> L = lists:seq(1, Max), Answer = lists:map(fun(X) -> count_triples(X) end, L), lists:foreach(fun(Result) -> display_result(Result) end, Answer). Output: Up to 10 ** 1 : 0 triples, 0 primitives Up to 10 ** 2 : 17 triples, 7 primitives Up to 10 ** 3 : 325 triples, 70 primitives Up to 10 ** 4 : 4858 triples, 703 primitives Up to 10 ** 5 : 64741 triples, 7026 primitives Up to 10 ** 6 : 808950 triples, 70229 primitives Up to 10 ** 7 : 9706567 triples, 702309 primitives Up to 10 ** 8 : 113236940 triples, 7023027 primitives Up to 10 ** 9 : 1294080089 triples, 70230484 primitives Up to 10 ** 10 : 14557915466 triples, 702304875 primitives Up to 10 ** 11 : 161750315680 triples, 7023049293 primitives  ## ERRE PROGRAM PIT BEGIN PRINT(CHR$(12);) !CLS  PRINT(TIME$) FOR POWER=1 TO 7 DO PLIMIT=10#^POWER UPPERBOUND=INT(1+PLIMIT^0.5) PRIMITIVES=0 TRIPLES=0 EXTRAS=0 ! will count the in-range multiples of any primitive FOR M=2 TO UPPERBOUND DO FOR N=1+(M MOD 2=1) TO M-1 STEP 2 DO TERM1=2*M*N TERM2=M*M-N*N TERM3=M*M+N*N PERIMETER=TERM1+TERM2+TERM3 IF PERIMETER<=PLIMIT THEN TRIPLES=TRIPLES+1 A=TERM1 B=TERM2 REPEAT R=A-B*INT(A/B) A=B B=R UNTIL R<=0 ! we've found a primitive triple if a = 1, since hcf =1. ! and it is inside perimeter range. Save it in an array IF (A=1) AND (PERIMETER<=PLIMIT) THEN PRIMITIVES=PRIMITIVES+1 !----------------------------------------------- !swap so in increasing order of side length !----------------------------------------------- IF TERM1>TERM2 THEN SWAP(TERM1,TERM2) !----------------------------------------------- !we have the primitive & removed any multiples. !Now calculate ALL the multiples in range. !----------------------------------------------- NEX=INT(PLIMIT/PERIMETER) EXTRAS=EXTRAS+NEX END IF !scan END FOR END FOR PRINT("Primit. with perimeter <=";10#^power;"is";primitives;"&";extras;"non-prim.triples.") PRINT(TIME$)  END FOR   PRINT PRINT("** End **")END PROGRAM
Output:
16:08:39
Primit. with perimeter <= 10 is 0 & 0 non-prim.triples.
16:08:39
Primit. with perimeter <= 100 is 7 & 17 non-prim.triples.
16:08:39
Primit. with perimeter <= 1000 is 70 & 325 non-prim.triples.
16:08:39
Primit. with perimeter <= 10000 is 703 & 4858 non-prim.triples.
16:08:39
Primit. with perimeter <= 100000 is 7026 & 64741 non-prim.triples.
16:08:41
Primit. with perimeter <= 1000000 is 70229 & 808950 non-prim.triples.
16:09:07
Primit. with perimeter <= 10000000 is 702309 & 9706567 non-prim.triples.
16:13:10

** End **

## Euphoria

Translation of: D
function tri(atom lim, sequence in)    sequence r    atom p    p = in[1] + in[2] + in[3]    if p > lim then        return {0, 0}    end if    r = {1, floor(lim / p)}    r += tri(lim, { in[1]-2*in[2]+2*in[3],  2*in[1]-in[2]+2*in[3],  2*in[1]-2*in[2]+3*in[3]})    r += tri(lim, { in[1]+2*in[2]+2*in[3],  2*in[1]+in[2]+2*in[3],  2*in[1]+2*in[2]+3*in[3]})    r += tri(lim, {-in[1]+2*in[2]+2*in[3], -2*in[1]+in[2]+2*in[3], -2*in[1]+2*in[2]+3*in[3]})    return rend function atom max_perimax_peri = 10while max_peri <= 100000000 do    printf(1,"%d: ", max_peri)    ? tri(max_peri, {3, 4, 5})    max_peri *= 10end while

Output:

10: {0,0}
100: {7,17}
1000: {70,325}
10000: {703,4858}
100000: {7026,64741}
1000000: {70229,808950}
10000000: {702309,9706567}
100000000: {7023027,113236940}


## F#

Translation of: OCaml
let isqrt n =    let rec iter t =        let d = n - t*t        if (0 <= d) && (d < t+t+1) // t*t <= n < (t+1)*(t+1)        then t else iter ((t+(n/t))/2)    iter 1 let rec gcd a b =    let t = a % b    if t = 0 then b else gcd b t let coprime a b = gcd a b = 1 let num_to ms =    let mutable ctr = 0    let mutable prim_ctr = 0    let max_m = isqrt (ms/2)    for m = 2 to max_m do        for j = 0 to (m/2) - 1 do            let n = m-(2*j+1)            if coprime m n then                let s = 2*m*(m+n)                if s <= ms then                    ctr <- ctr + (ms/s)                    prim_ctr <- prim_ctr + 1    (ctr, prim_ctr) let show i =    let s, p = num_to i in    printfn "For perimeters up to %d there are %d total and %d primitive" i s p;; List.iter show [ 100; 1000; 10000; 100000; 1000000; 10000000; 100000000 ]
Output:
For perimeters up to 100 there are 17 total and 7 primitive
For perimeters up to 1000 there are 325 total and 70 primitive
For perimeters up to 10000 there are 4858 total and 703 primitive
For perimeters up to 100000 there are 64741 total and 7026 primitive
For perimeters up to 1000000 there are 808950 total and 70229 primitive
For perimeters up to 10000000 there are 9706567 total and 702309 primitive
For perimeters up to 100000000 there are 113236940 total and 7023027 primitive

## Factor

Pretty slow (100 times slower than C)...

Output:
$jq -M -c -r -n -f Pythagorean_triples.jq10: [0,0]100: [17,7]1000: [325,70]10000: [4858,703]100000: [64741,7026]1000000: [808950,70229]10000000: [9706567,702309]100000000: [113236940,7023027]  ## Julia This solution uses the the Euclidian concept of m and n as generators of Pythagorean triplets. When m and n are coprime and have opposite parity, the generated triplets are primitive. It works reasonably well up to a limit of 10^10.  function primitiven{T<:Integer}(m::T) 1 < m || return T[] m != 2 || return T[1] !isprime(m) || return T[2:2:m-1] rp = trues(m-1) if isodd(m) rp[1:2:m-1] = false end for p in keys(factor(m)) rp[p:p:m-1] = false end T[1:m-1][rp]end function pythagoreantripcount{T<:Integer}(plim::T) primcnt = 0 fullcnt = 0 11 < plim || return (primcnt, fullcnt) for m in 2:plim p = 2m^2 p+2m <= plim || break for n in primitiven(m) q = p + 2m*n q <= plim || break primcnt += 1 fullcnt += div(plim, q) end end return (primcnt, fullcnt)end println("Counting Pythagorian Triplets within perimeter limits:")println(" Limit All Primitive")for om in 1:10 (pcnt, fcnt) = pythagoreantripcount(10^om) println(@sprintf " 10^%02d %11d %9d" om fcnt pcnt)end  Output: Counting Pythagorian Triplets within perimeter limits: Limit All Primitive 10^01 0 0 10^02 17 7 10^03 325 70 10^04 4858 703 10^05 64741 7026 10^06 808950 70229 10^07 9706567 702309 10^08 113236940 7023027 10^09 1294080089 70230484 10^10 14557915466 702304875  ## Kotlin Translation of: Go Due to deep recursion, I needed to increase the stack size to 4MB to get up to a maximum perimeter of 10 billion. Expect a run time of around 30 seconds on a typical laptop. // version 1.1.2 var total = 0Lvar prim = 0Lvar maxPeri = 0L fun newTri(s0: Long, s1: Long, s2: Long) { val p = s0 + s1 + s2 if (p <= maxPeri) { prim++ total += maxPeri / p newTri( s0 - 2 * s1 + 2 * s2, 2 * s0 - s1 + 2 * s2, 2 * s0 - 2 * s1 + 3 * s2) newTri( s0 + 2 * s1 + 2 * s2, 2 * s0 + s1 + 2 * s2, 2 * s0 + 2 * s1 + 3 * s2) newTri(-s0 + 2 * s1 + 2 * s2, -2 * s0 + s1 + 2 * s2, -2 * s0 + 2 * s1 + 3 * s2) }} fun main(args: Array<String>) { maxPeri = 100 while (maxPeri <= 10_000_000_000L) { prim = 0 total = 0 newTri(3, 4, 5) println("Up to$maxPeri: $total triples,$prim primatives")        maxPeri *= 10    }}
Output:
Up to 100: 17 triples, 7 primatives
Up to 1000: 325 triples, 70 primatives
Up to 10000: 4858 triples, 703 primatives
Up to 100000: 64741 triples, 7026 primatives
Up to 1000000: 808950 triples, 70229 primatives
Up to 10000000: 9706567 triples, 702309 primatives
Up to 100000000: 113236940 triples, 7023027 primatives
Up to 1000000000: 1294080089 triples, 70230484 primatives
Up to 10000000000: 14557915466 triples, 702304875 primatives


## Lasso

// Brute Force: Too slow for large numbersdefine num_pythagorean_triples(max_perimeter::integer) => {   local(max_b) = (#max_perimeter / 3)*2    return (      with a in 1 to (#max_b - 1)      sum integer(         with b in (#a + 1) to #max_b         let c = math_sqrt(#a*#a + #b*#b)         where #c == integer(#c)         where #c > #b         where (#a+#b+#c) <= #max_perimeter         sum 1      )   )}stdout(Number of Pythagorean Triples in a Perimeter of 100: )stdoutnl(num_pythagorean_triples(100))

Output:

Number of Pythagorean Triples in a Perimeter of 100: 17


## Liberty BASIC

 print time$() for power =1 to 6 perimeterLimit =10^power upperBound =int( 1 +perimeterLimit^0.5) primitives =0 triples =0 extras =0 ' will count the in-range multiples of any primitive for m =2 to upperBound for n =1 +( m mod 2 =1) to m -1 step 2 term1 =2 *m *n term2 =m *m -n *n term3 =m *m +n *n perimeter =term1 +term2 +term3 if perimeter <=perimeterLimit then triples =triples +1 a =term1 b =term2 do r = a mod b a =b b =r loop until r <=0 if ( a =1) and ( perimeter <=perimeterLimit) then 'we've found a primitive triple if a =1, since hcf =1. And it is inside perimeter range. Save it in an array primitives =primitives +1 if term1 >term2 then temp =term1: term1 =term2: term2 =temp 'swap so in increasing order of side length nEx =int( perimeterLimit /perimeter) 'We have the primitive & removed any multiples. Now calculate ALL the multiples in range. extras =extras +nEx end if scan next n next m print " Number of primitives having perimeter below "; 10^power, " was "; primitives, " & "; extras, " non-primitive triples." print time$()next power print "End"end
17:59:34
Number of primitives having perimeter below 10 was 0 & 0 non-primitive triples.
17:59:34
Number of primitives having perimeter below 100 was 7 & 17 non-primitive triples.
17:59:34
Number of primitives having perimeter below 1000 was 70 & 325 non-primitive triples.
17:59:34
Number of primitives having perimeter below 10000 was 703 & 4858 non-primitive triples.
17:59:35
Number of primitives having perimeter below 100000 was 7026 & 64741 non-primitive triples.
17:59:42
Number of primitives having perimeter below 1000000 was 70229 & 808950 non-primitive triples.
18:01:30
End


## Mathematica

Short code but not a very scalable approach...

pythag[n_] := Block[{soln = Solve[{a^2 + b^2 == c^2, a + b + c <= n, 0 < a < b < c}, {a, b, c}, Integers]},        {Length[soln], Count[GCD[a, b] == GCD[b, c] == GCD[c, a] == 1 /. soln, True]}      ]
Output:
pythag[10]
{0,0}

pythag[100]
{17, 7}

pythag[1000]
{325, 70}

## MATLAB / Octave

N=  100;  a = 1:N;  b = a(ones(N,1),:).^2;  b = b+b';  b = sqrt(b);  [y,x]=find(b==fix(b)); % test   % here some alternative tests  % b = b.^(1/k); [y,x]=find(b==fix(b)); % test 2  % [y,x]=find(b==(fix(b.^(1/k)).^k));  % test 3  % b=b.^(1/k); [y,x]=find(abs(b - round(b)) <= 4*eps*b);   z  = sqrt(x.^2+y.^2);  ix = (z+x+y<100) & (x < y) & (y < z);   p  = find(gcd(x(ix),y(ix))==1);   % find primitive triples   printf('There are %i Pythagorean Triples and %i primitive triples with a perimeter smaller than %i.\n',...         sum(ix), length(p), N);

Output:

 There are 17 Pythagorean Triples and 7 primitive triples with a perimeter smaller than 100.

## Mercury

From List comprehensions:

 :- module comprehension.:- interface.:- import_module io.:- import_module int. :- type triple ---> triple(int, int, int). :- pred pythTrip(int::in,triple::out) is nondet.:- pred main(io::di, io::uo) is det. :- implementation.:- import_module solutions. pythTrip(Limit,triple(X,Y,Z)) :-    nondet_int_in_range(1,Limit,X),    nondet_int_in_range(X,Limit,Y),    nondet_int_in_range(Y,Limit,Z),    pow(Z,2) = pow(X,2) + pow(Y,2). main(!IO) :-    solutions((pred(Triple::out) is nondet :- pythTrip(20,Triple)),Result),    write(Result,!IO).

## Modula-3

Note that this code only works on 64bit machines (where INTEGER is 64 bits). Modula-3 provides a LONGINT type, which is 64 bits on 32 bit systems, but there is a bug in the implementation apparently.

MODULE PyTriple64 EXPORTS Main; IMPORT IO, Fmt; VAR tcnt, pcnt, max, i: INTEGER; PROCEDURE NewTriangle(a, b, c: INTEGER; VAR tcount, pcount: INTEGER) =  VAR perim := a + b + c;        BEGIN    IF perim <= max THEN      pcount := pcount + 1;      tcount := tcount + max DIV perim;      NewTriangle(a-2*b+2*c, 2*a-b+2*c, 2*a-2*b+3*c, tcount, pcount);      NewTriangle(a+2*b+2*c, 2*a+b+2*c, 2*a+2*b+3*c, tcount, pcount);      NewTriangle(2*b+2*c-a, b+2*c-2*a, 2*b+3*c-2*a, tcount, pcount);    END;  END NewTriangle; BEGIN  i := 100;   REPEAT    max := i;    tcnt := 0;    pcnt := 0;    NewTriangle(3, 4, 5, tcnt, pcnt);    IO.Put(Fmt.Int(i) & ": " & Fmt.Int(tcnt) & " Triples, " &      Fmt.Int(pcnt) & " Primitives\n");    i := i * 10;  UNTIL i = 10000000;END PyTriple64.

Output:

100: 17 Triples, 7 Primitives
1000: 325 Triples, 70 Primitives
10000: 4858 Triples, 703 Primitives
100000: 64741 Triples, 7026 Primitives
1000000: 808950 Triples, 70229 Primitives


## Nanoquery

### Brute force method

Translation of: Java
import math // a function to check if three numbers are a valid tripledef is_triple(a, b, c)    if not (a < b) and (b < c)        return false    end     return (a^2 + b^2) = c^2end // a function to check if the numbers are coprimedef is_coprime(a, b, c)    global math    return (math.gcd(a, b)=1) && (math.gcd(a, c)=1) && (math.gcd(b, c)=1)end // the maximum perimeter to checkperimeter  = 100perimeter2 = int(perimeter / 2) - 1perimeter3 = int(perimeter / 3) - 1 // loop though and look for pythagorean triplests = 0ps = 0for a in range(1, perimeter3)    for b in range(a + 1, perimeter2)        for c in range(b + 1, perimeter2)            if (a + b + c) <= perimeter                if is_triple(a,b,c)                    ts += 1                    print a + ", " + b + ", " + c                    if is_coprime(a,b,c)                        ps += 1                        print " primitive"                    end                    println                end            end        end    endend print   "Up to a perimeter of " + perimeter + ", there are " + tsprintln " triples, of which " + ps + " are primitive."
Output:
3, 4, 5 primitive
5, 12, 13 primitive
6, 8, 10
7, 24, 25 primitive
8, 15, 17 primitive
9, 12, 15
9, 40, 41 primitive
10, 24, 26
12, 16, 20
12, 35, 37 primitive
15, 20, 25
15, 36, 39
16, 30, 34
18, 24, 30
20, 21, 29 primitive
21, 28, 35
24, 32, 40
Up to a perimeter of 100, there are 17 triples, of which 7 are primitive.

## Nim

Translation of: C
const u = [[ 1, -2,  2,  2, -1,  2,  2, -2,  3],           [ 1,  2,  2,  2,  1,  2,  2,  2,  3],           [-1,  2,  2, -2,  1,  2, -2,  2,  3]] var  total, prim = 0  maxPeri = 10 proc newTri(ins: array[0..2, int]) =  var p = ins[0] + ins[1] + ins[2]  if p > maxPeri: return  inc(prim)  total += maxPeri div p   for i in 0..2:    newTri([u[i][0] * ins[0] + u[i][1] * ins[1] + u[i][2] * ins[2],            u[i][3] * ins[0] + u[i][4] * ins[1] + u[i][5] * ins[2],            u[i][6] * ins[0] + u[i][7] * ins[1] + u[i][8] * ins[2]]) while maxPeri <= 100_000_000:  total = 0  prim = 0  newTri([3, 4, 5])  echo "Up to ", maxPeri, ": ", total, " triples, ", prim, " primitives"  maxPeri *= 10

Output:

Up to 10: 0 triples, 0 primitives
Up to 100: 17 triples, 7 primitives
Up to 1000: 325 triples, 70 primitives
Up to 10000: 4858 triples, 703 primitives
Up to 100000: 64741 triples, 7026 primitives
Up to 1000000: 808950 triples, 70229 primitives
Up to 10000000: 9706567 triples, 702309 primitives
Up to 100000000: 113236940 triples, 7023027 primitives

## OCaml

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 ]

Output:

For perimeters up to 100 there are 17 total and 7 primitive
For perimeters up to 1000 there are 325 total and 70 primitive
For perimeters up to 10000 there are 4858 total and 703 primitive
For perimeters up to 100000 there are 64741 total and 7026 primitive
For perimeters up to 1000000 there are 808950 total and 70229 primitive
For perimeters up to 10000000 there are 9706567 total and 702309 primitive
For perimeters up to 100000000 there are 113236940 total and 7023027 primitive

## Ol

 ; triples generator based on Euclid's formula, creates lazy list(define (euclid-formula max)   (let loop ((a 3) (b 4) (c 5) (tail #null))      (if (<= (+ a b c) max)         (cons (tuple a b c) (lambda ()            (let ((d (- b)) (z (- a)))            (loop (+ a d d c c) (+ a a d c c) (+ a a d d c c c) (lambda ()               (loop (+ a b b c c) (+ a a b c c) (+ a a b b c c c) (lambda ()                  (loop (+ z b b c c) (+ z z b c c) (+ z z b b c c c) tail))))))))         tail))) ; let's do calculations(define (calculate max)   (let loop ((p 0) (t 0) (ll (euclid-formula max)))      (cond         ((null? ll)            (cons p t))         ((function? ll)            (loop p t (ll)))         (else            (let ((triple (car ll)))               (loop (+ p 1) (+ t (div max (apply + triple)))               (cdr ll))))))) ; print values for 10..100000(for-each (lambda (max)      (print max ": " (calculate max)))   (map (lambda (n) (expt 10 n)) (iota 6 1)))
Output:
10: (0 . 0)
100: (7 . 17)
1000: (70 . 325)
10000: (703 . 4858)
100000: (7026 . 64741)
1000000: (70229 . 808950)


## PARI/GP

This version is reasonably efficient and can handle inputs like a million quickly.

do(lim)={  my(prim,total,P);  lim\=1;  for(m=2,sqrtint(lim\2),    forstep(n=1+m%2,min(sqrtint(lim-m^2),m-1),2,      P=2*m*(m+n);      if(gcd(m,n)==1 && P<=lim,        prim++;        total+=lim\P      )    )  );  [prim,total]};do(100)

## Pascal

Program PythagoreanTriples (output); var  total, prim, maxPeri: int64; procedure newTri(s0, s1, s2: int64);  var    p: int64;  begin    p := s0 + s1 + s2;    if p <= maxPeri then    begin      inc(prim);      total := total + maxPeri div p;      newTri( s0 + 2*(-s1+s2),  2*( s0+s2) - s1,  2*( s0-s1+s2) + s2);      newTri( s0 + 2*( s1+s2),  2*( s0+s2) + s1,  2*( s0+s1+s2) + s2);      newTri(-s0 + 2*( s1+s2),  2*(-s0+s2) + s1,  2*(-s0+s1+s2) + s2);    end;  end; begin  maxPeri := 100;  while maxPeri <= 1e10 do  begin    prim := 0;    total := 0;    newTri(3, 4, 5);    writeln('Up to ', maxPeri, ': ', total, ' triples, ', prim, ' primitives.');    maxPeri := maxPeri * 10;  end;end.

Output (on Core2Duo 2GHz laptop):

time ./PythagoreanTriples
Up to 100: 17 triples, 7 primitives.
Up to 1000: 325 triples, 70 primitives.
Up to 10000: 4858 triples, 703 primitives.
Up to 100000: 64741 triples, 7026 primitives.
Up to 1000000: 808950 triples, 70229 primitives.
Up to 10000000: 9706567 triples, 702309 primitives.
Up to 100000000: 113236940 triples, 7023027 primitives.
Up to 1000000000: 1294080089 triples, 70230484 primitives.
Up to 10000000000: 14557915466 triples, 702304875 primitives.
109.694u 0.094s 1:50.43 99.4%   0+0k 0+0io 0pf+0w

## Perl

sub gcd {    my ($n,$m) = @_;    while($n){ my$t = $n;$n = $m %$n;        $m =$t;    }    return $m;} sub tripel { my$pmax  = shift;    my $prim = 0; my$count = 0;    my $nmax = sqrt($pmax)/2;    for( my $n=1;$n<=$nmax;$n++ ) {        for( my $m=$n+1; (my $p = 2*$m*($m+$n)) <= $pmax;$m+=2 ) {            next unless 1==gcd($m,$n);            $prim++;$count += int $pmax/$p;        }    }    printf "Max. perimeter: %d, Total: %d, Primitive: %d\n", $pmax,$count, $prim;} tripel 10**$_ for 1..8;
Output:
Max. perimeter: 10, Total: 0, Primitive: 0
Max. perimeter: 100, Total: 17, Primitive: 7
Max. perimeter: 1000, Total: 325, Primitive: 70
Max. perimeter: 10000, Total: 4858, Primitive: 703
Max. perimeter: 100000, Total: 64741, Primitive: 7026
Max. perimeter: 1000000, Total: 808950, Primitive: 70229
Max. perimeter: 10000000, Total: 9706567, Primitive: 702309
Max. perimeter: 100000000, Total: 113236940, Primitive: 7023027


## Perl 6

Works with: Rakudo version 2018.09

Here is a straight-forward, naive brute force implementation:

constant limit = 100; for [X] [^limit] xx 3 -> (\a, \b, \c) {    say [a, b, c] if a < b < c and a**2 + b**2 == c**2}
Output:
3 4 5
5 12 13
6 8 10
7 24 25
8 15 17
9 12 15
9 40 41
10 24 26
11 60 61
12 16 20
12 35 37
13 84 85
14 48 50
15 20 25
15 36 39
16 30 34
16 63 65
18 80 82
20 21 29
20 48 52
21 28 35
21 72 75
24 32 40
24 45 51
24 70 74
25 60 65
27 36 45
28 45 53
30 40 50
30 72 78
32 60 68
33 44 55
33 56 65
35 84 91
36 48 60
36 77 85
39 52 65
39 80 89
40 42 58
40 75 85
42 56 70
45 60 75
48 55 73
48 64 80
51 68 85
54 72 90
57 76 95
60 63 87
65 72 97

Here is a slightly less naive brute force implementation, but still not practical for large perimeter limits.

my $limit = 10000;my atomicint$i = 0;my @triples[$limit/2];(3 ..$limit/2).race.map: -> $c { for 1 ..$c -> $a { my$b = ($c *$c - $a *$a).sqrt;       last if $c +$a + $b >$limit;       last if $a >$b;       @triples[$i⚛++] = ([gcd]$c, $a,$b) > 1 ?? 0 !! 1 if $b ==$b.Int;   }} say my $result = "There are {[email protected]:{$_ !eqv Any}} Pythagorean Triples with a perimeter <= $limit," ~"\nof which {[+] @triples.grep: so *} are primitive."; Output: There are 4858 Pythagorean Triples with a perimeter <= 10000, of which 703 are primitive. Here's a much faster version. Hint, "oyako" is Japanese for "parent/child". :-) sub triples($limit) {    my $primitive = 0; my$civilized = 0;     sub oyako($a,$b, $c) { my$perim = $a +$b + $c; return if$perim > $limit; ++$primitive; $civilized +=$limit div $perim; oyako($a - 2*$b + 2*$c,  2*$a -$b + 2*$c, 2*$a - 2*$b + 3*$c);        oyako( $a + 2*$b + 2*$c, 2*$a + $b + 2*$c,  2*$a + 2*$b + 3*$c); oyako(-$a + 2*$b + 2*$c, -2*$a +$b + 2*$c, -2*$a + 2*$b + 3*$c);    }     oyako(3,4,5);    "$limit => ($primitive $civilized)";} for 10,100,1000 ... * ->$limit {    say triples $limit;} Output: 10 => (0 0) 100 => (7 17) 1000 => (70 325) 10000 => (703 4858) 100000 => (7026 64741) 1000000 => (70229 808950) 10000000 => (702309 9706567) 100000000 => (7023027 113236940) 1000000000 => (70230484 1294080089) ^C The geometric sequence of limits will continue on forever, so eventually when you get tired of waiting (about a billion on my computer), you can just stop it. Another efficiency trick of note: we avoid passing the limit in as a parameter to the inner helper routine, but instead supply the limit via the lexical scope. Likewise, the accumulators are referenced lexically, so only the triples themselves need to be passed downward, and nothing needs to be returned. Here is an alternate version that avoids naming any scalars that can be handled by vector processing instead. Using vectorized ops allows a bit more potential for parallelization in theory, but techniques like the use of complex numbers to add two numbers in parallel, and the use of gather/take generate so much overhead that this version runs 70-fold slower than the previous one. constant @coeff = [[+1, -2, +2], [+2, -1, +2], [+2, -2, +3]], [[+1, +2, +2], [+2, +1, +2], [+2, +2, +3]], [[-1, +2, +2], [-2, +1, +2], [-2, +2, +3]]; sub triples($limit) {     sub oyako(@trippy) {        my $perim = [+] @trippy; return if$perim > $limit; take (1 + ($limit div $perim)i); for @coeff -> @nine { oyako (map -> @three { [+] @three »*« @trippy }, @nine); } return; } my$complex = 0i + [+] gather oyako([3,4,5]);    "$limit => ({$complex.re, $complex.im})";} for 10, 100, 1000, 10000 ->$limit {    say triples $limit;} Output: 10 => (0 0) 100 => (7 17) 1000 => (70 325) 10000 => (703 4858) ## Phix Translation of: Pascal atom total, prim, maxPeri = 10 procedure tri(atom s0, s1, s2)atom p = s0 + s1 + s2 if p<=maxPeri then prim += 1 total += floor(maxPeri/p) tri( s0+2*(-s1+s2), 2*( s0+s2)-s1, 2*( s0-s1+s2)+s2); tri( s0+2*( s1+s2), 2*( s0+s2)+s1, 2*( s0+s1+s2)+s2); tri(-s0+2*( s1+s2), 2*(-s0+s2)+s1, 2*(-s0+s1+s2)+s2); end ifend procedure while maxPeri<=1e8 do prim := 0; total := 0; tri(3, 4, 5); printf(1,"Up to %d: %d triples, %d primitives.\n", {maxPeri,total,prim}) maxPeri *= 10;end while Output: Up to 10: 0 triples, 0 primitives. Up to 100: 17 triples, 7 primitives. Up to 1000: 325 triples, 70 primitives. Up to 10000: 4858 triples, 703 primitives. Up to 100000: 64741 triples, 7026 primitives. Up to 1000000: 808950 triples, 70229 primitives. Up to 10000000: 9706567 triples, 702309 primitives. Up to 100000000: 113236940 triples, 7023027 primitives.  ## PHP <?php function gcd($a, $b){ if ($a == 0)       return $b; if ($b == 0)       return $a; if($a == $b) return$a;    if($a >$b)        return gcd($a-$b, $b); return gcd($a, $b-$a);} $pytha = 0;$prim = 0;$max_p = 100; for ($a = 1; $a <=$max_p / 3; $a++) {$aa = $a**2; for ($b = $a + 1;$b < $max_p/2;$b++) {        $bb =$b**2;        for ($c =$b + 1; $c <$max_p/2; $c++) {$cc = $c**2; if ($aa + $bb <$cc) break;            if ($a +$b + $c >$max_p) break;             if ($aa +$bb == $cc) {$pytha++;                if (gcd($a,$b) == 1) $prim++; } } }} echo 'Up to ' .$max_p . ', there are ' . $pytha . ' triples, of which ' .$prim . ' are primitive.';
Output:
Up to 100, there are 17 triples, of which 7 are primitive.

## PicoLisp

Translation of: C
(for (Max 10  (>= 100000000 Max)  (* Max 10))   (let (Total 0  Prim 0  In (3 4 5))      (recur (In)         (let P (apply + In)            (when (>= Max P)               (inc 'Prim)               (inc 'Total (/ Max P))               (for Row                  (quote                     (( 1 -2 2) ( 2 -1 2) ( 2 -2 3))                     (( 1  2 2) ( 2  1 2) ( 2  2 3))                     ((-1  2 2) (-2  1 2) (-2  2 3)) )                  (recurse                     (mapcar '((U) (sum * U In)) Row) ) ) ) ) )      (prinl "Up to " Max ": " Total " triples, " Prim " primitives.") ) )

Output:

Up to 10: 0 triples, 0 primitives.
Up to 100: 17 triples, 7 primitives.
Up to 1000: 325 triples, 70 primitives.
Up to 10000: 4858 triples, 703 primitives.
Up to 100000: 64741 triples, 7026 primitives.
Up to 1000000: 808950 triples, 70229 primitives.
Up to 10000000: 9706567 triples, 702309 primitives.
Up to 100000000: 113236940 triples, 7023027 primitives.

## PL/I

Version 1

*process source attributes xref or(!); /********************************************************************* * REXX pgm counts number of  Pythagorean triples * that exist given a max perimeter of  N, * and also counts how many of them are primatives. * 05.05.2013 Walter Pachl  translated from REXX version 2 *********************************************************************/ pyt: Proc Options(main); Dcl sysprint Print; Dcl (addr,mod,right) Builtin; Dcl memn Bin Fixed(31) Init(0); Dcl mabca(300) Char(12); Dcl 1 mabc,      2 ma Dec fixed(7),      2 mb Dec fixed(7),      2 mc Dec fixed(7); Dcl mabce Char(12) Based(addr(mabc)); Dcl 1 abc,      2 a Dec fixed(7),      2 b Dec fixed(7),      2 c Dec fixed(7); Dcl abce Char(12) Based(addr(abc)); Dcl (prims,trips,m,n,aa,aabb,cc,aeven,ab) Dec Fixed(7); mabca=''; trips=0; prims=0; n=100; la: Do a=3 To n/3;    aa=a*a;                       /* limit side to 1/3 of perimeter.*/    aeven=mod(a,2)=0; lb:Do b=a+1 By 1+aeven;          /* triangle can't be isosceles.   */      ab=a+b;                     /* compute partial perimeter.     */      If ab>=n Then        Iterate la;               /* a+b>perimeter?  Try different A*/      aabb=aa+b*b;                /* compute sum of  a² + b² (cheat)*/      Do c=b+1 By 1;        cc=c*c;                   /* 3rd side:   also compute  c²   */        If aeven Then          If mod(c,2)=0 Then            Iterate;        If ab+c>n Then          Iterate la;              /* a+b+c > perimeter?  Try diff A.*/        If cc>aabb Then          Iterate lb;              /* c² >  a²+b² ?  Try different B.*/        If cc^=aabb Then          Iterate;                 /* c² ¬= a²+b² ?  Try different C.*/        If mema(abce) Then          Iterate;        trips=trips+1;             /* eureka.                        */        prims=prims+1;             /* count this primitive triple.   */        Put Edit(a,b,c,'   ',right(a**2+b**2,5),right(c**2,5),a+b+c)                (Skip,f(4),2(f(5)),a,2(f(6)),f(9));        Do m=2 By 1;          ma=a*m;          mb=b*m;          mc=c*m;                  /* gen non-primitives.            */          If ma+mb+mc>n Then            Leave;                                   /* is this multiple a triple ?    */          trips=trips+1;           /* yuppers, then we found another.*/          If mod(m,2)=1 Then       /* store as even multiple.        */            call mems(mabce);          Put Edit(ma,mb,mc,' * ',                          right(ma**2+mb**2,5),right(mc**2,5),ma+mb+mc)                (Skip,f(4),2(f(5)),a,2(f(6)),f(9));          End;                     /* m                              */        End;                       /* c                              */      End;                         /* b                              */    End;                           /* a                              */  Put Edit('max perimeter = ',n,   /* show a single line of output.  */           'Pythagorean triples =',trips,           'primitives =',prims)          (Skip,a,f(5),2(x(9),a,f(4)));  mems: Proc(e); Dcl e Char(12); memn+=1; mabca(memn)=e; End;  mema: Proc(e) Returns(bit(1)); Dcl e Char(12); Do memi=1 To memn;   If mabca(memi)=e Then Return('1'b);   End; Return('0'b); End;  End;

Version 2

 pythagorean: procedure options (main, reorder); /* 23 January 2014 */   declare (a, b, c) fixed (3);   declare (asquared, bsquared) fixed;   declare (triples, primitives) fixed binary(31) initial (0);    do a = 1 to 100;      asquared = a*a;      do b = a+1 to 100;         bsquared = b*b;         do c = b+1 to 100;            if a+b+c <= 100 then             if asquared + bsquared = c*c then               do;                  triples = triples + 1;                  if GCD(a,b) = 1 then primitives = primitives + 1;               end;         end;      end;   end;   put skip data (triples, primitives);  GCD: procedure (a, b) returns (fixed binary (31)) recursive;   declare (a, b) fixed binary (31);    if b = 0 then return (a);    return (GCD (b, mod(a, b)) ); end GCD; end pythagorean;

Output:

TRIPLES=            17  PRIMITIVES=             7;


Output for P=1000:

TRIPLES=           325  PRIMITIVES=            70;


## Scheme

Works with: Gauche Scheme
(use srfi-42) (define (py perim)  (define prim 0)  (values    (sum-ec      (: c perim) (: b c) (: a b)      (if (and (<= (+ a b c) perim)               (= (square c) (+ (square b) (square a)))))      (begin (when (= 1 (gcd a b)) (inc! prim)))      1)    prim))

Testing:

gosh> (py 100)
17
7


## Scratch

Scratch is a visual programming language. Click the link, then "see inside" to see the code.

Output: 17 Pythagorean triples with a perimeter less than 100, 7 of which are primitive.

## Seed7

Translation of: C efficient method

The example below uses bigInteger numbers:

$include "seed7_05.s7i"; include "bigint.s7i"; var bigInteger: total is 0_;var bigInteger: prim is 0_;var bigInteger: max_peri is 10_; const proc: new_tri (in bigInteger: a, in bigInteger: b, in bigInteger: c) is func local var bigInteger: p is 0_; begin p := a + b + c; if p <= max_peri then incr(prim); total +:= max_peri div p; new_tri( a - 2_*b + 2_*c, 2_*a - b + 2_*c, 2_*a - 2_*b + 3_*c); new_tri( a + 2_*b + 2_*c, 2_*a + b + 2_*c, 2_*a + 2_*b + 3_*c); new_tri(-a + 2_*b + 2_*c, -2_*a + b + 2_*c, -2_*a + 2_*b + 3_*c); end if; end func; const proc: main is func begin while max_peri <= 100000000_ do total := 0_; prim := 0_; new_tri(3_, 4_, 5_); writeln("Up to " <& max_peri <& ": " <& total <& " triples, " <& prim <& " primitives."); max_peri *:= 10_; end while; end func; Output: Up to 10: 0 triples, 0 primitives. Up to 100: 17 triples, 7 primitives. Up to 1000: 325 triples, 70 primitives. Up to 10000: 4858 triples, 703 primitives. Up to 100000: 64741 triples, 7026 primitives. Up to 1000000: 808950 triples, 70229 primitives. Up to 10000000: 9706567 triples, 702309 primitives. Up to 100000000: 113236940 triples, 7023027 primitives.  ## Sidef Translation of: Perl 6 func triples(limit) { var primitive = 0 var civilized = 0 func oyako(a, b, c) { (var perim = a+b+c) > limit || ( primitive++ civilized += int(limit / perim) oyako( a - 2*b + 2*c, 2*a - b + 2*c, 2*a - 2*b + 3*c) oyako( a + 2*b + 2*c, 2*a + b + 2*c, 2*a + 2*b + 3*c) oyako(-a + 2*b + 2*c, -2*a + b + 2*c, -2*a + 2*b + 3*c) ) } oyako(3,4,5) "#{limit} => (#{primitive} #{civilized})"} for n (1..Inf) { say triples(10**n)} Output: 10 => (0 0) 100 => (7 17) 1000 => (70 325) 10000 => (703 4858) 100000 => (7026 64741) 1000000 => (70229 808950) ^C  ## Swift Translation of: Pascal var total = 0var prim = 0var maxPeri = 100 func newTri(s0:Int, _ s1:Int, _ s2: Int) -> () { let p = s0 + s1 + s2 if p <= maxPeri { prim += 1 total += maxPeri / p newTri( s0 + 2*(-s1+s2), 2*( s0+s2) - s1, 2*( s0-s1+s2) + s2) newTri( s0 + 2*( s1+s2), 2*( s0+s2) + s1, 2*( s0+s1+s2) + s2) newTri(-s0 + 2*( s1+s2), 2*(-s0+s2) + s1, 2*(-s0+s1+s2) + s2) }} while maxPeri <= 100_000_000 { prim = 0 total = 0 newTri(3, 4, 5) print("Up to \(maxPeri) : \(total) triples \( prim) primitives.") maxPeri *= 10} Output: Up to 100 : 17 triples 7 primitives. Up to 1000 : 325 triples 70 primitives. Up to 10000 : 4858 triples 703 primitives. Up to 100000 : 64741 triples 7026 primitives. Up to 1000000 : 808950 triples 70229 primitives. Up to 10000000 : 9706567 triples 702309 primitives. Up to 100000000 : 113236940 triples 7023027 primitives.  ## Tcl Using the efficient method based off the Wikipedia article: proc countPythagoreanTriples {limit} { lappend q 3 4 5 set idx [set count [set prim 0]] while {$idx < [llength $q]} { set a [lindex$q $idx] set b [lindex$q [incr idx]]    set c [lindex $q [incr idx]] incr idx if {$a + $b +$c <= $limit} { incr prim for {set i 1} {$i*$a+$i*$b+$i*$c <=$limit} {incr i} {        incr count        }        lappend q \        [expr {$a + 2*($c-$b)}] [expr {2*($a+$c) -$b}] [expr {2*($a-$b) + 3*$c}] \ [expr {$a + 2*($b+$c)}] [expr {2*($a+$c) + $b}] [expr {2*($a+$b) + 3*$c}] \        [expr {2*($b+$c) - $a}] [expr {2*($c-$a) +$b}] [expr {2*($b-$a) + 3*$c}] } } return [list$count $prim]}for {set i 10} {$i <= 10000000} {set i [expr {$i*10}]} { lassign [countPythagoreanTriples$i] count primitive    puts "perimeter limit $i =>$count triples, \$primitive primitive"}

Output:

perimeter limit 10 => 0 triples, 0 primitive
perimeter limit 100 => 17 triples, 7 primitive
perimeter limit 1000 => 325 triples, 70 primitive
perimeter limit 10000 => 4858 triples, 703 primitive
perimeter limit 100000 => 64741 triples, 7026 primitive
perimeter limit 1000000 => 808950 triples, 70229 primitive
perimeter limit 10000000 => 9706567 triples, 702309 primitive


## VBA

Translation of: Pascal
Dim total As Variant, prim As Variant, maxPeri As VariantPrivate Sub newTri(s0 As Variant, s1 As Variant, s2 As Variant)    Dim p As Variant    p = CDec(s0) + CDec(s1) + CDec(s2)    If p <= maxPeri Then        prim = prim + 1        total = total + maxPeri \ p        newTri s0 + 2 * (-s1 + s2), 2 * (s0 + s2) - s1, 2 * (s0 - s1 + s2) + s2        newTri s0 + 2 * (s1 + s2), 2 * (s0 + s2) + s1, 2 * (s0 + s1 + s2) + s2        newTri -s0 + 2 * (s1 + s2), 2 * (-s0 + s2) + s1, 2 * (-s0 + s1 + s2) + s2      End IfEnd SubPublic Sub Program_PythagoreanTriples()    maxPeri = CDec(100)    Do While maxPeri <= 10000000#        prim = CDec(0)        total = CDec(0)        newTri 3, 4, 5        Debug.Print "Up to "; maxPeri; ": "; total; " triples, "; prim; " primitives."        maxPeri = maxPeri * 10    LoopEnd Sub
Output:
Up to  100 :  17  triples,  7  primitives.
Up to  1000 :  325  triples,  70  primitives.
Up to  10000 :  4858  triples,  703  primitives.
Up to  100000 :  64741  triples,  7026  primitives.
Up to  1000000 :  808950  triples,  70229  primitives.
Up to  10000000 :  9706567  triples,  702309  primitives.


## VBScript

Translation of: Perl
 For i=1 To 8	WScript.StdOut.WriteLine triples(10^i)Next Function triples(pmax)	prim=0 : count=0 : nmax=Sqr(pmax)/2 : n=1	Do While n <= nmax		m=n+1 : p=2*m*(m+n)		Do While p <= pmax			If gcd(m,n)=1 Then				prim=prim+1				count=count+Int(pmax/p)			End If			m=m+2			p=2*m*(m+n)		Loop		n=n+1	Loop 	triples = "Max Perimeter: " & pmax &_				", Total: " & count &_				", Primitive: " & primEnd Function Function gcd(a,b)	c = a : d = b	Do		If c Mod d > 0 Then			e = c Mod d			c = d			d = e		Else			gcd = d			Exit Do		End If	LoopEnd Function

## Visual Basic

Translation of: VBA
Works with: Visual Basic version 5
Works with: Visual Basic version 6
Works with: VBA version Access 97
Works with: VBA version 6.5
Works with: VBA version 7.1
Option Explicit Dim total As Long, prim As Long, maxPeri As Long Public Sub NewTri(ByVal s0 As Long, ByVal s1 As Long, ByVal s2 As Long)Dim p As Long, x1 As Long, x2 As Long    p = s0 + s1 + s2    If p <= maxPeri Then        prim = prim + 1        total = total + maxPeri \ p        x1 = s0 + s2        x2 = s1 + s2        NewTri s0 + 2 * (-s1 + s2), 2 * x1 - s1, 2 * (x1 - s1) + s2        NewTri s0 + 2 * x2, 2 * x1 + s1, 2 * (x1 + s1) + s2        NewTri -s0 + 2 * x2, 2 * (-s0 + s2) + s1, 2 * (-s0 + x2) + s2    End IfEnd Sub Public Sub Main()    maxPeri = 100    Do While maxPeri <= 10& ^ 8        prim = 0        total = 0        NewTri 3, 4, 5        Debug.Print "Up to "; maxPeri; ": "; total; " triples, "; prim; " primitives."        maxPeri = maxPeri * 10    LoopEnd Sub
Output:
Up to  100 :  17  triples,  7  primitives.
Up to  1000 :  325  triples,  70  primitives.
Up to  10000 :  4858  triples,  703  primitives.
Up to  100000 :  64741  triples,  7026  primitives.
Up to  1000000 :  808950  triples,  70229  primitives.
Up to  10000000 :  9706567  triples,  702309  primitives.
Up to  100000000 :  113236940  triples,  7023027  primitives.

## zkl

Translation of: D
fcn tri(lim,a=3,b=4,c=5){    p:=a + b + c;    if(p>lim) return(0,0);    T(1,lim/p).zipWith('+,       tri(lim,  a - 2*b + 2*c,  2*a - b + 2*c,  2*a - 2*b + 3*c),       tri(lim,  a + 2*b + 2*c,  2*a + b + 2*c,  2*a + 2*b + 3*c),       tri(lim, -a + 2*b + 2*c, -2*a + b + 2*c, -2*a + 2*b + 3*c)    );}
n:=10; do(10){ println("%,d: %s".fmt(n,tri(n).reverse())); n*=10; }
Output:
10: L(0,0)
100: L(17,7)
1,000: L(325,70)
10,000: L(4858,703)
100,000: L(64741,7026)
1,000,000: L(808950,70229)
10,000,000: L(9706567,702309)
VM#1 caught this unhandled exception:
AssertionError : That is one big stack, infinite recursion?
Stack trace for VM#1 ():
<repeats 3578 times>
...


Max stack size is arbitrary but not adjustable.

## ZX Spectrum Basic

ZX Spectrum: 8 bit microprocessor 3.5 Mhz doing all the work. In an effort to get some decent speed the program is made to be as fast as it can.

It takes about 90 seconds for limit = 10 000 and 17 minutes for limit=100 000 and 3.5 hours for limit = 1000 000.

To set the limits. Set in line nr: 1 L to the starting limit. Set in line nr: 11 IF L<=(last limit to calculate)

Ex. start as limit 100 and end on limit 1000. Set in line nr: 1 LET L=100. Set in line nr: 11 IF L<=1000 THEN GO TO 2

   1 LET Y=0: LET X=0: LET Z=0: LET V=0: LET U=0: LET L=10: LET T=0: LET P=0: LET N=4: LET M=0: PRINT "limit   trip.   prim."   2 FOR U=2 TO INT (SQR (L/2)): LET Y=U-INT (U/2)*2: LET N=N+4: LET M=U*U*2: IF Y=0 THEN LET M=M-U-U   3 FOR V=1+Y TO U-1 STEP 2: LET M=M+N: LET X=U: LET Y=V   4 LET Z=Y: LET Y=X-INT (X/Y)*Y: LET X=Z: IF Y<>0 THEN GO TO 4   5 IF X>1 THEN GO TO 8   6 IF M>L THEN GO TO 9   7 LET P=P+1: LET T=T+INT (L/M)   8 NEXT V   9 NEXT U  10 PRINT L;TAB 8;T;TAB 16;P  11 LET N=4: LET T=0: LET P=0: LET L=L*10: IF L<=100000 THEN GO TO 2
Output:
limit   trip.   prim.
10      0       0
100     17      7
1000    325     70
10000   4858    703
100000  64741   7026
1000000 808950  70229