Tenth Week Examples ---
   Monte Carlo methods.
  Arrays-Pointers
  Sieve of Eratosthenes

Following D. Milicic Tenth Week Examples and Eleventh Week Examples.

Monte Carlo Methods

The random number generator picks a new random number every time its called. (The numbers aren't really random. Actually the numbers are determined by some rule one after another, but they have good statistical properties and can be used for integration and to simulate randomness.) The numbers are uniformly distributed integers that fall between 0 and RAND_MAX which is 2147483647 in this system. To get a random number that is uniformly distributed between zero and one, we compute 

 
double x;
 . . .
x = rand () / 2147483648.0;

Program to compute Pi

Here we use the Monte Carlo method to compute the area of the part inside the unit circle relative to the area of the unit square. π is four times this fraction. 
/*****************************************************************************
Treibergs                                                        Feb. 24, 2006 

Monte Carlo Pi

montecarlopi.c
*****************************************************************************/
      

# include  <stdio.h>
# include  <stdlib.h>
# include  <math.h>
# include  <time.h>




int 
main ( void ) 
{
    double x, y, p ,s=0.0,  ranmaxpo;
    int i, j, k, n = 1000000;
    ranmaxpo = 1.0 + RAND_MAX;

    srand ( (unsigned int)time ( NULL ) );    /* seed for rand is time of day  */

    printf ( "  n\t\t\t\tMonte Carlo Pi             MCP - pi\n" );

    for ( k=1; k<= 10; k++ )
    {
        j=0;
        for ( i=1; i <= n; i=i+1 )
        {
            x = rand () / ranmaxpo;
            y = rand () / ranmaxpo;
            if ( x * x  +  y * y  <= 1 )
                    j++;
        }
        p = 4.0 * j / n;
        s = s + p;
        printf ( "%4d The approximate for pi is %21.15f %21.15f\n", k, p, p-M_PI );
     }

     p = s / 10.0;

     printf ( " Average approximate for pi is %21.15f %21.15f\n", p, p-M_PI );
     printf ( "                          pi = %21.15f \n", M_PI );

     return EXIT_SUCCESS;
}


/* Milicic

 Monte-Carlo method for calculating pi

montecarlo.c                                        */


#include <stdio.h>
#include <stdlib.h>
#define RANDOM_MAX 2147483647

int main() {
 
    int i=0, j=0;
    float x, y;

    while(1) {
      x = (double)random()/(double)RANDOM_MAX;
      y = (double)random()/(double)RANDOM_MAX;
    if ((x*x+y*y) < 1)
       j++;
    i++;
   if (i%1000000 == 1)
   printf("%f\n",4*(float)j/(float)i);
   }
}

Monte Carlo Method for the integral of f(x)

Here is a program that uses the Monte Carlo method to compute the integral of a function. We compute the area that is both under the curve  0 ≤ y ≤ f(x) and in the box  a ≤ x ≤ b  and  0 ≤ y ≤ c . 
/*****************************************************************************
Treibergs                                                         Mar. 6, 2006 

Monte Carlo integral. We assume that the function is given between a <= x <= b
We compute the area of the { (x,y) : a <= x <= b  and  0 <= y <= min(c,f(x)) }
by counting the relative number of random points in the rectangle [a,b]x[0,c]
that fall in the set.  

mo_ca_int.c
*****************************************************************************/
      

# include <stdio.h>
# include <stdlib.h>
# include <math.h>


double
f( double x );

int 
main ( void ) 
{
    double a, b, c, d, x, y, p, q, r, sum=0.0, ans, ranmaxpo;
    int i, j, k, m=15, n=2000000;
    ranmaxpo = 1.0 + RAND_MAX;
    
    printf ( " Enter the left, right and upper bounds : ");
    scanf ( "%lf %lf %lf", &a, &b, &c );
    a = fabs ( a );
    b = fabs ( b );
    c = fabs ( c );
    if ( b < a )
    {
           d = a;
           a = b;
           c = d;
     }
     printf ( "\n Monte Carlo Int of min(f(x),%f) over %f <= x <= %f\n", c, a, b); 
     
     d = pow( c, 1.0 / 3.0 );
     if ( d >= b)
         ans = ( b * b  +  a * a) * ( b * b  -  a * a ) / 4.0;
     else if ( d <= a )
     {
         ans = c * ( b  -  a );
         printf(" Note that  f(x) > %f  for  %f < x <= %f\n", c, a, b );
     }
     else
     {
         ans =  ( d * d  +  a * a ) * ( d * d  -  a * a ) / 4.0  +  ( b  -  d ) * c;
         printf(" Note that  f(x) > %f  for  %f < x <= %f\n", c, d, b );
     }
     
     
    printf( "\n\t n\tApproximate integral\t\tError\n");
    for(k=1;k<= m; k++)
    {
    j = 0;
    q = c / ranmaxpo;
    r = ( b  -  a ) / ranmaxpo;
        for ( i = 1; i <= n; i = i+1 )
        {
        x = a  +  r * rand ();
        y = q * rand ();
        if ( f(x) > y  )
               j++;
        } 
    p = ( b  -  a ) * c * j / n;
    
   printf ( "%12d%22.15f%22.15f\n", k, p, p  -  ans );
   sum = sum  +  p;
   }
   p = sum / m;
   
   printf ( "Average int =%21.15f%22.15f\n", p, ans  -  p );
   printf ( "Actual int = %21.15f", ans );
   printf ( "      Number of points is %ld\n", (long int)n * (long int)m );
   
   return EXIT_SUCCESS;
}

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

double
f( double x )
{
       return  x * x * x;
}

Arrays

An example of an array is a vector in three space, which is given by listing three coordinates. Here is a program that enters vectors and computes the distance between them.

Distance 


/* Milicic

 Distance between points in 3-space 

dist.c                                                                   */

#include <stdio.h>
#include <math.h>

int main(){
  int i;

  double  x[3], y[3], sum, dist;
  
  printf("The coordinates of the first point: ");
   for ( i = 0 ; i < 3  ; i++ )
     scanf("%lf",&x[i]);

  printf("\nThe coordinates of the second point: ");
   for ( i = 0 ; i < 3  ; i++ )
     scanf("%lf",&y[i]);
   
   sum = 0;

   for( i = 0 ; i < 3 ; i++)
     sum = sum + (x[i] - y[i])*(x[i] - y[i]);

   dist = sqrt(sum);
     
      printf("\nThe distance is %lf.\n",dist);

}

Sieve of Eratosthenes

Here is another application of arrays. This program has a few bells and whistles. The idea to find the primes is to systematically cross out all nonprime numbers from the list of all numbers. The basic operation is crossing off all numbers that are multiples (hence nonprime) in an integer array. We then print out the primes in two ways. The first just lists the uncrossed numbers. In the second, we tabulate all numbers from 1 to 200, but replace the composite (nonprime) numbers by dots. The array thus printed gives a sense of where the prime numbers occur among the first thousand numbers.

There are many websites that animate the sieve. Check out my colleague's, Peter Alfeld's, website http://www.math.utah.edu/~alfeld/Eratosthenes.html 

/*****************************************************************************
Treibergs                                                        Feb. 24, 2006 

Sieve of Eratosthenes

Starting at k=2, cross out all multiples of 2 {4,6,8,10,12,...} in the list of 
the first  n  numbers.  (n=1000 here.)
Then keep  3  and cross out the multiples {6,9,12,15,18,21,...}
Since four and its multiples have been crossed out already, we can skip four.
Then do fives and so on as long as k*k <= n. The remaining numbers are prime. 

Here we have an integer array m[1000] that is initialized as ones. Then 
we loop through the multiples, and "cross out the number" j  by setting 
m[j]=0.  

After all crossing out, we print out those  j  such that  m[j]=1.
We also print the numbers in a table, replacing the crossed numbers by dots.

erat.c
*****************************************************************************/
      

# include <stdio.h>
# include <stdlib.h>
# include <math.h>




int 
main ( void ) 
{
   
    int i=2,j,k,n=1000, m[1000];
    
    for (j=0; j<=n; j++ )
                  m[j]=1;
    
    printf( " Sieve of Eratosthenes\n\n  List of Primes between 1 and %d\n\n",n);
    
    while( i*i <= n)
    {
         if(m[i] == 1)
         {
             for(j=2*i ; j<=n ; j=j+i)
                   m[j]=0;
         }
         i++;
    }
    
    j = 0;             /*   Print the table.  j  counts the number of primes.  */
    
    for(i=2; i<=n; i++)
    {
          if( m[i] == 1)
                {
                 printf ( "%6d", i );
                 if ( ++j % 10 == 0 )
                            printf ( "\n" );
                 }
    }
    printf ( "\n\n           " );
   
    for ( j = 1; j <= 10; j++ )
          printf ( "%5d", j );
    printf ( "\n" );
    
    for ( k = 0; k <=  100; k = k  +  100 )
    {
          printf ( "         +---------------------------------------------------+\n" );
        for ( j = k; j <= k + 90; j = j + 10 )
        {
    
   
        printf ( "%8d | ", j );
            for ( i = j + 1; i <= j + 10; i++ )
            {
                if (  i  +  j  == 1)
                    printf ( "     " );
                else if ( m[i] == 1 )
                    printf ( "%5d", i );
                else
                    printf ( "  ..." );
            }
            printf ( "|\n" );
         }
    }
    printf("         +---------------------------------------------------+\n");
 
    return EXIT_SUCCESS;
}



/* Milicic

Sieve of Erathostenes 

sieve.c                                 */

#include <stdio.h>

#define TRUE  1
#define FALSE 0
#define MAX 1000

int erath[MAX];

int main() {
int i, j, k;

for (i=0 ; i < MAX ; i++) /* initialize array*/
  erath[i] = FALSE; 

for (j=2 ; j*j < MAX ; j++ )  /* do sieve */
   if (erath[j-1] == FALSE) /* j is not a composite */
      for (k=2 ; j*k <= MAX; k++ )
          erath[k*j-1] = TRUE; /* multiples of j are composites */


for (i = 1 ; i < MAX ; i++) 
   if (erath[i] == FALSE)
      printf("%d is a prime \n", i+1);
}


/*  Treibergs                                                    3-20-6



Random numbers: Monte Carlo  Pi

Note: the range of rand() varies from system to system.
On other systems use  rand() / 2147483648.0;

We initialize the seed for the random number by calling the time of day.

today.c                                                               */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>

int
main(void)
{
    int ik , k=0, n=10000;
    double x,y;

    srand(time(NULL));

    for(ik=1; ik <= n; ik++)
    {

       x = rand() / 32768.0;
       y = rand()  / 32768.0;
  

       if( x*x + y*y < 1.0 )
                       k++;


     }

     printf ( "Monte Carlo Pi = %f  \n", 4.0 * (double)k / n );



return EXIT_SUCCESS;
}


/* Treibergs                3-22-06

Sieve of Eratosthenes.

today.c                                                                      */



# include <stdlib.h>
# include <stdio.h>


int
main(void)
{
       int x[2001], i, j, k = 0, n = 1000;

       for ( i = 1; i <= n; i++ )    /*   Initialize the list   */
             x[i] = 1;

       for ( i = 2; i * i <= n; i++ )   /*    Sieve off composite numbers */
       {
               if( x[i] == 1 )
                     for ( j = 2 * i; j <= n; j = j  +  i )
                             x[j] = 0;
        }
     
       printf ( "Primes Less Than or Equal to %d\n", n);    
                                   /* Print the numbers left after sieving. */
       for ( i = 2;  i <= n ; i++ )
       {
               if ( x[i] == 1 )
               {
                      printf ( "%6d", i );
                      k++;
                      if ( k % 10 == 0)
                                printf ( "\n" );
               }
       }

       printf ( "\nThere were %d primes less than %d\n\n", k, n );



       return EXIT_SUCCESS;
}