int main()

        double x, pi, f, dv;

        double sum, sum_pi, sum_pi2;

        int itry, iseed, ntry;

        printf("Input the number of MC trials\n");

        scanf("%d",&ntry);

        srand((unsigned)time((long *)0));

        sum_pi = sum_pi2 = 0.0;

        for (iseed=0; iseed<NSEED; iseed++)

        {

                sum = 0.0;

                for (itry=0; itry<ntry; itry++)

                {

                        x = rand()/(double)RAND_MAX;

                        f = 4.0/(1.0 + x*x);

                        sum += f;

                }

                pi = sum/ntry;

                sum_pi += pi;

                sum_pi2 += pi*pi;

        }

        dv = sqrt((sum_pi2/NSEED) - (sum_pi/NSEED)*(sum_pi/NSEED));

        printf("MC estimate for PI = %f +- %e \n", sum_pi/NSEED, dv);

        return 0;