Brute-force perfect-number algorithm

Question

I whipped a reasonably fast brute-force perfect number finding function in Haskell, after seeing a similar thing in mathematica.

It works in under 3 seconds when I search up to 1000, but 10000 is too much.

isDivisor :: Integral a => a -> a -> Bool
isDivisor n a = mod n a == 0

sumDivisors :: Integral a => a -> a
sumDivisors a =
   foldr (\n -> if isDivisor a n then (+ n) else id) 0 [1..(a `div` 2)]

isPerfect :: Integral a => a -> Bool
isPerfect n = n == sumDivisors n

-- The main prime-finding function:
upUntil :: Integral a => a -> [a]
upUntil n = filter isPerfect [1..n]

It certainly isn't the fastest, or the most elegant, but I'm unsure how I can improve it.

I would like to concentrate on the isDivisor function in particular. Is there a more efficient way of checking for divisibility?

Answer 1

Based on this Wikipedia article, here's a function to calculate the first n perfect numbers:

Import Data.Numbers.Primes(primes, isPrime) -- https://hackage.haskell.org/package/primes-0.2.1.0/docs/src/Data-Numbers-Primes.html

-- A perfect number is a number of the form 2p−1× (2p − 1) 
-- where 2p − 1 is a Mersenne prime (http://en.wikipedia.org/wiki/List_of_perfect_numbers)
firstPerfectNumbers :: Int -> [Integer]
firstPerfectNumbers n = 
    take n [2 ^ (x - 1) * (2 ^ x - 1) | x <- primes, isPrime(2 ^ x - 1)]

This function returns the first 8 numbers in a second, but slows down after that.

Your original upUntil function could be implemented like this:

perfectNumbersUpTo :: Integer -> [Integer]
perfectNumbersUpTo n = 
    takeWhile (< n) [2 ^ (x - 1) * (2 ^ x - 1) | x <- primes, isPrime(2 ^ x - 1)]

But of course, seeing as there are only 35 perfect numbers with less than 1,000,000 digits and 48 known numbers, if I wanted to use a perfect number in a program I would probably use a precomputed table.

Also, if you aren't using your upUntil function to try to find new perfect numbers, or if your program isn't critically dependent on a perfect number of humongous size, then your isPerfect function would probably benefit from this:

isPerfect n | odd n = False

Because it's not known whether odd perfect numbers exist.

Answer 2

Use backticks

Backticks allow you to put the function name between the values instead of before:

isDivisor n a = mod n a == 0

becomes:

isDivisor n a = n `mod` a == 0

Give less general names

The name upUntil tells nothing about perfect numbers, in a tiny script like this this is no problem, but the habit of meningful names will play yo your advantage in larger scripts

Do not write wrong comments

Comments don't run, but you are not justified to write wrong comments

-- The main prime-finding function:

The function does not find primes.

Answer 3

I'd say it's hard to improve isDivisor, as it's just a call to the primitive function mod. However, there are other areas for improvement.

In particular, going through [1..(adiv2)] is a very inefficient method for listing divisors. A much more efficient method would be to factorize a and then compute all its divisors from that (see Divisor function, in particular the formula specialized for σ₁(n)).

So my suggestion would be:

• Generate all primes up to your upper bound.
• Use them to factorize each number in the range.
• Use the formula for σ₁(n) to compute the sum of divisors.

See also package arithmoi which implements a lot of such functions.

Answer 4

The Answer based on mersenne primes is unbeatable in speed, but if you want to extend your code for abundant numbers sumDivisors n > n or deficient numbers ... < n, there are some options left.

Low-Level optimizations

Since you are asking for speed improvements, here is your code with

main = print $ upUntil 10000

$ time ./Main
[6,28,496,8128]

real 0m3.193s
user 0m0.000s
sys 0m0.046s

Just by using Int like this, time drops to 0.920s.

main = print $ upUntil (10000 :: Int)

mod handles negative values nicely, here rem sufficient and yields another drop to 0.710s:

isDivisor n a = rem n a == 0

Changing all functions' type signatures to Int does not provide any boost - GHC is already specializing the functions.

Algorithmic

You can gain a lot by taking also b = n / a if n / a is integer like this, where you just need to iterate up to the square root:

sumDivisors a =
   foldr (\n -> let (q,r) = a `quotRem` n in
       if r==0 then (+ (n+q)) else id) 1
       [2..(floor . sqrt . fromIntegral $ a)]

(This does not handle perfect squares properly)