Project Euler solutions in Clojure

euler.clj

;;; euler.clj -- Project Euler solutions in Clojure
 
;; by Robert Campbell, http://www.robert-campbell.com/
 
;; Most of my implemented solutions are crude, utilizing only brute force. 
;; I'm more interested in exploring functional programming techniques in 
;; Clojure than in mathematical optimizations. I do, however, keep to the 
;; standard one minute rule, running all programs on a 2.66GHz Core2 Quad.

;; Problems requiring a sequence of primes were originally solved using 
;; the Fermat Test implemented below. Unfortunately this test suffers from 
;; Carmichael numbers, resulting in the occasional false positive. I've 
;; therefore switched to using clojure.contrib.seq-utils/primes even though 
;; it feels like cheating :-)
 
(ns euler.core
  (:import [java.util Calendar])
  (:use [clojure.string :only [split-lines]]
	clojure.contrib.math
	clojure.contrib.lazy-seqs
	clojure.contrib.combinatorics))

(declare square trial-division triangle-numbers pentagonal-numbers 
	 hexagonal-numbers divides? divisors fast-prime?)

; Elapsed time: 3 msecs
(defn p1 [] {:post [(= 233168 %)]}
  (reduce #(if (or (divides? 3 %2) (divides? 5 %2)) (+ %1 %2) %1) 
	  (range 1000)))

