get "libhdr"

// Find the sum of the digits of N
let digsum(n) =
    n<10 -> n,
    n rem 10 + digsum(n/10)

// Factorize N
let factors(n, facs) = valof
$(  let count = 0 and fac = 3

    // Powers of 2
    while n>0 & (n & 1)=0
    $(  n := n >> 1
        facs!count := 2
        count := count + 1
    $)

    // Odd factors
    while fac <= n
    $(  while n rem fac=0
        $(  n := n / fac
            facs!count := fac
            count := count + 1
        $)
        fac := fac + 2
    $)

    resultis count
$)

// Is N a Smith number?
let smith(n) = valof
$(  let facs = vec 32
    let nfacs = factors(n, facs)
    let facsum = 0
    if nfacs<=1 resultis false   // primes are not Smith numbers
    for fac = 0 to nfacs-1 do
        facsum := facsum + digsum(facs!fac)
    resultis digsum(n) = facsum
$)

// Count and print Smith numbers below 10,000
let start() be
$(  let count = 0
    for i = 2 to 9999 if smith(i)
    $(  writed(i, 5)
        count := count + 1
        if count rem 16 = 0 then wrch('*N')
    $)
    writef("*NFound %N Smith numbers.*N", count)
$)