<style>
h2.part{color:#0099B0;}
h3.part{color:#D92424;}
h4.part{color:#005BB0;}
h5.part{color:#FD6F0A;}
h6.part{color:#4400B0;}
</style>

# 2016q3 Homework1 (Compute Pi)
###### tags: `sysprog`
### 開發環境

#### Ubuntu 14.04 LTS

$ lscpu

- CPU: Intel® Core™ i7-2600 CPU @ 3.40GHz × 8
- Mem: 8 GiB
- Cache:
L1d 快取： 32K
L1i 快取： 32K
L2 快取： 256K
L3 快取： 8192K

## 案例分析 - Compute Pi
### Baseline

#### 1. Formula of Pi


double compute_pi_baseline(size_t N)
{
    double pi = 0.0;
    double dt = 1.0 / N;                // dt = (b-a)/N, b = 1, a = 0
    for (size_t i = 0; i < N; i++) {
        double x = (double) i / N;      // x = ti = a+(b-a)*i/N = i/N
        pi += dt / (1.0 + x * x);       // integrate 1/(1+x^2), i = 0....N
    }
    return pi * 4.0;
}

time ./time_test_baseline
N = 400000000 , pi = 3.141593
38.15user 0.00system 0:26.11elapsed 146%CPU (0avgtext+0avgdata 1664maxresident)k
0inputs+0outputs (0major+91minor)pagefaults 0swaps

### OpenMP

double compute_pi_openmp(size_t N, int threads)
{
    double pi = 0.0;
    double dt = 1.0 / N;
    double x;
    #pragma omp parallel num_threads(threads)
    {
        #pragma omp for private(x) reduction(+:pi)
        for (size_t i = 0; i < N; i++) {
            x = (double) i / N;
            pi += dt / (1.0 + x * x);
        }
    }
    return pi * 4.0;
}

#### 2 threads

time ./time_test_openmp_2
N = 400000000 , pi = 3.141593
38.30user 0.00system 0:26.26elapsed 145%CPU (0avgtext+0avgdata 1664maxresident)k
0inputs+0outputs (0major+92minor)pagefaults 0swaps

#### 4 threads

time ./time_test_openmp_4
N = 400000000 , pi = 3.141593
38.08user 0.00system 0:25.97elapsed 146%CPU (0avgtext+0avgdata 1628maxresident)k
0inputs+0outputs (0major+89minor)pagefaults 0swaps

### Monte Carlo
隨機取點(x,y), 判斷是否再圓內

[參考資料](http://mathfaculty.fullerton.edu/mathews/n2003/montecarlopimod.html)


double monte_carlo_pi(size_t N){
	double pi=0.0;
	const unsigned max = 10000;
	unsigned long i = 0;
	double x=0.0, y=0.0;
	unsigned long in_times=0;

	srand(time(NULL)); 
	for(i=0; i<N; i++){
		x=rand()%max;
		y=rand()%max;

		if(x*x + y*y < (double)(max*max))
			in_times++;
	}
	return (4.0 * (double)(in_times)/(double)(N));
}

time ./time_test_monte
N = 400000000 , pi = 3.141593
38.38user 0.00system 0:26.75elapsed 143%CPU (0avgtext+0avgdata 1664maxresident)k
0inputs+0outputs (0major+90minor)pagefaults 0swaps

### Leibniz

#### Leibniz formula for π



double compute_pi_leibniz(size_t N)
{
    double pi = 0.0;
    for(int i=0; i<N; i++)
    {
        int temp = (i%2) ? -1 : 1;
        pi += (double) temp / (2*i+1);        
    }
    return pi * 4.0;
}

##### 執行成果