1package main
 2
 3import (
 4	"fmt"
 5	"github.com/montanaflynn/stats"
 6	"math"
 7	"math/big"
 8	"math/rand"
 9	"time"
10)
11
12type Trial struct {
13	A       *big.Int
14	B       *big.Int
15	Coprime bool
16}
17
18func main() {
19	var i uint64
20
21	fmt.Printf("n trials, max rand, variance, min, q1, q2, q3, max\n")
22	for i = 1; i < 100; i++ {
23		var max_rand = big.NewInt(120)
24		var num_trials = 2 << i
25		var exp_pis []float64
26		for i := 1; i <= 4096; i++ {
27			exp_pis = append(exp_pis, Experiment(max_rand, num_trials))
28		}
29		max, _ := stats.Max(exp_pis)
30		min, _ := stats.Min(exp_pis)
31		mean, _ := stats.Mean(exp_pis)
32		variance, _ := stats.Variance(exp_pis)
33		q, _ := stats.Quartile(exp_pis)
34		fmt.Printf("%d, %d, %f, %f, %f, %f, %f, %f\n", num_trials, max_rand, variance, min, q.Q1, mean, q.Q3, max)
35	}
36
37	fmt.Printf("n trials, max rand, variance, min, q1, q2, q3, max\n")
38	for i = 1; i < 100; i++ {
39		var max_rand = big.NewInt(2 << i)
40		var num_trials = 500
41		var exp_pis []float64
42		for i := 1; i <= 4096; i++ {
43			exp_pis = append(exp_pis, Experiment(max_rand, num_trials))
44		}
45		max, _ := stats.Max(exp_pis)
46		min, _ := stats.Min(exp_pis)
47		mean, _ := stats.Mean(exp_pis)
48		variance, _ := stats.Variance(exp_pis)
49		q, _ := stats.Quartile(exp_pis)
50		fmt.Printf("%d, %d, %f, %f, %f, %f, %f, %f\n", num_trials, max_rand, variance, min, q.Q1, mean, q.Q3, max)
51	}
52}
53
54func Experiment(max_rand *big.Int, num_trials int) (exp_pi float64) {
55	trials := make([]Trial, num_trials)
56	r := rand.New(rand.NewSource(time.Now().UnixNano()))
57	cc := 0 // coprime count
58	for i, t := range trials {
59		var C *big.Int
60		t.A, t.B, C = big.NewInt(0), big.NewInt(0), big.NewInt(0)
61		trials[i].A = t.A.Rand(r, max_rand)
62		trials[i].B = t.B.Rand(r, max_rand)
63		C = C.GCD(nil, nil, t.A, t.B)
64		if C.Cmp(big.NewInt(1)) == 0 {
65			trials[i].Coprime = true
66			cc += 1
67		}
68	}
69	//for _, t := range trials {
70	//	fmt.Printf("%s, %s, %s, %t\n", t.A.String(), t.B.String(), t.C.String(), t.Coprime)
71	//}
72	return math.Sqrt(6 / (float64(cc) / float64(num_trials)))
73}