From 8ff55698125efe22c05ada9a0bc127dc316f3383 Mon Sep 17 00:00:00 2001
From: "Steinar H. Gunderson" <sesse@debian.org>
Date: Sun, 21 Oct 2007 14:06:08 +0200
Subject: [PATCH] Minor cleanups and optimizations.

---
 foosrank.cpp | 72 +++++++++++++++++++---------------------------------
 1 file changed, 26 insertions(+), 46 deletions(-)

diff --git a/foosrank.cpp b/foosrank.cpp
index cc66076..45c08cc 100644
--- a/foosrank.cpp
+++ b/foosrank.cpp
@@ -16,9 +16,9 @@ static const double rating_constant = 455.0;
 
 using namespace std;
 
-double prob_score(int k, double a, double rd);
-double prob_score_real(int k, double a, double binomial, double rd_norm);
-double prodai(int k, double a);
+double prob_score(int k, int a, double rd);
+double prob_score_real(int k, int a, double binomial, double rd_norm);
+double prodai(int k, int a);
 double fac(int x);
 
 
@@ -38,23 +38,40 @@ double fac(int x);
 // Glicko/Bradley-Terry assumption that a player rated 400 points over
 // his/her opponent will win with a probability of 10/11 =~ 0.90909. 
 //
-double prob_score(int k, double a, double rd)
+double prob_score(int k, int a, double rd)
 {
 	return prob_score_real(k, a, prodai(k, a) / fac(k-1), rd/rating_constant);
 }
 
+// computes x^a, probably more efficiently than pow(x, a) (but requires that a
+// is n unsigned integer)
+double intpow(double x, unsigned a)
+{
+	double result = 1.0;
+
+	while (a > 0) {
+		if (a & 1) {
+			result *= x;
+		}
+		a >>= 1;
+		x *= x;
+	}
+
+	return result;
+}
+
 // Same, but takes in binomial(a+k-1, k-1) as an argument in
 // addition to a. Faster if you already have that precomputed, and assumes rd
 // is already divided by 455.
-double prob_score_real(int k, double a, double binomial, double rd_norm)
+double prob_score_real(int k, int a, double binomial, double rd_norm)
 {
-	double nom = binomial * pow(2.0, rd_norm * a); 
-	double denom = pow(1.0 + pow(2.0, rd_norm), k+a);
+	double nom = binomial * intpow(pow(2.0, rd_norm), a); 
+	double denom = intpow(1.0 + pow(2.0, rd_norm), k+a);
 	return nom/denom;
 }
 
 // Calculates Product(a+i, i=1..k-1) (see above).
-double prodai(int k, double a)
+double prodai(int k, int a)
 {
 	double prod = 1.0;
 	for (int i = 1; i < k; ++i)
@@ -70,48 +87,11 @@ double fac(int x)
 	return prod;
 }
 
-// 
-// Computes the integral
-//
-//   +inf
-//    /
-//    |
-//    | ProbScore[a] (r1-r2) Gaussian[mu2, sigma2] (r2) dr2
-//    |
-//   /
-// -inf
-//
-// For practical reasons, -inf and +inf are replaced by 0 and 3000, which
-// is reasonable in the this context.
-//
-// The Gaussian is not normalized.
-//
-// Set the last parameter to 1.0 if player 1 won, or -1.0 if player 2 won.
-// In the latter case, ProbScore will be given (r2-r1) instead of (r1-r2).
-//
-class ProbScoreEvaluator {
-private:
-	int k;
-	double a;
-	double binomial_precompute, r1, mu2, sigma2, winfac;
-
-public:
-	ProbScoreEvaluator(int k, double a, double binomial_precompute, double r1, double mu2, double sigma2, double winfac)
-		: k(k), a(a), binomial_precompute(binomial_precompute), r1(r1), mu2(mu2), sigma2(sigma2), winfac(winfac) {}
-	inline double operator() (double x) const
-	{
-		double probscore = prob_score_real(k, a, binomial_precompute, (r1 - x)*winfac);
-		double z = (x - mu2)/sigma2;
-		double gaussian = exp(-(z*z/2.0));
-		return probscore * gaussian;
-	}
-};
-
 void convolve(int size)
 {
 }
 
-void compute_opponent_rating_pdf(int k, double a, double mu2, double sigma2, double winfac, vector<pair<double, double> > &result)
+void compute_opponent_rating_pdf(int k, int a, double mu2, double sigma2, double winfac, vector<pair<double, double> > &result)
 {
 	double binomial_precompute = prodai(k, a) / fac(k-1);
 	winfac /= rating_constant;
-- 
2.39.2