init: v1.0.0

x/ope/hypergeometric.go

@@ -0,0 +1,115 @@
+package ope
+
+import (
+	"math"
+)
+
+// uniformRand returns a random number in [0, 1]
+// len(coins) >= 32
+func uniformRand(coins []byte) float64 {
+	b := int64(0)
+	for i := 0; i < 32; i++ {
+		b <<= 1
+		b |= int64(coins[i])
+	}
+	return float64(b) / float64(uint64(1)<<32-1)
+}
+
+// hypergeometric10 returns a sample in hypergeometric distribution.
+// Pr[X=x] = choose(y, x) * choose(N-y, M-x) / choose(N, M)
+func hypergeometricSmall(M int64, N int64, y int64, coins []byte) int64 {
+	d1 := N - y
+	d2 := min(M, N-M)
+
+	Y := d2
+	K := y
+	for Y > 0.0 && K != 0 {
+		U := uniformRand(coins)
+		Y -= int64(U + float64(Y)/float64(d1+K))
+		K--
+	}
+	Z := d2 - Y
+	if N-M < M {
+		Z = y - Z
+	}
+	return Z
+}
+
+func loggam[T int | int64 | float32 | float64](x T) float64 {
+	return math.Log(math.Gamma(float64(x)))
+}
+
+func hypergeometric(M int64, N int64, y int64, coins []byte) int64 {
+	if y > 10 {
+		return hypergeometricSmall(M, N, y, coins)
+	}
+	return hypergeometricBig(M, N, y, coins)
+}
+
+func hypergeometricBig(M int64, N int64, y int64, coins []byte) int64 {
+	D1 := 1.7155277699214135
+	D2 := 0.8989161620588988
+	// # long mingoodbad, maxgoodbad, popsize, m, d9;
+	// # double d4, d5, d6, d7, d8, d10, d11;
+	// # long Z;
+	// # double T, W, X, Y;
+	good := M
+	bad := N - M
+	sample := y
+
+	mingoodbad := min(good, bad)
+	popsize := N
+	maxgoodbad := max(good, bad)
+	m := min(sample, popsize-sample)
+	d4 := float64(mingoodbad) / float64(popsize)
+	d5 := 1.0 - d4
+	d6 := float64(m)*d4 + 0.5
+	d7 := math.Sqrt(float64(popsize-m)*float64(sample)*d4*d5/float64(popsize-1) + 0.5)
+	d8 := D1*d7 + D2
+	d9 := (m + 1) * (mingoodbad + 1) / (popsize + 2)
+	d10 := loggam(d9+1) + loggam(mingoodbad-d9+1) + loggam(m-d9+1) + loggam(maxgoodbad-m+d9+1)
+	d11 := min(float64(min(m, mingoodbad)+1), math.Floor(d6+16*d7))
+	// # 16 for 16-decimal-digit precision in D1 and D2
+
+	var Z int64
+	for {
+		X := uniformRand(coins)
+		Y := uniformRand(coins)
+		W := d6 + d8*(Y-0.5)/X
+
+		// # fast rejection:
+		if W < 0.0 || W >= d11 {
+			continue
+		}
+
+		Z = int64(math.Floor(W))
+		T := d10 - (loggam(Z+1) + loggam(mingoodbad-Z+1) + loggam(m-Z+1) + loggam(maxgoodbad-m+Z+1))
+
+		// # fast acceptance:
+		if (X*(4.0-X) - 3.0) <= T {
+			break
+		}
+
+		// # fast rejection:
+		if X*(X-T) >= 1 {
+			continue
+		}
+
+		// # acceptance:
+		if 2.0*math.Log(X) <= T {
+			break
+		}
+	}
+
+	// # this is a correction to HRUA* by Ivan Frohne in rv.py
+	if good > bad {
+		Z = m - Z
+	}
+
+	// # another fix from rv.py to allow sample to exceed popsize/2
+	if m < sample {
+		Z = good - Z
+	}
+
+	return Z
+}