67 lines
2.5 KiB
Common Lisp
67 lines
2.5 KiB
Common Lisp
;;;; 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))))))
|