% Get all digits of a number
digits = iter (n: int) yields (int)
    while n > 0 do
        yield(n // 10)
        n := n / 10
    end
end digits

% Get all prime factors of a number
prime_factors = iter (n: int) yields (int)
    % Take factors of 2 out first (the compiler should optimize)
    while n // 2 = 0 do yield(2) n := n/2 end

    % Next try odd factors
    fac: int := 3
    while fac <= n do
        while n // fac = 0 do
            yield(fac)
            n := n/fac
        end
        fac := fac + 2
    end
end prime_factors

% See if a number is a Smith number
smith = proc (n: int) returns (bool)
    dsum: int := 0
    fac_dsum: int := 0

    % Find the sum of the digits
    for d: int in digits(n) do dsum := dsum + d end

    % Find the sum of the digits of all factors
    nfac: int := 0
    for fac: int in prime_factors(n) do
        nfac := nfac + 1
        for d: int in digits(fac) do fac_dsum := fac_dsum + d end
    end

    % The number is a Smith number if these two are equal,
    % and the number is not prime (has more than one factor)
    return(fac_dsum = dsum cand nfac > 1)
end smith

% Yield all Smith numbers up to a limit
smiths = iter (max: int) yields (int)
    for i: int in int$from_to(1, max-1) do
        if smith(i) then yield(i) end
    end
end smiths

% Display all Smith numbers below 10,000
start_up = proc ()
    po: stream := stream$primary_output()
    count: int := 0
    for s: int in smiths(10000) do
        stream$putright(po, int$unparse(s), 5)
        count := count + 1
        if count // 16 = 0 then stream$putl(po, "") end
    end
    stream$putl(po, "\nFound " || int$unparse(count) || " Smith numbers.")
end start_up