title: Monte Carlo integration of the Gaussian function f(x) = exp(-x^2/2) in C incorrect output
body: <p>I'm writing a short program to approximate the definite integral of the gaussian function f(x) = exp(-x^2/2), and my codes are as follows: </p>
  
  <pre><code>#include &lt;stdio.h&gt;
  #include &lt;stdlib.h&gt;
  #include &lt;math.h&gt;
  
  double gaussian(double x) {
      return exp((-pow(x,2))/2);
  }
  
  int main(void) {
      srand(0);
      double valIntegral, yReal = 0, xRand, yRand, yBound;
      int xMin, xMax, numTrials, countY = 0;
  
      do {
          printf("Please enter the number of trials (n): ");
          scanf("%d", &amp;numTrials);
          if (numTrials &lt; 1) {
              printf("Exiting.\n");
              return 0;
          }  
          printf("Enter the interval of integration (a b): ");
          scanf("%d %d", &amp;xMin, &amp;xMax);      
          while (xMin &gt; xMax) { //keeps looping until a valid interval is entered
              printf("Invalid interval!\n");
              printf("Enter the interval of integration (a b): ");
              scanf("%d %d", &amp;xMin, &amp;xMax);
          }
          //check real y upper bound
          if (gaussian((double)xMax) &gt; gaussian((double)xMin))
              yBound = gaussian((double)xMax);
          else 
              yBound = gaussian((double)xMin);
          for (int i = 0; i &lt; numTrials; i++) {
              xRand = (rand()% ((xMax-xMin)*1000 + 1))/1000.00 + xMin; //generate random x value between xMin and xMax to 3 decimal places             
              yRand = (rand()% (int)(yBound*1000 + 1))/1000.00; //generate random y value between 0 and yBound to 3 decimal places
              yReal = gaussian(xRand);
              if (yRand &lt; yReal) 
                  countY++;
          }
          valIntegral = (xMax-xMin)*((double)countY/numTrials);
          printf("Integral of exp(-x^2/2) on [%.3lf, %.3lf] with n = %d trials is: %.3lf\n\n", (double)xMin, (double)xMax, numTrials, valIntegral);
  
          countY = 0; //reset countY to 0 for the next run
      } while (numTrials &gt;= 1);
  
      return 0;
  }
  </code></pre>
  
  <p>However, the outputs from my code doesn't match the solutions. I tried to debug and print out all xRand, yRand and yReal values for 100 trials (and checked yReal value with particular xRand values with Matlab, in case I had any typos), and those values didn't seem to be out of range in any way... I don't know where my mistake is. </p>
  
  <p>The correct output for # of trials = 100 on [0, 1] is 0.810, and mine is 0.880; correct output for # of trials = 50 on [-1, 0] is 0.900, and mine was 0.940. Can anyone find where I did wrong? Thanks a lot. </p>
  
  <p>Another question is, I can't find a reference to the use of following code: </p>
  
  <pre><code>double randomNumber = rand() / (double) RAND MAX;
  </code></pre>
  
  <p>but it was provided by the instructor and he said it would generate a random number from 0 to 1. Why did he use <code>'/'</code> instead of <code>'%'</code> after <code>"rand()"</code>?</p>