Skip to content

Instantly share code, notes, and snippets.

@nvanderw
Created November 20, 2012 05:13
Show Gist options
  • Save nvanderw/4116173 to your computer and use it in GitHub Desktop.
Save nvanderw/4116173 to your computer and use it in GitHub Desktop.
Project Euler #23
#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))))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment