d8750d0cda0f

Problem 357
[view raw] [browse files]
author Steve Losh <steve@stevelosh.com>
date Wed, 09 Aug 2017 14:31:30 -0400
parents d08ee014a398
children 9230e81d302c
branches/tags (none)
files src/euler.lisp src/primes.lisp src/utils.lisp

Changes

--- 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)))