# Brute-force perfect-number algorithm

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?

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.

## 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.

• Whoops, I just noticed that comment! Thanks!
– AJF
Apr 13, 2015 at 6:32

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.

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)