Source code in the pl/i programming language

*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;

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;