; Elapsed time: 1 msec
(defn p2 [] {:post [(= 4613732 %)]}
  (reduce #(if (even? %2) (+ %1 %2) %1) 
	  (take-while #(< % 4000000) (fibs))))

; Elapsed time: 387 msecs
(defn p3 [] {:post [(= 6857 %)]}
  (loop [i 775147 n 600851475143 times 3]
    (if (and (divides? i n) (fast-prime? i times)) i
      (recur (dec i) n times))))

; Elapsed time: 16692 msecs
(defn p4 [] {:post [(= 906609 %)]}
  (let [palindrome? (fn [x] (= (seq (str x))
			       (reverse (str x))))
	vals (reverse (range 1 1000))]
    (apply max (filter palindrome? 
		       (for [a vals b vals] (* a b))))))

; Elapsed time: 53618 msecs
(defn p5 [] {:post [(= 232792560 %)]}
  (loop [n 20]
    (if (reduce #(and (divides? %2 n) %1) true (range 11 20)) n
	(recur (+ n 20)))))

; Elapsed time: 1 msec
(defn p6 [] {:post [(= 25164150 %)]}
  (let [vals (range 1 101)]
	(- (square (reduce + vals)) 
	   (reduce + (map square vals)))))

; Elapsed time: 505 msecs
(defn p7 [] {:post [(= 104743 %)]}
  (nth primes 10001))

; Elapsed time: 25 msecs
(defn p8 [] {:post [(= 40824 %)]}
  (loop [num (map #(Integer. (str %))
		  (filter #(Character/isDigit %) (slurp "p8.dat")))
	 max 0]
    (if (empty? num) max
	(let [prod (reduce * (take 5 num))]
	  (recur (rest num) (if (> prod max) prod max))))))

; *** In Progress ***
(defn p9 []
  (for [a (iterate inc 1)
	b (iterate inc 1) :when (< a b)
	c (iterate inc 1) :when (< b c)]
    (when (and (= (+ (square a) (square b)) (square c))
	       (= (+ a b c) 1000)) [a b c])))

; Elapsed time: 22633 msecs
(defn p10 [] {:post [(= 142913828922 %)]}
  (loop [p primes sum 0]
    (if (< (first p) 2000000) 
      (recur (rest p) (+ (first p) sum)) sum)))

; Elapsed time: 54 msecs
(defn p11 [] {:post [(= 70600674 %)]}
  (let [int-row (fn [row] (map #(Integer. %) (.split row "\\s+")))
	matrix (map int-row (.split (slurp "p11.dat") "\\r\\n"))
	find-greatest-in-row (fn [row]
			       (loop [segment row max 0]
				 (if (empty? segment) max
				     (let [prod (reduce * (take 4 segment))]
				       (recur (rest segment) 
					      (if (> prod max) prod max))))))
	find-greatest (fn [matrix] (apply max (map find-greatest-in-row matrix)))
	rotate (fn [matrix] 
		 (loop [idx 0 coll ()]
		   (if (= idx (count matrix)) coll 
		       (recur (inc idx) 
			      (cons (for [row matrix] (nth row idx)) coll)))))
	diagonal (fn [matrix init-x init-y] 
		   (loop [x init-x y init-y coll ()] 
		     (if (or (>= x (count matrix)) 
			     (>= y (count matrix))) coll 
			     (recur (inc x) (inc y) 
				    (cons (nth (nth matrix x) y) coll)))))
	tilt (fn [matrix] 
	       (set (concat (map #(diagonal matrix 0 %) 
				 (range (count matrix))) 
			    (map #(diagonal matrix % 0) 
				 (range (count (first matrix)))))))]
      (max 
       (find-greatest matrix)
       (find-greatest (rotate matrix))
       (find-greatest (tilt matrix))
       (find-greatest (tilt (rotate matrix))))))    

; Elapsed time: 28513 msecs
(defn p12 [] {:post [(= 76576500 %)]}
  (loop [i 1]
    (if (> (count (divisors (nth triangle-numbers i))) 500) 
      (nth triangle-numbers i) 
      (recur (inc i)))))

; Elapsed time: 2 msecs
(defn p13 [] {:post [(= 5537376230 %)]}
  (let [vals (map #(BigInteger. %) 
		  (.split (slurp "p13.dat") "\\s+"))]
    (.substring (str (reduce + vals)) 0 10)))

; *** In Progress ***
(defn p14 []
  (letfn [(collatz [n] (loop [seq (vector n)]
			 (let [n (last seq)]
			   (cond (= n 1) seq
				 (even? n) (recur (conj seq (/ n 2)))
				 :else (recur (conj seq (inc (* n 3))))))))]
    (let [counts (apply merge (map #(array-map (count %) (first %))
				   (pmap collatz (range 1 100000))))]
      (counts (apply max (keys counts))))))
	  
; Elapsed time: 1 msec
(defn p16 [] {:post [(= 1366 %)]}
  (reduce + (map #(Integer. (str %)) 
		 (seq (str (.pow 2M 1000))))))

; Elapsed time: 76 msecs
(defn p18 
  ([] {:post [(= 1074 %)]}
     (p18 (map (fn [row] (vec (map #(Integer/valueOf %)
				   (.split row " "))))
	       (split-lines (slurp "p18.dat"))) 0 0))
  ([frontier index sum]
     (if (empty? frontier) sum
	 (let [row (first frontier)]
	   (max (p18 (rest frontier) index (+ (row index) sum))
		(if (= (inc index) (count row)) 0
		    (p18 (rest frontier) (inc index)
			 (+ (row (inc index)) sum))))))))

; Elapsed time: 219 msecs
(defn p19 [] {:post [(= 171 %)]}
  (let [calendar (doto (Calendar/getInstance) (.set 1900 12 31))
	end (doto (Calendar/getInstance) (.set 2000 12 31))]
    (loop [sundays 0]
      (if (= calendar end) sundays
	  (let [day (.get calendar Calendar/DAY_OF_WEEK)
		date (.get calendar Calendar/DATE)]
	    (.add calendar Calendar/DATE 1)
	    (recur (if (and (= day Calendar/SUNDAY) (= 1 date))
		     (inc sundays) sundays)))))))

; Elapsed time: 1100 msecs
(defn p21 [] {:post [(= 31626 %)]}
  (letfn [(sum-of-divisors [x] (reduce + (trial-division x)))]
    (loop [n 10000 sum 0]
      (let [a (sum-of-divisors n)
	    b (sum-of-divisors a)]
	(cond
	  (zero? n) sum
	  (and (= b n) (not (= a b))) (recur (dec n) (+ n sum))
	  :else (recur (dec n) sum))))))

; Elapsed time: 5052 msecs
(defn p24 [] {:post [(= "2783915460" %)]}
  (apply str (nth (permutations (range 10)) 999999)))


; --- Helper functions -----------

(defn square [x]
  (* x x))

(defn trial-division [n]
  (loop [s 2 factors #{1}]
    (cond
      (> s (ceil (sqrt n))) factors
      (divides? s n) (recur (inc s) (set (cons (/ n s) (cons s factors))))
      :else (recur (inc s) factors))))

(def triangle-numbers (map #(/ (* % (+ % 1)) 2) (iterate inc 1)))
(def pentagonal-numbers (map #(/ (* % (- (* 3 %) 1)) 2) (iterate inc 1)))
(def hexagonal-numbers (map #(* % (- (* 2 %) 1)) (iterate inc 1)))


; --- Divisor search -------------

(defn divides? [a b]
  (zero? (mod b a)))

(defn find-divisor [n test-divisor]
  (cond
    (> (square test-divisor) n) n
    (divides? test-divisor n) test-divisor
    :else (recur n (inc test-divisor))))

(defn divisors [n]
  (loop [i 1 coll []]
    (let [div (find-divisor n i)]
      (if (= n div) 
	coll
	(recur (if (= i div) (inc i) (inc div))
	       (cons div (cons (/ n div) coll)))))))

(defn smallest-divisor [n]
  (find-divisor n 2))

(defn prime? [n]
  (= n (smallest-divisor n)))


; --- Fermat Test -----------------

(defn exp-mod [base exp m]
  (cond
    (zero? exp) 1
    (even? exp) (mod (square (exp-mod base (/ exp 2) m)) m)
    :else (mod (* base (exp-mod base (dec exp) m)) m)))

(defn fermat-test 
  ([n a] (= (exp-mod a n n) a))
  ([n] (fermat-test n (inc (rand-int (dec n))))))

(defn fast-prime? [n times]
  (cond
    (zero? times) true
    (fermat-test n) (recur n (dec times))
    :else false))