Monte Carlo methods.

Arrays-Pointers

Sieve of Eratosthenes

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

/***************************************************************************** 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); } }

/***************************************************************************** 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; }

/* 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);## }

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; }