# HG changeset patch # User Steve Losh # Date 1502303490 14400 # Node ID d8750d0cda0f9414085ff5cfdbb43dd6c7b9d6b7 # Parent d08ee014a3983333c5285b8e2566329d144721a2 Problem 357 diff -r d08ee014a398 -r d8750d0cda0f src/euler.lisp --- 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)) diff -r d08ee014a398 -r d8750d0cda0f src/primes.lisp --- 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) diff -r d08ee014a398 -r d8750d0cda0f src/utils.lisp --- 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)))