A little scheme programming challenge

Question

I am learning scheme to get something new from a programming language and this code below is the solution to Project Euler question 21 but the code runs 10x slower than the listed Python code when I use Chibi-Scheme. Can any Scheme craftsman refactor the code into the Scheme idiom for efficiency.

Scheme Code:

(define (sum-of-amicable-pairs n)
  (let ((sums (make-vector n)))
        (for-each
          (lambda (i) 
            (vector-set! sums i
                         (reduce + 0 
                                 (filter (lambda (j) (= (remainder i j) 0)) 
                                         (iota (+ 1 (quotient i 2)) 1 1))))) 
            (iota n 0 1))

        (let loop ((res (make-vector n 0))
                   (i 0))
          (cond
            ((= i n) (reduce + 0 (vector->list res)))
            ((and (< (vector-ref sums i) n) (not (= (vector-ref sums i) i)) 
                  (= (vector-ref sums (vector-ref sums i)) i))
             (begin
               (vector-set! res i i)
               (vector-set! res (vector-ref sums i) (vector-ref sums i))
               (loop res (+ 1 i))))
            (else
              (loop res (+ i 1)))))))

(display (sum-of-amicable-pairs 10000))
(newline)

Python code:

def amicable_pairs(n):
    """returns sum of all amicable pairs under n. See project euler for 
    definition of an amicable pair"""
    div_sum = [None]*n
    amicable_pairs_set = set()
    for i in range(n):
        div_sum[i] = sum([j for j in range(1, i/2 + 1) if i%j == 0])
    for j in range(n):
        if div_sum[j] < n and div_sum[div_sum[j]] == j and div_sum[j] != j:
            amicable_pairs_set.add(j)
            amicable_pairs_set.add(div_sum[j])
    #print sum(amicable_pairs_set)
    return sum(list(amicable_pairs_set))

Answer 1

I would optimize the python code first and then translate it to the scheme (sorry, can't find Scheme compiler for Windows). Here is python solution that 100-200 times faster (the idea based on the Sieve of Eratosthenes algorithm):

python25\python 21.py
amicable_pairs      : 31626  Time (s): 2.313000
amicable_pairs_fast : 31626  Time (s): 0.015000

python26\python 21.py
amicable_pairs      : 31626  Time (s): 1.969000
amicable_pairs_fast : 31626  Time (s): 0.016000

python27\python 21.py
amicable_pairs      : 31626  Time (s): 1.782000
amicable_pairs_fast : 31626  Time (s): 0.000000

python32\python 21.py
amicable_pairs      : 31626  Time (s): 3.157000
amicable_pairs_fast : 31626  Time (s): 0.015000

pypy\pypy 21.py
amicable_pairs      : 31626  Time (s): 0.157000
amicable_pairs_fast : 31626  Time (s): 0.015000

.

 import time

"""

    http://projecteuler.net/problem=21
    Let d(n) be defined as the sum of proper divisors of n (numbers less than n which divide evenly into n).
    If d(a) = b and d(b) = a, where a !=  b, then a and b are an amicable pair and each of a and b are called amicable numbers.

    For example, the proper divisors of 220 are 1, 2, 4, 5, 10, 11, 20, 22, 44, 55 and 110; therefore d(220) = 284.
    The proper divisors of 284 are 1, 2, 4, 71 and 142; so d(284) = 220.
    Evaluate the sum of all the amicable numbers under 10000.

    http://codereview.stackexchange.com/questions/13476/a-little-scheme-programming-challenge

"""

def amicable_pairs(n):
    """returns sum of all amicable pairs under n. See project euler for
    definition of an amicable pair"""
    div_sum = [None]*n
    amicable_pairs_set = set()
    for i in range(n):
        div_sum[i] = sum([j for j in range(1, int(i/2) + 1) if i%j == 0])
    for j in range(n):
        if div_sum[j] < n and div_sum[div_sum[j]] == j and div_sum[j] != j:
            amicable_pairs_set.add(j)
            amicable_pairs_set.add(div_sum[j])
    #print sum(amicable_pairs_set)
    return sum(list(amicable_pairs_set))

def amicable_pairs_fast(n):

    div_sum = [1]*n
    for i in range(2, int(n/2) + 1):
        for j in range(i, n, i):
            if j != i:
                div_sum[j] += i
    total = 0
    for i in range(1, n):
        s = div_sum[i]
        if s < n:
           s1 = div_sum[s]
           if  s1 == i and s != i:
              total += s1
    return total

t = time.time()
print("amicable_pairs      : %d " % amicable_pairs(10000) + " Time (s): %f " % (time.time()-t))

t1 = time.time()
print("amicable_pairs_fast : %d " % amicable_pairs_fast(10000) + " Time (s): %f " % (time.time()-t1))

Answer 2

Well, for one thing, the python code is using a set for its amicable_pair_set, where-as you're using a large vector and setting the n'th element to n when you need to add to the set. It's a reasonable way to imitate a set, if your scheme doesn't have a set library; however, this situation doesn't need true set semantics. You can use a simple list instead, so your named let becomes:

(let loop ((res '())
           (i 0))
      (cond
       ((= i n) (reduce + 0 res))
       ((and (< (vector-ref sums i) n) (not (= (vector-ref sums i) i)) 
         (= (vector-ref sums (vector-ref sums i)) i))
    (loop (cons i res) (+ 1 i)))
       (else
    (loop res (+ i 1)))))

This keeps it close to the python code, but you could go a step further and have res be an integer accumulator, adding to it instead of consing to it, getting rid of the need to add the list at the end. Of course, the same optimization could be made to the python version as well.

Unfortunately, I've optimized the part of the code that runs quick :-). The real work is being done above, when the sums vector is being populated. Since that is a pretty direct translation, I would chalk it up to chibi's implementation of scheme versus whatever implementation of python you're using (PyPy, for instance, is probably going to be faster). I use racket-scheme which is jit compiled. Calling your code returns an answer in ~680ms for me. Calling it with the default python 2.7.3 compiler gives me a result in ~1400ms.

Here's a version which only calculates the sum of its proper divisors for the numbers that are needed, and then stores them for future use. Basically memoizing the results. The first run is slightly faster for me ~30ms, and then subsequent runs return in less than 1ms.

(define d
  (let ((cache (make-vector 10000 #f)))
    (lambda (n)
      (or (vector-ref cache n)
      (let ((val (reduce + 0 (filter (lambda (d) (= 0 (remainder n d)))
                     (iota (quotient n 2) 1)))))
        (vector-set! cache n val)
        val)))))

(define (sum-of-amicable-pairs n)
  (let loop ((a 0)
         (sum 0))
    (if (< a n)
      (let ((b (d a)))
        (loop (+ a 1)
          (if (and (< a b n) (= (d b) a))
              (+ sum a b)
              sum)))
      sum)))

(time (sum-of-amicable-pairs 10000))