--- a/src/euler.lisp Tue Aug 08 16:57:14 2017 -0400
+++ b/src/euler.lisp Wed Aug 09 14:31:30 2017 -0400
@@ -1714,6 +1714,71 @@
(when reversible
(sum (if (= i j) 1 2)))))))
+(defun problem-357 ()
+ ;; Consider the divisors of 30: 1,2,3,5,6,10,15,30. It can be seen that for
+ ;; every divisor d of 30, d+30/d is prime.
+ ;;
+ ;; Find the sum of all positive integers n not exceeding 100 000 000 such that
+ ;; for every divisor d of n, d+n/d is prime.
+ (labels ((check-divisor (n d)
+ (primep (+ d (truncate n d))))
+ (prime-generating-integer-p (n)
+ (declare (optimize speed)
+ (type fixnum n)
+ (inline divisors-up-to-square-root))
+ (every (curry #'check-divisor n)
+ (divisors-up-to-square-root n))))
+ ;; Observations about the candidate numbers, from various places around the
+ ;; web, with my notes for humans:
+ ;;
+ ;; * n+1 must be prime.
+ ;;
+ ;; Every number has 1 has a factor, which means one of
+ ;; the tests will be to see if 1+(n/1) is prime.
+ ;;
+ ;; * n must be even (except the edge case of 1).
+ ;;
+ ;; We know this because n+1 must be prime, and therefore odd, so n itself
+ ;; must be even.
+ ;;
+ ;; * 2+(n/2) must be prime.
+ ;;
+ ;; Because all candidates are even, they all have 2 as a divisor (see
+ ;; above), and so we can do this check before finding all the divisors.
+ ;;
+ ;; * n must be squarefree.
+ ;;
+ ;; Consider when n is squareful: then there is some prime that occurs more
+ ;; than once in its factorization. Choosing this prime as the divisor for
+ ;; the formula gives us d+(n/d). We know that n/d will still be divisible
+ ;; by d, because we chose a d that occurs multiple times in the
+ ;; factorization. Obviously d itself is divisible by d. Thus our entire
+ ;; formula is divisible by d, and so not prime.
+ ;;
+ ;; Unfortunately this doesn't really help us much, because there's no
+ ;; efficient way to tell if a number is squarefree (see
+ ;; http://mathworld.wolfram.com/Squarefree.html).
+ ;;
+ ;; * We only have to check d <= sqrt(n).
+ ;;
+ ;; For each divisor d of n we know there's a twin divisor d' such that
+ ;; d * d' = n (that's what it MEANS for d to be a divisor of n).
+ ;;
+ ;; If we plug d into the formula we have d + n/d.
+ ;; We know that n/d = d', and so we have d + d'.
+ ;;
+ ;; If we plug d' into the formula we have d' + n/d'.
+ ;; We know that n/d' = d, and so we have d' + d.
+ ;;
+ ;; This means that plugging d or d' into the formula both result in the
+ ;; same number, so we only need to bother checking one of them.
+ (1+ (iterate
+ ;; edge case: skip 2 (candidiate 1), we'll add it at the end
+ (for prime :in-vector (sieve (1+ 100000000)) :from 1)
+ (for candidate = (1- prime))
+ (when (and (check-divisor candidate 2)
+ (prime-generating-integer-p candidate))
+ (summing candidate))))))
;;;; Tests --------------------------------------------------------------------
(def-suite :euler)
@@ -1788,6 +1853,8 @@
(test p79 (is (= 73162890 (problem-79))))
(test p92 (is (= 8581146 (problem-92))))
(test p145 (is (= 608720 (problem-145))))
+(test p357 (is (= 1739023853137 (problem-357))))
-;; (run! :euler)
+(defun run-tests ()
+ (run! :euler))
--- a/src/primes.lisp Tue Aug 08 16:57:14 2017 -0400
+++ b/src/primes.lisp Wed Aug 09 14:31:30 2017 -0400
@@ -1,5 +1,32 @@
(in-package :euler)
+(declaim (ftype (function ((integer 0 #.array-dimension-limit))
+ (simple-array bit (*)))
+ sieve%))
+
+(eval-dammit
+ (defun sieve% (limit)
+ "Return a bit vector of primality for all numbers below `limit`."
+ (iterate
+ (with numbers = (make-array limit :initial-element 1 :element-type 'bit))
+ (for bit :in-vector numbers :with-index n :from 2)
+ (when (= 1 bit)
+ (iterate (for composite :from (* 2 n) :by n :below limit)
+ (setf (aref numbers composite) 0)))
+ (finally
+ (setf (aref numbers 0) 0
+ (aref numbers 1) 0)
+ (return numbers)))))
+
+(defun sieve (limit)
+ "Return a vector of all primes below `limit`."
+ (declare (optimize speed))
+ (iterate (declare (iterate:declare-variables))
+ (for bit :in-vector (sieve% limit) :with-index n :from 2)
+ (when (= 1 bit)
+ (collect n :result-type vector))))
+
+
(defun prime-factorization (n)
"Return the prime factors of `n`."
;; from http://www.geeksforgeeks.org/print-all-prime-factors-of-a-given-number/
@@ -106,6 +133,29 @@
(finally (return t))))))
+;;;; Precomputation -----------------------------------------------------------
+;;; We precompute a bit vector of the primality of the first hundred million
+;;; numbers to make checking primes faster. It'll take up about 12 mb of memory
+;;; and take a few seconds to compute, but saves a lot of brute forcing later.
+
+(defconstant +precomputed-primality-limit+ 100000000)
+
+(defparameter *precomputed-primality-bit-vector*
+ (sieve% +precomputed-primality-limit+))
+
+(deftype precomputed-primality-bit-vector ()
+ `(simple-array bit (,+precomputed-primality-limit+)))
+
+(deftype integer-with-precomputed-primality ()
+ `(integer 0 (,+precomputed-primality-limit+)))
+
+(defun-inline precomputed-prime-p% (n)
+ (declare (optimize speed (debug 0) (safety 1))
+ (type integer-with-precomputed-primality n)
+ (type precomputed-primality-bit-vector *precomputed-primality-bit-vector*))
+ (not (zerop (aref *precomputed-primality-bit-vector* n))))
+
+
(defun primep (n)
"Return (less slowly) whether `n` is prime."
(cond
@@ -113,7 +163,7 @@
((< n 2) nil)
((= n 2) t)
((evenp n) nil)
- ((< n 100000) (brute-force-prime-p n))
+ ((< n +precomputed-primality-limit+) (precomputed-prime-p% n))
(t (miller-rabin-prime-p n))))
(defun unitp (n)
@@ -157,19 +207,6 @@
(primes% (1+ min) max))
-(defun sieve (limit)
- "Return a vector of all primes below `limit`."
- (check-type limit (integer 0 #.array-dimension-limit))
- (iterate
- (with numbers = (make-array limit :initial-element 1 :element-type 'bit))
- (for bit :in-vector numbers :with-index n :from 2)
- (when (= 1 bit)
- (collect n :result-type vector)
- (iterate (for composite :from (* 2 n) :by n)
- (while (< composite limit))
- (setf (aref numbers composite) 0)))))
-
-
(defun nth-prime (n)
"Return the `n`th prime number."
(if (= n 1)
--- a/src/utils.lisp Tue Aug 08 16:57:14 2017 -0400
+++ b/src/utils.lisp Wed Aug 09 14:31:30 2017 -0400
@@ -112,14 +112,22 @@
(sort sequence #'<))
+(defun unsorted-divisors (n)
+ (iterate (for i :from 1 :to (sqrt n))
+ (for (values divisor remainder) = (truncate n i))
+ (when (zerop remainder)
+ (collect i)
+ ;; don't collect the square root twice
+ (unless (= i divisor)
+ (collect divisor)))))
+
+(defun-inlineable divisors-up-to-square-root (n)
+ (loop :for i :from 1 :to (floor (sqrt n))
+ :when (zerop (rem n i))
+ :collect i))
+
(defun divisors (n)
- (sort< (iterate (for i :from 1 :to (sqrt n))
- (when (dividesp n i)
- (collect i)
- (let ((j (/ n i)))
- ;; don't collect the square root twice
- (unless (= i j)
- (collect j)))))))
+ (sort< (unsorted-divisors n)))
(defun proper-divisors (n)
(remove n (divisors n)))