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