// ////////////////////////////////////////////////////////
// # Title
// Unfair wager
//
// # URL
// https://projecteuler.net/problem=436
// http://euler.stephan-brumme.com/436/
//
// # Problem
// Julie proposes the following wager to her sister Louise.
// She suggests they play a game of chance to determine who will wash the dishes.
// For this game, they shall use a generator of independent random numbers uniformly distributed between 0 and 1.
// The game starts with `S = 0`.
// The first player, Louise, adds to S different random numbers from the generator until `S > 1` and records her last random number `x`.
// The second player, Julie, continues adding to S different random numbers from the generator until `S > 2` and records her last random number `y`.
// The player with the highest number wins and the loser washes the dishes, i.e. if `y > x` the second player wins.
//
// For example, if the first player draws `0.62` and `0.44`, the first player turn ends since `0.62+0.44 > 1` and `x = 0.44`.
// If the second players draws `0.1`, `0.27` and `0.91`, the second player turn ends since `0.62+0.44+0.1+0.27+0.91 > 2` and `y = 0.91`. Since `y > x`, the second player wins.
//
// Louise thinks about it for a second, and objects: "That's not fair".
// What is the probability that the second player wins?
// Give your answer rounded to 10 places behind the decimal point in the form 0.abcdefghij
//
// # Monte-Carlo simulation written by
// Stephan Brumme
// November 2017
//
// # Idea
// I wrote a simple Monte-Carlo simulation which finds the first 4 digits in a few seconds: 0.5276
//
// My current thinking is that I could implement a better numerical approximation of the probability of drawing `x` as the last number (that is, some `p(x)`).
// Julie's number would be `p(0)` and Louise's `p(p(0))`. The result is the probability that `p(0) < p(p(0))`.
//
// However, I'm not sure whether that would be enough to compute the first 10 digits reliably.

#include <iostream>

// a simple pseudo-random number generator, returning a value between 0 and 1 (exclusive)
// (produces the same result no matter what compiler you have - unlike rand() from math.h)
double myrand()
{
  // based on code from problem 227, modified to return a value between 0 and 1 (exclusive)
  static unsigned long long seed = 0;
  seed = 6364136223846793005ULL * seed + 1;

  const unsigned int ValidBits = 31; // maybe not all bits are random, but it's not a high-quality algorithm anyway
  const unsigned int Scaling   = 1 << ValidBits;
  auto result = ((seed >> 30) % Scaling) / double(Scaling);
  return result;
}

double monteCarlo(unsigned int iterations)
{
  // count how often player 2 wins
  unsigned int twoWins = 0;

  for (unsigned int i = 0; i < iterations; i++)
  {
    // sum of all random numbers
    double total = 0;
    // last random number of player 1
    double lastOne = 0;
    while (total < 1)
    {
      lastOne = myrand();
      total += lastOne;
    }

    // last random number of player 1
    double lastTwo = 0;
    while (total < 2)
    {
      lastTwo = myrand();
      total += lastTwo;
    }

    if (lastTwo > lastOne)
      twoWins++;
  }

  return twoWins / double(iterations);
}

int main()
{
  // the first 4 digits should be correct - but I have to find the first 10 digits, which is impossible with Monte-Carlo
  while (true)
    std::cout << monteCarlo(100000000) << std::endl;
  return 0;
}