cb34c7bc34fc

Problem 346
[view raw] [browse files]
author Steve Losh <steve@stevelosh.com>
date Fri, 13 Oct 2017 01:32:55 -0400
parents 393d1f8dd754
children 813d3e887794
branches/tags (none)
files src/primes.lisp src/problems.lisp src/utils.lisp

Changes

--- a/src/primes.lisp	Thu Oct 12 21:29:54 2017 -0400
+++ b/src/primes.lisp	Fri Oct 13 01:32:55 2017 -0400
@@ -28,7 +28,15 @@
 
 
 (defun prime-factors (n)
-  "Return the prime factors of `n`."
+  "Return the prime factors of `n`.
+
+  The result will be a set of primes.  For example:
+
+    (prime-factors 60)
+    ; =>
+    (2 3 5)
+
+  "
   (let ((result (list)))
     (iterate (while (evenp n)) ; handle 2, the only even prime factor
              (if-first-time
@@ -45,7 +53,15 @@
     (nreverse result)))
 
 (defun prime-factorization (n)
-  "Return the prime factorization of `n`."
+  "Return the prime factorization of `n` as a flat list.
+
+  The result will be a list of primes, with duplicates.  For example:
+
+    (prime-factorization 60)
+    ; =>
+    (2 2 3 5)
+
+  "
   ;; from http://www.geeksforgeeks.org/print-all-prime-factors-of-a-given-number/
   (let ((result (list)))
     (iterate (while (evenp n)) ; handle 2, the only even prime factor
@@ -60,6 +76,33 @@
       (push n result))
     (nreverse result)))
 
+(defun prime-factorization-pairs (n)
+  "Return the prime factorization of `n` as a list of prime/exponent conses.
+
+  The result will be a list of `(prime . exponent)` conses.  For example:
+
+    (prime-factorization-pairs 60)
+    ; =>
+    ((2 . 2) (3 . 1) (5 . 1))
+
+  "
+  (let ((result (list)))
+    (iterate (while (evenp n)) ; handle 2, the only even prime factor
+             (if-first-time
+               (push (cons 2 1) result)
+               (incf (cdar result)))
+             (setf n (/ n 2)))
+    (iterate
+      (for i :from 3 :to (sqrt n) :by 2) ; handle odd (prime) divisors
+      (iterate (while (dividesp n i))
+               (if-first-time
+                 (push (cons i 1) result)
+                 (incf (cdar result)))
+               (setf n (/ n i))))
+    (when (> n 2) ; final check in case we ended up with a prime
+      (push (cons n 1) result))
+    (nreverse result)))
+
 
 (defun expmod (base exp m)
   "Return base^exp % m quickly."
--- a/src/problems.lisp	Thu Oct 12 21:29:54 2017 -0400
+++ b/src/problems.lisp	Fri Oct 13 01:32:55 2017 -0400
@@ -2084,6 +2084,57 @@
       (for (row . col) :in (euler.hungarian:find-minimal-assignment matrix))
       (summing (- (aref matrix row col))))))
 
+(defun problem-346 ()
+  ;; The number 7 is special, because 7 is 111 written in base 2, and 11 written
+  ;; in base 6 (i.e. 7_10 = 11_6 = 111_2). In other words, 7 is a repunit in at
+  ;; least two bases b > 1.
+  ;;
+  ;; We shall call a positive integer with this property a strong repunit. It
+  ;; can be verified that there are 8 strong repunits below 50:
+  ;; {1,7,13,15,21,31,40,43}.  Furthermore, the sum of all strong repunits below
+  ;; 1000 equals 15864.
+  ;;
+  ;; Find the sum of all strong repunits below 10^12.
+  (iterate
+    ;; Let's start with a few observations.  First, 1 is a repunit in ALL bases,
+    ;; and 2 isn't a repunit in any base.
+    ;;
+    ;; Next, EVERY number `n` greater than 2 is a repunit in base `n - 1`,
+    ;; because it will be `11` in that base.  So every number already fulfills
+    ;; half the requirements of the problem.  This means we only need to find
+    ;; and sum all the numbers are repunits in at least one more base, and
+    ;; they'll be of the form `11[1]+` in that base.
+    ;;
+    ;; Instead of iterating through each number trying to search for a base that
+    ;; it happens to be a repunit in, it's faster if we iterate through the
+    ;; possible BASES instead and generate all their repunits below the limit.
+    ;; Observe how new repunits can be iteratively formed for a given base `b`:
+    ;;
+    ;;   R₁ = 1    = 1
+    ;;   R₂ = 11   = 1 + b¹
+    ;;   R₃ = 111  = 1 + b¹ + b²
+    ;;   R₄ = 1111 = 1 + b¹ + b² + b³
+    ;;   ...
+    ;;
+    ;; We can use a set to store the results to avoid overcounting numbers that
+    ;; are repunits in more than two bases.
+    (with limit = (expt 10 12))
+    (with strong-repunits = (make-hash-set :initial-contents '(1)))
+    (for base :from 2 :to (sqrt limit))
+    (iterate
+      ;; To check a particular base we'll start at number 111 in that base,
+      ;; because we've already accounted for 1 and 11 in our previous
+      ;; assumptions.
+      ;;
+      ;; This is why we only iterate up to √b in the outer loop: the first
+      ;; number we're going to be checking is `111 = 1 + b¹ + b²`, and so any
+      ;; base higher than √b would immediately exceed the limit.  We could
+      ;; probably narrow the range a bit further, but it's plenty fast already.
+      (for n :first (+ 1 base (expt base 2)) :then (+ 1 (* base n)))
+      (while (< n limit))
+      (hset-insert! strong-repunits n))
+    (finally (return (hset-reduce strong-repunits #'+)))))
+
 (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.
@@ -2301,6 +2352,7 @@
 (test p315 (is (= 13625242 (problem-315))))
 (test p323 (is (= 6.3551758451d0 (problem-323))))
 (test p345 (is (= 13938 (problem-345))))
+(test p346 (is (= 336108797689259276 (problem-346))))
 (test p357 (is (= 1739023853137 (problem-357))))
 (test p387 (is (= 696067597313468 (problem-387))))
 
--- a/src/utils.lisp	Thu Oct 12 21:29:54 2017 -0400
+++ b/src/utils.lisp	Fri Oct 13 01:32:55 2017 -0400
@@ -290,11 +290,20 @@
   (remove n (divisors n)))
 
 (defun count-divisors (n)
-  (+ (* 2 (iterate (for i :from 1 :below (sqrt n))
-                   (counting (dividesp n i))))
-     (if (squarep n)
-       1
-       0)))
+  ;; From *Recreations in the Theory of Numbers: The Queen of Mathematics
+  ;; Entertains* by Albert Beiler.
+  ;;
+  ;; To find the number of divisors of n, we first find its prime factorization:
+  ;;
+  ;;   p₁^k₁ × p₂^k₂ × ... × pₓ^kₓ
+  ;;
+  ;; The divisors of n are each formed by some combination of these numbers.  To
+  ;; find the total number of divisors we take ∏(k+1) for each exponent k in the
+  ;; factorization.  This is because there are k₁+1 options for how many of the
+  ;; first prime to include (because 0 is an option), k₂+1 options for the
+  ;; second prime, etc.
+  (iterate (for (nil . exponent) :in (prime-factorization-pairs n))
+           (multiplying (1+ exponent))))
 
 
 (defmacro-driver (FOR var IN-COLLATZ n)