Project Euler #23

p23.rkt

#lang racket

; Solution to Project Euler problem #23
 
(require racket/set)

; Predicate to check if m | n
(define (divides? m n)
  (= 0 (modulo n m)))

(define (list-split lst null-case pair-case)
  (if (null? lst)
      (null-case)
      (pair-case (car lst) (cdr lst))))

; Get a suitable list of prime numbers using the Sieve of Eratosthenes
(define (sieve n)
  (letrec ((vec-len (- n 1))

           ; We can stop once we finish with the sqrt(n)th element
           (max-search (inexact->exact
                         (ceiling
                           (sqrt n))))

           (vec (make-vector vec-len))

           ; Loop to initialize the vector, setting the ith element to 2 + i,
           ; so that the vector contains [2,3,...,n]. Should be passed the
           ; starting value for i, which is in this case 0.
           (init-vec
             (lambda (i)
               (if (< i vec-len)
                   (begin
                     (vector-set! vec i
                               (+ i 2))
                     (init-vec (+ i 1)))
                   (void))))

           ; We use two loops here, the "prime loop," which is the outer loop,
           ; and the "composite loop," which is the inner loop. The prime loop
           ; starts at the first element, 2, and calls the composite loop to
           ; "cross off" all elements [4, 6, ..]. We eliminate these elements
           ; by setting them to 0.
           
           (composite-loop
             (lambda (i stride)
               (if (< i vec-len)
                   (begin
                     (vector-set! vec i 0)
                     (composite-loop (+ i stride) stride))
                   (void))))

           (prime-loop
             (lambda (i)
               (if (>= i max-search)
                   (void) ; We're done sieving

                   ; The element we're looking at must be either 0 or a prime
                   ; number. If it's 0, we skip it. Otherwise, we use
                   ; composite-loop to strike off numbers which are multiples
                   ; of this prime.
                   (let ((current-prime
                           (vector-ref vec i)))
                     (begin
                       (if (= current-prime 0)
                           (void)
                           (composite-loop (+ i current-prime) current-prime))
                       (prime-loop (+ i 1)))))))

           ; When we're done, we use this loop to iterate over the vector,
           ; skipping zeroes and making a cons list of the primes.
           (collect
             (lambda (i)
               (if (>= i vec-len)
                   null
                   (let ((head (vector-ref vec i))
                         (tail (collect
                                 (+ i 1))))
                     (if (not (= 0 head))
                         (cons head tail)
                         tail))))))
                         

    (init-vec 0)
    (prime-loop 0)
    (collect 0)))

(define factor-base (sieve 30000))

; Given an ascending list of prime numbers, try to factorize n.
(define (factorize prime-base n)
  (if (= n 1)
      null
      (list-split
        prime-base
        (lambda ()
          n)
        (lambda (head tail)
          (if (divides? head n)
              (cons head
                    (factorize prime-base (/ n head)))
              (factorize tail n))))))

; Given an ascending factorization like '(2 2 3 5 5), make a list of unique
; primes raised to powers, like '((2 2) (3 1) (5 2))
(define (group-powers factorization)
  (foldr
    (lambda (elem acc)
      (list-split
        acc
        (lambda () ; Null case
          (list
            (list elem 1)))
        (lambda (head tail)
          (if (= elem (first head))
              (cons
                (list (first head) (+ 1 (second head)))
                tail)
              (cons
                (list elem 1)
                acc)))))
    null
    factorization))

; Sigma function of an integer. This is a multiplicative function, so that
; sigma(m*n) = sigma(m) * sigma(n) if m and n are coprime.
;
; Hence we factor a number into its prime power components, find their sigma
; functions, and then multiply those.
(define (sigma n)
  (foldr * 1
    (map
      (lambda (pair)
        (let ((p (first pair))
              (a (second pair)))
          (/
            (- 
              (expt p
                    (+ a 1))
              1)
            (- p 1))))
      (group-powers (factorize factor-base n)))))

; Define an abundant number to be one where sigma(n) > 2 n
(define (abundant? n)
        (> (sigma n)
           (* 2 n)))

; Creates ranges.
(define (seq start step stop)
  (if (< start stop)
      (cons start
            (seq (+ start step) step stop))
      null))

; Initialize abundant-sums to be an empty set. We will use this to record the
; sums of two abundant numbers.
(define abundant-sums (set))

; Get an ascending list of abundant numbers
(define abundant-numbers
  (filter
    abundant?
    (seq 12 1 28123)))

; Iterate through pairs of these numbers, adding them and writing them to the set.
(for-each
  (lambda (x)
    (let/cc loop-break
      (for-each
        (lambda (y)
          (set!
            abundant-sums
            (let ((sum (+ x y)))
              (if (or
                    (> sum 28123)
                    (> y x))
                  (loop-break (void))
                  (set-add abundant-sums (+ x y))))))
        abundant-numbers)))
  abundant-numbers)

(displayln 
  (foldl + 0
    (filter
      (lambda (elem)
        (not (set-member? abundant-sums elem)))
      (seq 1 1 28123))))