import Data.List (sort, group, union)
import Data.Array

pairwise :: (a -> a -> a) -> [a] -> [a]
pairwise f (xs:ys:t) = f xs ys : pairwise f t
pairwise _ t = t

primes :: [Int]
primes = 2 : _Y ((3 :) . gaps 5 . _U . map (\p-> [p*p, p*p+2*p..]))
    where
        _Y g = g (_Y g)                      -- recursion, Y combinator
        _U ((x:xs):t) = x : (union xs . _U . pairwise union) t   -- ~= nub.sort.concat
        gaps k s@(x:xs) 
            | k < x     = k : gaps (k+2) s    -- ~= [k,k+2..]\\s, when
            | otherwise =     gaps (k+2) xs   --  k <= head s && null(s\\[k,k+2..])

factorize :: Int -> [Int]
factorize n = primeFactors n primes where
    primeFactors 1 _ = []
    primeFactors _ [] = []
    primeFactors m (p:ps) | m < p * p = [m]
                          | r == 0 = p : primeFactors q (p:ps)
                          | otherwise = primeFactors m ps
                          where (q, r) = quotRem m p

primePowers :: Int -> [(Int, Int)]
primePowers n = [(head x, length x) | x <- group $ factorize n]

divisors :: Int -> [Int]
divisors n = filter (<n) $ map product $ sequence
    [take (k+1) $ iterate (p*) 1 | (p, k) <- primePowers n]

upperBound :: Int
upperBound = 20161

abundant :: Int -> Bool
abundant n = (sum . divisors) n > n

abundantsArray :: Array Int Bool
abundantsArray = listArray (1, upperBound) $ map abundant [1..upperBound]

abundants :: [Int]
abundants = filter (abundantsArray !) [1..upperBound]

remainders :: Int -> [Int]
remainders x = map (x-) $ takeWhile (<= x `quot` 2) abundants

sums :: [Int]
sums = filter (any (abundantsArray !) . remainders) [1..upperBound]

main :: IO ()
main = print $ sum [1..upperBound] - sum sums

$ ghc -O2 -o abundant abundant.hs
$ time ./abundant
real   0m0.034s
user   0m0.034s
sys    0m0.000s

#!/usr/bin/env python3
import math
from collections import defaultdict
from itertools import *
from functools import reduce
def factorize(n):
    if n < 1:
        raise ValueError('fact() argument should be >= 1')
    if n == 1:
        return []  # special case
    res = []
    # iterate over all even numbers first.
    while n % 2 == 0:
        res.append(2)
        n //= 2
    # try odd numbers up to sqrt(n)
    limit = math.sqrt(n+1)
    i = 3
    while i <= limit:
        if n % i == 0:
            res.append(i)
            n //= i
            limit = math.sqrt(n+i)
        else:
            i += 2
    if n != 1:
        res.append(n)
    factors = sorted(res)
    histogram = defaultdict(int)
    for factor in factors:
        histogram[factor] += 1

    return list(histogram.items())

def divisors(n):
    factors = factorize(n)
    nfactors = len(factors)
    f = [0] * nfactors
    while True:
        yield reduce(lambda x, y: x*y, [factors[x][0]**f[x] for x in range(nfactors)], 1)
        i = 0
        while True:
            f[i] += 1
            if f[i] <= factors[i][1]:
                break
            f[i] = 0
            i += 1
            if i >= nfactors:
                return

def proper_divisors(n):
    return list(divisors(n))[:-1]

def classify(n):
    total = sum(proper_divisors(n))
    if total == n:
        # perfect
        return 0
    elif total > n:
        # abundant
        return 1
    else:
        # deficient
        return -1

def main():
    abundant = set(number for number in range(2, 30000) if classify(number) == 1)
    sums = sorted(set(sum(c) for c in combinations_with_replacement(abundant, 2)))
    print((sum(number for number in range(1,30000) if number not in sums)))
    
if __name__ == "__main__":
    main()

$ time python3 abundant.py
real   0m16.594s
user   0m16.593s
sys    0m0.000s

#!/usr/bin/env ruby
require 'mathn' 

class Integer 
  def divisors
    return [1] if self == 1
    primes, powers = self.prime_division.transpose 
    exponents = powers.map{|i| (0..i).to_a} 
    divisors = exponents.shift.product(*exponents).map do |powers| 
      primes.zip(powers).map{|prime, power| prime ** power}.inject(:*) 
    end 
    divisors.take divisors.length - 1
  end

  def abundant?
    self.divisors.reduce(:+) > self
  end
end

abundants = (1..28213).select { |n| n.abundant? }
i = 0
sums = []
abundants.each do |x|
  abundants[i..abundants.length].each do |y|
    sum = x + y
    sums << sum unless sum > 28213
  end
  i += 1
end
sums.uniq!
puts (1..28213).reject { |n| sums.include? n }.reduce(:+)

$ time ruby abundant.rb
real   0m4.263s
user   0m4.231s
sys    0m0.032s

fn sum_divisors(n: usize) -> usize {
    let mut result = 0;
    let max = 1 + (n as f64).sqrt() as usize;
    for i in 2..max {
        if n % i == 0 {
            let x = n / i;
            if x == i {
                result += i;
            } else {
                result += i + x;
            }
        }
    }
    1 + result
}

fn main() {
    let max = 28123;
    let abundant: Vec<usize> = (2..max + 1).filter(|&n| sum_divisors(n) > n).collect();
    let mut abundant_sums = vec![false; 2 * max + 1];
    for i in 0..abundant.len() {
        for j in i..abundant.len() {
            abundant_sums[abundant[i] + abundant[j]] = true;
        }
    }
    let sum: usize = (1..max + 1).filter(|&i| !abundant_sums[i]).sum();
    println!("{}", sum);
}

$ rustc -C target-cpu=native -C opt-level=3 -o abundant abundant.rs
$ time ./abundant
real   0m0.030s
user   0m0.030s
sys    0m0.000s