Solution to NWERC 2000, Problem 4: Smith Numbers.

> import ReaderMonad
> import Int
> import ST
> import MArray
> type Number = Int32

List of prime numbers (and their sqaures) computed by the sieve of
Eratosthenes.  To make it as fast as possible, we use a mutable array
in a state thread monad.

> bound :: Number
> bound = 32000

> primes :: [(Number,Number)]
> primes = runST statesieve

> statesieve :: ST s [(Number,Number)]
> statesieve = do sieve <- marray (2,bound) :: ST s (STArray s Number Bool)
>                 sieve // zip [2..bound] (repeat True)
>                 let select n = do b <- get sieve n
>                                   if b then (sieve // comps n)
>                                      else return ()
>                 mapM_ select [2..ceiling (sqrt (fromIntegral bound))]
>                 l <- assocs sieve
>                 return [ (x,x*x) | (x,True) <- l ]
>   where comps x = tail [ (y,False) | y <- [x,x+x..bound] ]

Factorize a number by trial division.

> factor :: Number -> [Number]
> factor n = f n primes
>   where f n l@((p,p2):ps) | n < p2 = [n]   -- "n<=p" is too slow; I've tried.
>                           | r == 0 = p : f q l
>                           | otherwise = f n ps
>           where (q,r) = divMod n p
>         f n [] = [n]

Just doing the above in C/C++/Java/Pascal will take twice as much code
and forever to debug.

The decimal digits of a positive number, least significant ones first.

> digits :: Number -> [Number]
> digits n | n < 10 = [n]
>          | otherwise = r : digits q
>   where (q,r) = divMod n 10

To determine if a number is a Smith number:

> isSmith :: Number -> Bool
> isSmith n = let fs = factor n
>                 in length fs > 1 &&
>                    sum (digits n) == sum (concatMap digits (factor n))

Even the main program is short and sweet.

> main :: IO ()
> main = do n <- readLn
>           sequence_ (replicate n (do i <- readLn
>                                      print (head (filter isSmith [i+1..]))))