/* Arithmetic numbers */

arithm_nums:
  procedure options(main);
  declare
    (n, arithm_cnt, comp_cnt) fixed binary(15),
    (dv, dv_cnt, sum, quot) fixed binary(15);
  n = 1;
  arithm_cnt = 0;
  comp_cnt = 0;
  put list('The first 100 arithmetic numbers are:');
  put skip;
  do while (arithm_cnt < 1001);
    dv = 1;
    dv_cnt = 0;
    sum = 0;
    do while ('1'b);
      quot = n / dv;
      if quot < dv then
        goto done;
      if quot = dv & mod(n, dv) = 0 then do; /* n is a square */
        sum = sum + quot;
        dv_cnt = dv_cnt + 1;
        goto done;
      end;
      if mod(n, dv) = 0 then do;
        sum = sum + dv + quot;
        dv_cnt = dv_cnt + 2;
      end;
      dv = dv + 1;
    end;
    done:
    if mod(sum, dv_cnt) = 0 then do; /* n is arithmetic */
      arithm_cnt = arithm_cnt + 1;
      if arithm_cnt <= 100 then do;
        put edit(n)(f(4));
        if mod(arithm_cnt, 10) = 0 then
          put skip;
      end;
      if dv_cnt > 2 then
        comp_cnt = comp_cnt + 1;
      if arithm_cnt = 1000 then do;
        put skip edit('The ', arithm_cnt, 'th arithmetic number is ',
            n, ' up to which ', comp_cnt, ' are composite.')
            (a, f(4), a, f(4), a, f(3), a);
        put skip;
      end;
    end;
    n = n + 1;
  end;
end;