# HG changeset patch
# User Steve Losh <steve@stevelosh.com>
# Date 1533053913 0
# Node ID 99fa06464617ebc7d2cad78be1ba63fb4907f86c
# Parent  9a8362bd95aeedb8712c4b1515b7b5e7320e3b1f
More random stuff

diff -r 9a8362bd95ae -r 99fa06464617 src/primes.lisp
--- a/src/primes.lisp	Tue Jul 24 15:11:16 2018 +0000
+++ b/src/primes.lisp	Tue Jul 31 16:18:33 2018 +0000
@@ -6,25 +6,34 @@
 
 (eval-dammit
   (defun sieve% (limit)
-    "Return a bit vector of primality for all numbers below `limit`."
+    "Return a bit vector of primality for all odd 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)
+      (with length = (truncate limit 2))
+      (with numbers = (make-array length :initial-element 1 :element-type 'bit))
+      (for bit :in-vector numbers :with-index n :from 1)
       (when (= 1 bit)
-        (iterate (for composite :from (* 2 n) :by n :below limit)
-                 (setf (aref numbers composite) 0)))
+        (iterate
+          (with step = (1+ (* 2 n)))
+          (for composite :from (+ n step) :by step :below length)
+          (setf (aref numbers composite) 0)))
       (finally
-        (setf (aref numbers 0) 0
-              (aref numbers 1) 0)
+        (setf (aref numbers 0) 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))))
+  (check-type limit (integer 0 #.array-dimension-limit))
+  (if (< limit 2)
+    (vector)
+    (iterate
+      (declare (iterate:declare-variables))
+      (if-first-time
+        (collect 2 :into result :result-type vector))
+      (for bit :in-vector (sieve% limit) :with-index n :from 1)
+      (when (= 1 bit)
+        (collect (1+ (* 2 n)) :into result :result-type vector))
+      (finally (return result)))))
 
 
 (defun prime-factors (n)
@@ -201,17 +210,17 @@
 
 
 ;;;; Precomputation -----------------------------------------------------------
-;;; We precompute a bit vector of the primality of the first hundred million
+;;; We precompute a bit vector of the primality of the first hundred million odd
 ;;; 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)
+(defconstant +precomputed-primality-limit+ 200000000)
 
 (defparameter *precomputed-primality-bit-vector*
   (sieve% +precomputed-primality-limit+))
 
 (deftype precomputed-primality-bit-vector ()
-  `(simple-array bit (,+precomputed-primality-limit+)))
+  `(simple-array bit (,(truncate +precomputed-primality-limit+ 2))))
 
 (deftype integer-with-precomputed-primality ()
   `(integer 0 (,+precomputed-primality-limit+)))
@@ -220,7 +229,7 @@
   (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))))
+  (not (zerop (aref *precomputed-primality-bit-vector* (truncate n 2)))))
 
 
 (defun primep (n)
diff -r 9a8362bd95ae -r 99fa06464617 src/utils.lisp
--- a/src/utils.lisp	Tue Jul 24 15:11:16 2018 +0000
+++ b/src/utils.lisp	Tue Jul 31 16:18:33 2018 +0000
@@ -540,6 +540,78 @@
            (counting t)))
 
 
+;;;; Matrices and Vectors -----------------------------------------------------
+(defun count-rows (matrix)
+  (array-dimension matrix 0))
+
+(defun count-cols (matrix)
+  (array-dimension matrix 1))
+
+(defun mat* (m n)
+  (assert (= (count-cols m) (count-rows n))
+      () "Cannot multiply incompatibly-sized matrices.")
+  (let ((rows (count-rows m))
+        (cols (count-cols n)))
+    (iterate
+      (with numbers = (count-cols m))
+      (with result = (make-array (list rows cols) :element-type 'number))
+      (for-nested ((row :from 0 :below rows)
+                   (col :from 0 :below cols)))
+      (setf (aref result row col)
+            (iterate (for i :below numbers)
+                     (summing (* (aref m row i)
+                                 (aref n i col)))))
+      (finally (return result)))))
+
+(defun mv* (matrix vector)
+  (iterate
+    (with (rows cols) = (array-dimensions matrix))
+    (initially (assert (= cols (length vector))))
+    (with result = (make-array rows :initial-element 0))
+    (for row :from 0 :below rows)
+    (iterate (for col :from 0 :below cols)
+             (for v = (aref vector col))
+             (for a = (aref matrix row col))
+             (incf (aref result row)
+                   (* v a)))
+    (finally (return result))))
+
+
+(defun mat2 (a b c d)
+  (let ((result (make-array (list 2 2) :element-type 'number)))
+    (setf (aref result 0 0) a
+          (aref result 0 1) b
+          (aref result 1 0) c
+          (aref result 1 1) d)
+    result))
+
+
+(defun vec2 (x y)
+  (vector x y))
+
+(defun vx (vec2)
+  (aref vec2 0))
+
+(defun vy (vec2)
+  (aref vec2 1))
+
+(defun vec2+ (a b)
+  (vec2 (+ (vx a) (vx b))
+        (+ (vy a) (vy b))))
+
+(defun vec2- (a b)
+  (vec2 (- (vx a) (vx b))
+        (- (vy a) (vy b))))
+
+(defun vec2-dot (a b)
+  (+ (* (vx a) (vx b))
+     (* (vy a) (vy b))))
+
+(defun vec2= (a b)
+  (and (= (vx a) (vx b))
+       (= (vy a) (vy b))))
+
+
 ;;;; Fibonacci ----------------------------------------------------------------
 (defmacro-driver (FOR var IN-FIBONACCI _)
   (declare (ignore _))
@@ -558,6 +630,24 @@
            (for i :in-fibonacci t)
            (collect i)))
 
+(defun nth-fibonacci (n)
+  ;; https://blog.paulhankin.net/fibonacci/
+  (labels
+      ((mexpt (matrix exponent)
+         (cond ((zerop exponent) (mat2 1 0 0 1))
+               ((= 1 exponent) matrix)
+               ((evenp exponent)
+                (mexpt (mat* matrix matrix)
+                       (truncate exponent 2)))
+               (t (mat* matrix
+                        (mexpt (mat* matrix matrix)
+                               (truncate exponent 2)))))))
+    (-<> (mat2 1 1
+               1 0)
+      (mexpt <> n)
+      (mv* <> (vec2 1 1))
+      (aref <> 1))))
+
 
 ;;;; Factorial ----------------------------------------------------------------
 (defun factorial% (n)
@@ -691,18 +781,6 @@
 (define-modify-macro adjoinf (item &rest keyword-args) adjoin%)
 
 
-(defun mv* (matrix vector)
-  (iterate
-    (with (rows cols) = (array-dimensions matrix))
-    (initially (assert (= cols (length vector))))
-    (with result = (make-array rows :initial-element 0))
-    (for row :from 0 :below rows)
-    (iterate (for col :from 0 :below cols)
-             (for v = (aref vector col))
-             (for a = (aref matrix row col))
-             (incf (aref result row)
-                   (* v a)))
-    (finally (return result))))
 
 
 (defun pythagorean-triplet-p (a b c)
@@ -896,32 +974,6 @@
   (* precision (round number precision)))
 
 
-;;;; Yet Another Vec2 Implementation™ -----------------------------------------
-(defun vec2 (x y)
-  (vector x y))
-
-(defun vx (vec2)
-  (aref vec2 0))
-
-(defun vy (vec2)
-  (aref vec2 1))
-
-(defun vec2+ (a b)
-  (vec2 (+ (vx a) (vx b))
-        (+ (vy a) (vy b))))
-
-(defun vec2- (a b)
-  (vec2 (- (vx a) (vx b))
-        (- (vy a) (vy b))))
-
-(defun vec2-dot (a b)
-  (+ (* (vx a) (vx b))
-     (* (vy a) (vy b))))
-
-(defun vec2= (a b)
-  (and (= (vx a) (vx b))
-       (= (vy a) (vy b))))
-
 
 ;;;; A* Search ----------------------------------------------------------------
 (defstruct path