(subscriptrange):
Hamming: procedure options (main); /* 14 November 2013 with fixes 2021 */
   declare (H(2000), p2, p3, p5, twoTo31, Hm, tenP(11)) decimal(12)fixed;
   declare (i, j, k, m, d, w) fixed binary;

   /* Quicksorts in-place the array of integers H, from lb to ub */
   quicksortH: procedure( lb, ub ) recursive;
      declare ( lb,   ub    )binary(15)fixed;
      declare ( left, right )binary(15)fixed;
      declare ( pivot, swap )decimal(12)fixed;
      declare sorting        bit(1);
      if ub > lb then do
         /* more than one element, so must sort */
         left    = lb;
         right   = ub;
         /* choosing the middle element of the array as the pivot */
         pivot   = H( left + ( ( right + 1 ) - left ) / 2 );
         sorting = '1'b;
         do while( sorting );
            do while( left  <= ub & H( left  ) < pivot ); left  = left  + 1; end;
            do while( right >= lb & H( right ) > pivot ); right = right - 1; end;
            sorting = ( left <= right );
            if sorting then do;
               swap       = H( left  );
               H( left  ) = H( right );
               H( right ) = swap;
               left       = left  + 1;
               right      = right - 1;
            end;
         end;
         call quicksortH( lb,   right );
         call quicksortH( left, ub    );
      end;
   end quicksortH ;

   /* find 2^31 - the limit for Hamming numbers we need to find */
   twoTo31 = 2;
   do i = 2 to 31;
      twoTo31 = twoTo31 * 2;
   end;
   /* calculate powers of 10 so we can check the number of digits */
   /* the numbers will have */
   tenP( 1 ) = 10;
   do i = 2 to 11;
      tenP( i ) = 10 * tenP( i - 1 );
   end;

   /* find the numbers */
   m = 0;
   p5 = 1;
   do k = 0 to 13;
      p3 = 1;
      do j = 0 to 19;
         Hm = 0;
         p2 = 1;
         do i = 0 to 31 while( Hm < twoTo31 );
            /* count the number of digits p2 * p3 * p5 will have */
            d = 0;
            do w = 1 to 11 while( tenP(w) < p2 ); d = d + 1; end;
            do w = 1 to 11 while( tenP(w) < p3 ); d = d + 1; end;
            do w = 1 to 11 while( tenP(w) < p5 ); d = d + 1; end;
            if d < 11 then do;
               /* the product will be small enough */
               Hm = p2 * p3 * p5;
               if Hm < twoTo31 then do;
                  m = m + 1;
                  H(m) = Hm;
               end;
            end;
            p2 = p2 * 2;
         end;
         p3 = p3 * 3;
      end;
      p5 = p5 * 5;
   end;

   /* sort the numbers */
   call quicksortH( 1, m );

   put skip list( 'The first 20 Hamming numbers:' );
   do i = 1 to 20;
      put skip list (H(i));
   end;
   put skip list( 'Hamming number 1691:' );
   put skip list (H(1691));

end Hamming;