smith: procedure options(main);
    /* find the digit sum of N */
    digitSum: procedure(nn) returns(fixed);
        declare (n, nn, s) fixed;
        s = 0;
        do n=nn repeat(n/10) while(n>0);
            s = s + mod(n,10);
        end;
        return(s);
    end digitSum;

    /* find and count factors of N */
    factors: procedure(nn, facs) returns(fixed);
        declare (n, nn, cnt, fac, facs(16)) fixed;
        cnt = 0;
        if nn<=1 then return(0);

        /* factors of two */
        do n=nn repeat(n/2) while(mod(n,2)=0);
            cnt = cnt + 1;
            facs(cnt) = 2;
        end;

        /* take out odd factors */
        do fac=3 repeat(fac+2) while(fac <= n);
            do n=n repeat(n/fac) while(mod(n,fac) = 0);
                cnt = cnt + 1;
                facs(cnt) = fac;
            end;
        end;

        return(cnt);
    end factors;

    /* see if a number is a Smith number */
    smith: procedure(n) returns(bit);
        declare (n, nfacs, facsum, i, facs(16)) fixed;
        nfacs = factors(n, facs);
        if nfacs <= 1 then
            return('0'b); /* primes are not Smith numbers */

        facsum = 0;
        do i=1 to nfacs;
            facsum = facsum + digitSum(facs(i));
        end;
        return(facsum = digitSum(n));
    end smith;

    /* print all Smith numbers up to 10000 */
    declare (i, cnt) fixed;
    cnt = 0;
    do i=2 to 9999;
        if smith(i) then do;
            put edit(i) (F(5));
            cnt = cnt + 1;
            if mod(cnt,16) = 0 then put skip;
        end;
    end;
    put skip list('Found', cnt, 'Smith numbers.');
end smith;