[Haskell-cafe] Zumkeller numbers

[Daniel Fischer]

Am Montag 07 Dezember 2009 22:33:32 schrieb Frank Buss:
> Anyone interested in writing some lines of Haskell code for generating the
> Zumkeller numbers?
>
> http://www.luschny.de/math/seq/ZumkellerNumbers.html
>
> My C/C# solutions looks clumsy (but is fast). I think this can be done much
> more elegant in Haskell with lazy evaluation.

A fairly naive but not too slow version:
---------------------------------------------------------------------------------
module Main (main) where

import Data.List (sort)
import Control.Monad (liftM2)
import System.Environment (getArgs)

main = do
    args <- getArgs
    let bd = case args of
                (a:_) -> read a
                _ -> 1000
    mapM_ print $ filter isZumkeller [2 .. bd]

trialDivPrime :: [Int] -> Int -> Bool
trialDivPrime (p:ps) n = (p*p > n) || (n `mod` p /= 0) && trialDivPrime ps n

oprs = 5:7:filter (trialDivPrime oprs) (scanl (+) 11 $ cycle [2,4])

primes :: [Int]
primes = 2:3:oprs

factors :: Int -> [(Int,Int)]
factors n
    | n < 0     = factors (-n)
    | n == 0    = [(0,1)]
    | otherwise = go n primes
      where
        go 1 _ = []
        go m (p:ps)
            | p*p > m   = [(m,1)]
            | otherwise =
              case m `quotRem` p of
                (q,0) -> case countFactors q p of
                            (m',k) -> (p,k+1):go m' ps
                _ -> go m ps

countFactors :: Int -> Int -> (Int,Int)
countFactors n p = go n 0
      where
        go m k = case m `quotRem` p of
                    (q,0) -> go q (k+1)
                    _ -> (m,k)

divisors :: Int -> [Int]
divisors n = sort $ foldr (liftM2 (*)) [1] ppds
      where
        ppds = map (\(p,k) -> take (k+1) $ iterate (*p) 1) $ factors n

partition :: Int -> [Int] -> [[Int]]
partition 0 _   = [[]]
partition n []  = []
partition n (p:ps)
    | n < p     = partition n ps
    | otherwise = [p:parts | parts <- partition (n-p) ps] ++ partition n ps

isZumkeller :: Int -> Bool
isZumkeller n
    | sigma < 2*n   = False
    | odd sigma     = False
    | otherwise     = not . null $ partition target candidates
      where
        divs = divisors n
        sigma = sum divs
        target = (sigma-n) `quot` 2
        candidates = reverse $ takeWhile (<= target) divs
----------------------------------------------------------------------------------------------------

Of course, with sieving, you'd be much faster.

>
> Not related to Haskell, but do you think semi-Zumkeller numbers are
> semi-perfect numbers?

The site you linked to says so. I've not investigated.

> Maybe some Haskell code for testing it for some numbers?