# HG changeset patch # User Steve Losh # Date 1502988867 14400 # Node ID 540f1cb1093af835403a710cd42c7d5369d13b05 # Parent 252d1614e2d920debfc423dbf9dcd52d343122d2 Problem 323 diff -r 252d1614e2d9 -r 540f1cb1093a euler.asd --- a/euler.asd Thu Aug 10 20:22:09 2017 -0400 +++ b/euler.asd Thu Aug 17 12:54:27 2017 -0400 @@ -10,9 +10,10 @@ :depends-on ( + :1am :anaphora + :cl-pcg :cl-strings - :fiveam :iterate :local-time :losh diff -r 252d1614e2d9 -r 540f1cb1093a package.lisp --- a/package.lisp Thu Aug 10 20:22:09 2017 -0400 +++ b/package.lisp Thu Aug 17 12:54:27 2017 -0400 @@ -3,9 +3,9 @@ :cl :iterate :losh - :5am :euler.quickutils) - (:shadowing-import-from :5am :test)) + (:import-from :1am :is) + (:shadowing-import-from :1am :test)) (defpackage :euler.poker (:use diff -r 252d1614e2d9 -r 540f1cb1093a src/problems.lisp --- a/src/problems.lisp Thu Aug 10 20:22:09 2017 -0400 +++ b/src/problems.lisp Thu Aug 17 12:54:27 2017 -0400 @@ -1714,6 +1714,81 @@ (when reversible (sum (if (= i j) 1 2))))))) +(defun problem-323 () + ;; Let y0, y1, y2,... be a sequence of random unsigned 32 bit integers (i.e. + ;; 0 ≤ yi < 2^32, every value equally likely). + ;; + ;; For the sequence xi the following recursion is given: + ;; + ;; x0 = 0 and + ;; xi = x_i-1 | y_i-1, for i > 0. (| is the bitwise-OR operator) + ;; + ;; It can be seen that eventually there will be an index N such that + ;; xi = 2^32 - 1 (a bit-pattern of all ones) for all i ≥ N. + ;; + ;; Find the expected value of N. + ;; + ;; Give your answer rounded to 10 digits after the decimal point. + (labels + ;; Assuming a perfectly uniform random number generator, each bit of each + ;; randomly-generated bit will be independent of the others. So we can + ;; treat this problem as a set of 32 independent Bernoulli trials, and are + ;; trying to find the expected numbers of iterations required for all 32 + ;; trials to succeed at least once. + ;; + ;; We can start by actually doing the experiments. This will be + ;; impractical for large numbers of parallel trials, but gives us a way to + ;; sanity check our work later. + ((empirical-test () + (iterate (counting t) (until (randomp)))) + (empirical-tests (runners) + (iterate (repeat runners) + (maximizing (empirical-test)))) + (run-empirical-tests (runners) + (coerce (iterate (repeat 1000000) + (averaging (empirical-tests runners))) + 'double-float)) + ;; Running the experiments for 32 parallel trials enough times to get 10 + ;; digits of accuracy is computationally infeasible, so we'll need to use + ;; some math. + ;; + ;; The expected value the problem is looking for is: + ;; + ;; (1 * P(finishing on turn 1)) + + ;; (2 * P(finishing on turn 2)) + + ;; ... + + ;; (n * P(finishing on turn n)) + ;; + ;; We solve this by finding a closed form solution for P(finishing on + ;; turn n), and then adding enough terms to get our ten digits of + ;; precision. + (probability-of-finishing-by (runners n) + ;; The probability of a single run finishing on or before turn N is + ;; given by the cumulative distribution function of the geometric + ;; distribution. + ;; + ;; Because all the trials are independent, we can multiply their + ;; probabilities to find the probability for the group as a whole. + (expt (geometric-cdf 0.5d0 n) runners)) + (probability-of-finishing-on (runners n) + ;; The probability of finishing EXACTLY on turn N can be found by + ;; taking N's cumulative probability and subtracting N-1's cumulative + ;; probability. + (- (probability-of-finishing-by runners n) + (probability-of-finishing-by runners (1- n)))) + (expected-value (runners) + (iterate + (for n :from 1 :below 1000) + (for prob = (probability-of-finishing-on runners n)) + (summing (* n prob)))) + (sanity-check (runners) + (assert (= (round-to (run-empirical-tests runners) 1) + (round-to (expected-value runners) 1))))) + (when nil + (iterate (for i :from 1 :to 6) + (sanity-check i))) + (round-to (expected-value 32) 10))) + (defun problem-345 () ;; We define the Matrix Sum of a matrix as the maximum sum of matrix elements ;; with each element being the only one in his row and column. For example, @@ -1820,9 +1895,6 @@ ;;;; Tests -------------------------------------------------------------------- -(def-suite :euler) -(in-suite :euler) - (test p1 (is (= 233168 (problem-1)))) (test p2 (is (= 4613732 (problem-2)))) (test p3 (is (= 6857 (problem-3)))) @@ -1892,8 +1964,10 @@ (test p79 (is (= 73162890 (problem-79)))) (test p92 (is (= 8581146 (problem-92)))) (test p145 (is (= 608720 (problem-145)))) +(test p323 (is (= 6.3551758451d0 (problem-323)))) +(test p345 (is (= 13938 (problem-345)))) (test p357 (is (= 1739023853137 (problem-357)))) (defun run-tests () - (run! :euler)) + (1am:run)) diff -r 252d1614e2d9 -r 540f1cb1093a src/utils.lisp --- a/src/utils.lisp Thu Aug 10 20:22:09 2017 -0400 +++ b/src/utils.lisp Thu Aug 17 12:54:27 2017 -0400 @@ -509,3 +509,29 @@ (aref matrix source-row source-col))) (finally (return result)))) + +(defun geometric-pmf (success-probability trials) + "The probability mass function of the geometric distribution. + + Returns the probability that exactly `trials` trials will be required to see + the first success in a series of Bernoulli trials with `success-probability`. + + " + (* (expt (- 1 success-probability) (1- trials)) + success-probability)) + +(defun geometric-cdf (success-probability trials) + "The cumulative distribution function of the geometric distribution. + + Returns the probability that `trials` or fewer trials will be required to see + the first success in a series of Bernoulli trials with `success-probability`. + + " + (- 1 (expt (- 1 success-probability) trials))) + + +(defun round-to (number precision &optional (rounder #'round)) + ;; http://www.codecodex.com/wiki/Round_a_number_to_a_specific_decimal_place#Common_Lisp + (coerce (let ((div (expt 10 precision))) + (/ (funcall rounder (* number div)) div)) + (type-of number)))