Probability Simulation Error not Converging

Monte Carlo methods tend to converge slowly --- the error you expect after n simulations is proportional to 1/sqrt(n). Five digits of accuracy after 10^9 trials seems about right, actually. There's no numerical voodoo going on here.

If the interviewer was talking about straight-up Monte Carlo sampling as you've done, it's...implausible that he could get nine digits of accuracy after a million trials.

because the probabilities are given as rational numbers (with small integers in the denominator), you can view the possible situations as a cube of dimensions 7x5x4 (which makes 140 (product of the denominators) subcubes). Rather than jumping around randomly, you can visit each subcube explicitly as follows and get the exact number in 140 iterations:

#include <cstdio>
#include <cctype>
#include <ctime>
#include <random>

int main()
{
  std::size_t total_success = 0, num_trials = 0;

  for (unsigned a = 1; a <= 7; ++a)
  {
    unsigned success_a = 0;

    if (a <= 6)
      // a hits 6 out of 7 times
      success_a = 1;

    for (unsigned b = 1; b <= 5; ++b)
    {
      unsigned success_b = 0;

      if (b <= 4)
        // b hits 4 out of 5 times
        success_b = 1;

        for (unsigned c = 1; c <= 4; ++c)
        {
          unsigned success_c = 0;

          // c hits 3 out of 4 times
          if (c <= 3)
            success_c = 1;

          // count cases where exactly two of them hit
          if (success_a + success_b + success_c == 2)
            ++total_success;

          ++num_trials;

        } // loop over c
    } // loop over b
  } // loop over a

  double prob = (total_success * 1.0) / num_trials;

  printf("Pr(...) = %12.10f  error:%15.13f\n",
         prob,
         std::abs((27.0/70.0) - prob));

   return 0;
}

FWIW the following Java seems to be converging on the predicted answer from above at about the rate you would expect (it calculates the standard deviation of the worst case error)

import java.util.Random;
import java.security.SecureRandom;
/** from question in Stack Overflow */
public class SoProb
{
  public static void main(String[] s)
  {
    long seed = 42;


/*
In an interview, I was given the following problem to solve initially using pen/paper, then via a program to verify the result.

The question is as follows:

There are three people A,B and C. Each person is capable of hitting a target with a probability of 6/7, 4/5 and 3/4 respectively. What is the probability that if they were to each fire one shot that exactly two of them will hit the target?

The answer is:

P(...) = P(A)*P(B)*(1-P(C)) +
         P(B)*P(C)*(1-P(A)) +
         P(C)*P(A)*(1-P(B))
       = 27.0/70.0
       = 38.57142857142857142857142857142857142857....%

Below is my solution to the problem:
*/

/*
int main()
{
   std::mt19937 engine(time(0));
*/

   Random r = new Random(seed);
   // Random r = new SecureRandom(new byte[] {(byte)seed});
   // std::uniform_real_distribution<double> uniform_real(0.0,1.0);

   double prA = (6.0 / 7.0);
   double prB = (4.0 / 5.0);
   double prC = (3.0 / 4.0);
   // double prB = (6.0 / 7.0);
   // double prC = (4.0 / 5.0);
   // double prA = (3.0 / 4.0);

   double pp = prA*prB*(1-prC) +
         prB*prC*(1-prA) +
         prC*prA*(1-prB);
   System.out.println("Pp " + pp);
   System.out.println("2870 " + (27.0 / 70.0));

   // std::size_t trails = 4000000000;
   int trails = Integer.MAX_VALUE;
   // std::size_t total_success = 0;
   int total_success = 0;

   int aCount = 0;
   int bCount = 0;
   int cCount = 0;

   int pat3 = 0; // A, B
   int pat5 = 0; // A, C
   int pat6 = 0; // B, C
   double pat3Prob = prA * prB * (1.0 - prC);
   double pat5Prob = prA * prC * (1.0 - prB);
   double pat6Prob = prC * prB * (1.0 - prA);
   System.