content: {- 
   - Solution to Project Euler problem 323
   - Copyright (c) Project Nayuki. All rights reserved.
   - 
   - https://www.nayuki.io/page/project-euler-solutions
   - https://github.com/nayuki/Project-Euler-solutions
   -}
  
  import Numeric (showFFloat)
  
  
  {- 
   - Suppose that n 32-bit integers have been OR'd together.
   - For any arbitrary bit:
   -   The probability that it is 0 is 1/2^n.
   -   The probability that it is 1 is 1 - 1/2^n.
   - Thus for the entire number:
   -   The probability that all bits are 1 is (1 - 1/2^n)^32.
   -     This is the cumulative distribution function that we want.
   -   The probability that some bit is 0 is 1 - (1 - 1/2^n)^32.
   - 
   - The probability density function is simply pdf(n) = cdf(n) - cdf(n-1).
   - So the expected value of the index where the number becomes all 1's is
   - sum(n * pdf(n) for n = 0 to infinity).
   - 
   - This program computes an approximate answer using floating-point, not guaranteed to be correct.
   - Moreover, it truncates the infinite series by ignoring insignificant contributions to the sum.
   - However, the Mathematica version of the solution is exact.
   -}
  main = putStrLn (show ans)
  ans = (showFFloat (Just 10) expectedValue) ""
  expectedValue = sum $ zipWith (*) (takeWhile (> 1e-20) (tail pdf)) [1..]
  
  -- Probability density function
  pdf = 0 : (zipWith (-) (tail cdf) cdf)
  
  -- Cumulative distribution function
  cdf = [(1 - 2**(-n))^32 | n <- [0..]] :: [Double]
  