Made it to 100

p100.lisp

@@ -0,0 +1,66 @@
+;;;; Using B for the number of blue discs and N for the total number
+;;;; of discs, starting from P(BB) as given in the problem, it's
+;;;; possible to rearrange terms into a negative Pell's equation with
+;;;; X being (- (* 2 B) 1) and Y being (- (* 2 N) 1) and the equation
+;;;; being (= (- (* X X) (* 2 Y Y)) -1).
+;;;;
+;;;; Negative Pell's equation has a known solution from squaring both
+;;;; sides of the equation.
+
+(defun correct-probability-p (blues total)
+  "Sanity check for solutions"
+  (= 1/2 (* (/ blues total) (/ (1- blues) (1- total)))))
+
+(defun find-next-blues (&optional (start 2))
+  "Naive implemenation that is only used for seeding the first couple of solutions"
+  (loop for total upfrom start
+	for right = (+ 1/4 (/ (- (* total total) total) 2))
+	for root = (sqrt right)
+	for blues = (round (+ 1/2 root))
+	until (= 1/2 (* (/ blues total) (/ (1- blues) (1- total))))
+	finally (return (values blues total))))
+
+(defun x->total (x)
+  "Convert X (or Y) to the count of blue (or total) discs, as described in the opening comment"
+  (/ (1+ x) 2))
+
+(defun total->x (total)
+  "Inverse of X->TOTAL"
+  (1- (* 2 total)))
+
+(defun square (n)
+  (* n n))
+
+(defun compute-negative-pell-solution (n x y k)
+  "Computes the Kth solution to the negative Pell equation with constant N and preceding values in X and Y, returning (LIST XK YK)"
+  (let ((x1 (aref x 1))
+	(y1 (aref y 1))
+	(xk-2 (aref x (- k 2)))
+	(yk-2 (aref y (- k 2))))
+    (list (+ (* xk-2 (square x1))
+	     (* n xk-2 (square y1))
+	     (* 2 n yk-2 y1 x1))
+	  (+ (* yk-2 (square x1))
+	     (* n yk-2 (square y1))
+	     (* 2 xk-2 y1 x1)))))
+
+(defun arranged-probability (&optional (min 1000000000000))
+  (let ((x (make-array 10000 :adjustable t :fill-pointer 1))
+	(y (make-array 10000 :adjustable t :fill-pointer 1)))
+    (flet ((add-solution (blues total)
+	     (vector-push-extend (total->x total) x)
+	     (vector-push-extend (total->x blues) y)
+	     (format t "(~a) ~a / ~a~%" (correct-probability-p blues total) blues total)))
+      ;; Seed initial solutions (note: the problem doesn't mention the
+      ;; initial solution of 3 out of 4)
+      (multiple-value-bind (blues total) (find-next-blues)
+	(add-solution blues total)
+	(multiple-value-bind (blues total) (find-next-blues (1+ total))
+	  (add-solution blues total)))
+      (loop for k upfrom 3
+	    for (xk yk) = (compute-negative-pell-solution 2 x y k)
+	    for total = (x->total xk)
+	    for blues = (x->total yk)
+	    do (add-solution blues total)
+	       (when (> total min)
+		 (return blues))))))