Divisors

February 14, 2012

We sieve the list of divisors, their count and their sum all in the same double-nested loops:

(define (divisors n)
  (let ((divs (make-vector (+ n 1) (list)))
        (cnts (make-vector (+ n 1) 0))
        (sums (make-vector (+ n 1) 0)))
    (do ((i 1 (+ i 1))) ((<= n i) (values divs cnts sums))
      (do ((j i (+ j i))) ((<= n j))
        (vector-set! divs j (cons i (vector-ref divs j)))
        (vector-set! cnts j (+ 1 (vector-ref cnts j)))
        (vector-set! sums j (+ i (vector-ref sums j)))))))

We will store the lists of divisors, counts and sums in global variables that are updated to the maximum size desired:

(define divs #f)
(define cnts #f)
(define sums #f)
(define max-d 0)

(define (update n)
  (call-with-values
    (lambda () (divisors n))
    (lambda (d c s)
      (set! divs d)
      (set! cnts c)
      (set! sums s)))
  (set! max-d n))

Thus, (vector-ref divs 360) is a list of the divisors of 360, (vector-ref cnts 360) is their count, and (vector-ref sums 360) is their sum:

> (reverse (vector-ref divs 360))
(1 2 3 4 5 6 8 9 10 12 15 18 20 24 30 36 40 45 60 72 90 120 180 360)
> (vector-ref cnts 360)
24
> (vector-ref sums 360)
1170

Note that the sum includes the number itself. Often, what is wanted is the sum excluding the number itself, which is computed as (- (vector-ref sums n) n).

Perfect numbers, which were studied by the ancient Greek mathematician Pythagoras circa 500 B.C., are those with the sum of the divisors, less the number itself, equal to the number itself; here we compute the perfect numbers less than n:

(define (perfect n)
  (define (s n) (- (vector-ref sums n) n))
  (when (< max-d n) (update n))
  (let ((ps (list)))
    (do ((p 1 (+ p 1)))
        ((< n p) (reverse ps))
      (when (= p (s p))
        (set! ps (cons p ps))))))

> (perfect 10000)
(6 28 496 8128)

Amicable pairs are pairs of numbers that are each the sum of the other’s divisors (excluding the numbers themselves). For instance, the sum of the divisors of 220 is 1 + 2 + 4 + 5 + 10 + 11 + 20 + 22 + 44 + 55 + 110 = 284, and the sum of the divisors of 284 is 1 + 2 + 4 + 71 + 142 = 220, so the numbers 220 and 284 form an amicable pair. Here we compute the amicable pairs less than n, showing only the smaller of the pair:

(define (amicable n)
  (define (s n) (- (vector-ref sums n) n))
  (when (< max-d (* 5 n)) (update (* 5 n)))
  (let loop ((a 1) (as (list)))
    (cond ((< n a) (reverse as))
          ((and (= a (s (s a))) (< a (s a)))
            (loop (+ a 1) (cons a as)))
          (else (loop (+ a 1) as)))))

> (display (amicable 10000))
(220 1184 2620 5020 6232)

We multiply by 5 because sometimes the other half of an amicable pair is larger than the starting half, and 5 seems to give sufficient slop, as determined by experiment.

You can run the program at http://programmingpraxis.codepad.org/3y4c9Ybl. And if you wish, you can now easily solve Problem 21 at Project Euler.

About these ads

Pages: 1 2

11 Responses to “Divisors”

  1. Paul G. said

    First time posting something to this website. Here’s a perl solution to the problem – http://pastebin.com/4gaWFTzJ

  2. Abe C. said

    Paul,

    If I’m not mistaken, it doesn’t look like your perl solution is actually seiving in order to determine the divisors of each number, which was the point of the exercise. It looks instead like you are performing modular arithmetic on each number $n and checking for a remainder against each possible divisor, $_. Am I mistaken?

    Here’s a java solution: http://pastebin.com/6xL9BwS0

    Thanks.

  3. Jussi Piitulainen said

    Comprehensions :) With multiple filters. (Ignore the element at index 0.)

    def divisors(m):
        ds = [ [] for k in range(m) ]
        for k in range(1, m):
            for n in range(k, m, k):
                ds[n].append(k)
        return ds
    
    def divisorsums(m):
        ss = [ 0 for k in range(m) ]
        for k in range(1, m):
            for n in range(k, m, k):
                ss[n] += k
        return ss
    
    def perfect(sums):
        return [ n for n, s in enumerate(sums)
                 if n > 0
                 if n + n == s ]
    
    def amicable(sums):
        return [ (n, s - n) for n, s in enumerate(sums)
                 if s - n < len(sums)
                 if n < s - n
                 if n == sums[s - n] - (s - n) ]
    

    The divisor lists are redundant. Their sums do the work. In amicable, the condition to exclude double numbers and duplicates also excludes the double zero. In perfect, the 0 is excluded by its very own condition.

    >>> perfect(divisorsums(10000))
    [6, 28, 496, 8128]
    >>> amicable(divisorsums(10000))
    [(220, 284), (1184, 1210), (2620, 2924), (5020, 5564), (6232, 6368)]
    
  4. Axio said
    (defun make-divisors-sieve (n)
      (let ((a (make-array (1+ n) :initial-element nil)))
        (loop for i from 1 to n
              do (loop for j from 1 to (truncate n i)
                    do (push i (aref a (* i j)))))
        a))
    
  5. ardnew said

    This site seems to have attracted mainly a crowd of functional programmers with crazy powerful language facilities at their disposal.

    We need some more love given to the gritty implementation details of these solutions (particularly when the posed problem is simple).

    But enough trying to justify my approach..

    No fancy data structures, no analytical computations.. purely imperative number crunching:

    #include <stdio.h>
    #include <stdlib.h>
    #include <string.h>
    
    // linked list storing divisors of an integer
    struct divisor_node
    {
      int divisor;
      struct divisor_node *next;
    };
    typedef struct divisor_node divisor;
    
    // command-line usage
    void usage(char *name)
    {
      printf("\nusage:\n\t%s <max>\n\n", name);
    }
    
    // allocate and initialize memory for our table
    divisor **alloc_table(int n)
    {
      divisor **table;
      int b = n * sizeof table;
    
      memset(table = malloc(b), 0, b);
      return table;
    }
    
    // sieve the divisors of all n elements
    void build_table(divisor **table, int n)
    {
      int i, j;
      divisor *curr;
    
      for (i = 0; i < n; ++i)
      {
        for (j = i + 1; j <= n; j += i + 1)
        {
          curr = (divisor *) malloc(sizeof(divisor));
          curr->divisor = i + 1;
          curr->next = table[j - 1];
          table[j - 1] = curr;
        }
      }
    }
    
    // print the linked list of all n table entries
    void print_table(divisor **table, int n)
    {
      int i;
      divisor *curr;
    
      for (i = 0; i < n; ++i)
      {
        printf("%d: ", i + 1);
        curr = table[i];
    
        while (curr)
        {
          printf("%d", curr->divisor);
          printf("%s",
            (curr = curr->next) ? ", " : "\n");
        }
      }
    }
    
    // free the memory we allocated for our table
    void free_table(divisor **table, int n)
    {
      divisor *next;
    
      while (n--) // each linked list
      {
        while (table[n]) // each element
        {
          next = table[n]->next;
          free(table[n]); table[n] = next;
        }
      }
      
      free(table);
    }
    
    // sum all children of (and including) node
    int list_sum(divisor *node)
    {
      int sum = 0;
      
      while (node)
      {
        sum += node->divisor;
        node = node->next;
      }
      
      return sum;
    }
    
    // a perfect number is a positive integer that is 
    // equal to the sum of its proper positive divisors
    void print_perfect(divisor **table, int n)
    {
      while (n--)
        if (list_sum(table[n]->next) == n + 1)
        {
          printf("%d\n", n + 1);
        }
    }
    
    // amicable numbers are two different numbers so 
    // related that the sum of the proper divisors of 
    // each is equal to the other number.
    void print_amicable(divisor **table, int n)
    {
      int s;
      
      while (n--)
        if ((s = list_sum(table[n]->next)) && s <= n &&
            n + 1 == list_sum(table[s - 1]->next))
        {
          printf("%d, %d\n", n + 1, s);
        }
    }
    
    int main(int argc, char *argv[])
    {
      int n;
      divisor **table;
    
      // validate user input
      if (argc > 1 && (n = atoi(argv[1])) > 0)
      {
        table = alloc_table(n);   
        build_table(table, n);
    
        //printf("table of divisors:\n");
        //print_table(table, n);
    
        printf("\namicable pairs:\n");
        print_amicable(table, n);
    
        printf("\nperfect numbers:\n");
        print_perfect(table, n);
    
        free_table(table, n);
    
        return 0;
      }
      // handle invalid command line parameters
      else
      {
        usage(*argv);
    
        return 1;
      }
    }
    
    

    I would be interested in seeing some of the functional solutions’ performance!

    $ time ./divisors.exe 10000
    
    amicable pairs:
    6368, 6232
    5564, 5020
    2924, 2620
    1210, 1184
    284, 220
    
    perfect numbers:
    8128
    496
    28
    6
    
    real	0m0.104s
    user	0m0.109s
    sys	0m0.015s
    
  6. programmingpraxis said

    Ardnew: I added timing code at http://programmingpraxis.codepad.org/WlLh3Gfx. Time was 10ms to calculate the perfect numbers and 90ms to calculate the amicable pairs, giving a total 100ms which barely beats your 104ms. I ran out of fingers before I finished counting your lines of code, which might say something about functional programmers with crazy powerful language facilities at their disposal.

  7. ardnew said

    haha, you’ll need all the fingers and toes of everyone on the block to count my lines of code

  8. ardnew said

    And I would like to point out the timing result I posted is for process invocation, building the divisors table, finding both perfects and amicable pairs, and then printing to screen.

    Timing only the amicable pair and perfect number calls (as your scheme implementation does) results in significantly less time:

    amicable pairs:
    6368, 6232
    5564, 5020
    2924, 2620
    1210, 1184
    284, 220
    time elapsed = 16.0 ms
    
    perfect numbers:
    8128
    496
    28
    6
    time elapsed = 15.0 ms
    
  9. programmingpraxis said

    Actually, both perfect and amicable call update, so the divisor table is built twice in the timings I gave. A version that extracts the calculation of the divisors table is given at http://programmingpraxis.codepad.org/TLMWMAkg. The time to calculate 50000 entries in the divisors tables is 80ms, and the time to compute the perfect numbers and amicable pairs is 0.

  10. ardnew said

    Very nice, cache hits always help runtime

  11. Manoj Senguttuvan said

    And, that’s it! :)

    <? $n=$argv[1];
    for($count=0,$j=0,$i=1;$i<=$n;$i++,$j=0,$count=0)
                    while($j++<=$i)
                            $arr[$i][$count]=$i%$j==0?$j*(++$count/$count):$arr[$i][$count];
    echo "Perfect: ";
    for($sums=array(),$sum=0,$i=1;$i<=$n;$i++,$sum=0)
    {
            foreach( $arr[$i] as $j => $val)
                    $sum+=$val;
             echo $sum==2*$i?$i." ":"";
            $sums[$i]=$sum-$i;
    }
    print "\nAmicable Pairs: ";
    for($i=2;$i<=$n;$i++)
    {
            $x=$sums[$i];
            if($x<$i && $sums[$x]==$i)
                    echo "($i,$x) ";
    } ?>
    

    Output:

    Perfect: 6 28 496 8128
    Amicable Pairs: (284,220) (1210,1184) (2924,2620) (5564,5020) (6368,6232)
    

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 609 other followers

%d bloggers like this: