178 lines
5.0 KiB
C#
178 lines
5.0 KiB
C#
using System;
|
|
using System.Collections.Generic;
|
|
using System.Globalization;
|
|
using System.Numerics;
|
|
using System.Text;
|
|
|
|
namespace MinesweeperSolver {
|
|
|
|
|
|
public class Binomial {
|
|
|
|
private static readonly double LOG10 = Math.Log(10);
|
|
|
|
private readonly int max;
|
|
private readonly PrimeSieve ps;
|
|
|
|
private readonly BigInteger[,] binomialLookup;
|
|
private readonly int lookupLimit;
|
|
|
|
public Binomial(int max, int lookup) {
|
|
|
|
this.max = max;
|
|
|
|
ps = new PrimeSieve(this.max);
|
|
|
|
if (lookup < 10) {
|
|
lookup = 10;
|
|
}
|
|
this.lookupLimit = lookup;
|
|
|
|
int lookup2 = lookup / 2;
|
|
|
|
binomialLookup = new BigInteger[lookup + 1, lookup2 + 1];
|
|
|
|
for (int total = 1; total <= lookup; total++) {
|
|
for (int choose = 0; choose <= total / 2; choose++) {
|
|
//try {
|
|
binomialLookup[total, choose] = Generate(choose, total);
|
|
//Console.WriteLine("Binomial " + total + " choose " + choose + " is " + binomialLookup[total, choose]);
|
|
//} catch (Exception e) {
|
|
// throw e;
|
|
//}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
public BigInteger Generate(int k, int n) {
|
|
|
|
if (n == 0 && k == 0) {
|
|
return BigInteger.One;
|
|
}
|
|
|
|
if (n < 1) {
|
|
throw new Exception("Binomial: 1 <= n required, but n was " + n);
|
|
}
|
|
|
|
if (0 > k || k > n) {
|
|
throw new Exception("Binomial: 0 <= k and k <= n required, but n was " + n + " and k was " + k);
|
|
}
|
|
|
|
int choose = Math.Min(k, n - k);
|
|
|
|
if (n <= lookupLimit && binomialLookup[n, choose] != 0) { // it is zero when it hasn't been built yet
|
|
return binomialLookup[n, choose];
|
|
} else if (choose < 125) {
|
|
return Combination(choose, n);
|
|
} else if (n <= max) {
|
|
return CombinationLarge(choose, n);
|
|
} else {
|
|
return CombinationApprox(choose, n);
|
|
}
|
|
|
|
}
|
|
|
|
private static BigInteger Combination(int mines, int squares) {
|
|
|
|
BigInteger top = BigInteger.One;
|
|
BigInteger bot = BigInteger.One;
|
|
|
|
int range = Math.Min(mines, squares - mines);
|
|
|
|
// calculate the combination.
|
|
for (int i = 0; i < range; i++) {
|
|
top = top * new BigInteger(squares - i);
|
|
bot = bot * new BigInteger(i + 1);
|
|
}
|
|
|
|
BigInteger result = top / bot;
|
|
|
|
return result;
|
|
|
|
}
|
|
|
|
|
|
private BigInteger CombinationLarge(int k, int n) {
|
|
|
|
if ((k == 0) || (k == n)) return BigInteger.One;
|
|
|
|
int n2 = n / 2;
|
|
|
|
if (k > n2) {
|
|
k = n - k;
|
|
}
|
|
|
|
int nk = n - k;
|
|
|
|
int rootN = (int)Math.Floor(Math.Sqrt(n));
|
|
|
|
BigInteger result = BigInteger.One;
|
|
|
|
foreach (int prime in ps.getPrimesIterable(2, n)) {
|
|
|
|
if (prime > nk) {
|
|
result = result * new BigInteger(prime);
|
|
continue;
|
|
}
|
|
|
|
if (prime > n2) {
|
|
continue;
|
|
}
|
|
|
|
if (prime > rootN) {
|
|
if (n % prime < k % prime) {
|
|
result = result * new BigInteger(prime);
|
|
}
|
|
continue;
|
|
}
|
|
|
|
int r = 0, N = n, K = k, p = 1;
|
|
|
|
while (N > 0) {
|
|
r = (N % prime) < (K % prime + r) ? 1 : 0;
|
|
if (r == 1) {
|
|
p *= prime;
|
|
}
|
|
N /= prime;
|
|
K /= prime;
|
|
}
|
|
if (p > 1) {
|
|
result = result * new BigInteger(p);
|
|
}
|
|
}
|
|
|
|
return result;
|
|
}
|
|
|
|
// use the stirling approximation for factorials to create an approximate combinatorial
|
|
public BigInteger CombinationApprox(int k, int n) {
|
|
|
|
double logComb = (LogFactorialApprox(n) - LogFactorialApprox(k) - LogFactorialApprox(n - k));
|
|
|
|
int power = (int) Math.Floor(logComb);
|
|
|
|
int dp = Math.Min(6, power);
|
|
|
|
//double value = Math.Exp((logComb - power + dp) * LOG10);
|
|
//BigInteger result = new BigInteger(sigDigit) * BigInteger.Pow(new BigInteger(10), power - dp);
|
|
|
|
// find the significant digits
|
|
int sigDigits = (int) Math.Round(Math.Pow(10, logComb - power + dp));
|
|
|
|
String scientific = "" + sigDigits + "E" + (power - dp);
|
|
|
|
BigInteger result = BigInteger.Parse(scientific, NumberStyles.AllowExponent);
|
|
|
|
return result;
|
|
|
|
}
|
|
|
|
|
|
// returns an approximation for Log(n!)
|
|
private double LogFactorialApprox(int n) {
|
|
return n * Math.Log10(n) + 0.5d * Math.Log10(2 * Math.PI * n);
|
|
}
|
|
}
|
|
}
|