```--- a/src/modules/bugs/distributions/DHyper.cc
+++ b/src/modules/bugs/distributions/DHyper.cc
@@ -73,15 +73,13 @@
return true;
}

-/* Calculates the mode of the hypergeometric distribution
-
-   We solve the equation p(x) = p(x-1) for continuous x and then take
-   the floor of x.
-
-   This reduces to solving a quadratic equation in x, by the
-   recurrence relation embodied in the function "rfunction", i.e.
-
-   p(x) = p(x - 1) * rfunction(n1, n2, m1, psi, x)
+/*
+ * Calculates the mode of the hypergeometric distribution
+ *
+ * We solve the equation p(x) = p(x-1) for continuous x. Using the
+ * recurrence relation given in rfunction (see below) this reduces to
+ * a quadratic equation. We then take the floor of x to find the mode
+ * of the distribution.
*/
static int modeCompute(int n1, int n2, int m1, double psi)
{
@@ -106,26 +104,90 @@
}

/*
-   The recurrence relation
-
-   p(x) = p(x - 1) * rfunction(n1, n2, m1, psi, x)
-
-   allows us to calculate the hypergeometric probabilities
-   recursively, avoiding combinatoric problems.
- */
-double rfunction(int n1, int n2, int m1, double psi, int i)
+ * Recurrence relation that allows us to calculate hypergeometric
+ * probabilities without running into combinatoric problems:
+ *
+ * p(x) = p(x - 1) * rfunction(n1, n2, m1, psi, x)
+ */
+static inline double rfunction(int n1, int n2, int m1, double psi, int i)
{
return psi * (n1 - i + 1) * (m1 - i + 1)/(i * (n2 - m1 + i));
}

/**
+ * Calculates the unnormalized density of x. The mode has reference
+ * density value 1.
+ */
+static double density_unnormalized(int x, int n1, int n2, int m1, double psi)
+{
+    int ll = max(0U, m1 - n2);
+    int uu = min(n1, m1);
+
+    if (x < ll || x > uu) return 0;
+
+    // Density at mode has reference value 1
+    int mode = modeCompute(n1, n2, m1, psi);
+    double r = 1;
+
+    if (x > mode) {
+	for (int i = mode + 1; i <= x; ++i) {
+	    r *= rfunction(n1, n2, m1, psi, i);
+	}
+    }
+    else if (x < mode) {
+	for (int i = mode - 1; i >= x; --i) {
+	    r /= rfunction(n1, n2, m1, psi, i + 1);
+	}
+    }
+
+    return r;
+}
+
+/*
+ * Calculates the normalized density of x
+ */
+static double density_normalized(int x, int n1, int n2, int m1, double psi)
+{
+    int ll = max(0U, m1 - n2);
+    int uu = min(n1, m1);
+
+    if (x < ll || x > uu) return 0;
+
+    // Density at mode has reference value 1
+    int mode = modeCompute(n1, n2, m1, psi);
+    double psum = 1, px = 1;
+
+    // Calculate relative density above the mode
+    if (mode < uu) {
+	double r = 1;
+	for (int i = mode + 1; i <= uu; ++i) {
+	    r *= rfunction(n1, n2, m1, psi, i);
+	    if (i == x) px = r;
+	    psum += r;
+	}
+    }
+    // Calculate relative density below the node
+    if (mode > ll) {
+	double r = 1;
+	for (int i = mode - 1; i >= ll; --i) {
+	    r /= rfunction(n1, n2, m1, psi, i + 1);
+	    if (i == x) px = r;
+	    psum += r;
+	}
+    }
+
+    // Normalize
+    return px/psum;
+}
+
+/*
* Returns a vector of normalized probabilities for a hypergeometric
* distribution. If the returned vector p is of length N, then p[0]
* represents the probability of the lowest possible value ll = max(0,
* m1 - n2) and p[N - 1] represents the probability of the largest
- * possible value uu = min(n1, m1)
- */
-static vector<double> density(int n1, int n2, int m1, double psi)
+ * possible value uu = min(n1, m1).
+ */
+static vector<double> density_full(int n1, int n2, int m1, double psi)
{
int ll = max(0U, m1 - n2);
int uu = min(n1, m1);
@@ -135,6 +197,7 @@
// Density at mode has reference value 1
int mode = modeCompute(n1, n2, m1, psi);
p[mode - ll] = 1;
+    double sump = 1;

// Calculate relative density above the mode
if (mode < uu) {
@@ -142,6 +205,7 @@
for (int i = mode + 1; i <= uu; ++i) {
r *= rfunction(n1, n2, m1, psi, i);
p[i - ll] = r;
+	    sump += r;
}
}
// Calculate relative density below the node
@@ -150,14 +214,11 @@
for (int i = mode - 1; i >= ll; --i) {
r /= rfunction(n1, n2, m1, psi, i + 1);
p[i - ll] = r;
+	    sump += r;
}
}

//Normalize
-    double sump = 0;
-    for (int i = 0; i < N; ++i) {
-	sump += p[i];
-    }
for (int i = 0; i < N; ++i) {
p[i] /= sump;
}
@@ -165,48 +226,37 @@
}

/*
- * Sample from a unimodal discrete distribution given a vector of normalized
- * probabilities pi[] and an estimate of the mode. U should be drawn from a
- * uniform(0,1) distribution.
+ * Sample from a discrete unimodal distribution given a vector of
+ * normalized probabilities pi[] and an estimate of the mode. U should
+ * be drawn from a uniform(0,1) distribution. We spend our probability
+ * allocation U as rapidly as possible, starting at the mode and then
+ * stepping out to the left and the right.
+ *
+ * This function is quite general and could be used for other discrete
+ * distributions when we have a good guess of the mode. The function
+ * also works if "mode" is not actually the mode, but is
+ * correspondingly less efficient.
*/
static int sampleWithMode(int mode, vector<double> const &pi, double U)
{
-    //We try to spend our probability allocation U as quickly as
-    //possible, starting with the mode.
-    U -= pi[mode];
-    if (U <= 0) return mode;
-
-    //Then we work our way outwards from the mode, one step at a time.
-    int lower = mode - 1;
-    int upper = mode + 1;
-    while (true) {
-	if (lower < 0) {
-	    for ( ; upper < pi.size() - 1; ++upper) {
-		U -= p[upper];
-		if (U <= 0) break;
-	    }
-	    return upper;
-	}
-	else if (upper >= pi.size()) {
-	    for ( ; lower > 0; --lower) {
-		U -= p[lower];
-		if (U <= 0) break;
-	    }
-	    return lower;
-	}
-	else if (p[upper] > p[lower]) {
-	    U -= p[upper];
-	    if (U <= 0) return upper;
-	    else ++upper;
+    int N = pi.size();
+    int left = mode - 1;
+    int right = mode;
+    while (left >= 0 || right < N) {
+	if (right < N && (left < 0 || p[right] > p[left])) {
+	    U -= pi[right];
+	    if (U <= 0) return right;
+	    else ++right;
}
else {
-	    U -= p[lower];
-	    if (U <= 0) return lower;
-	    else --lower;
-	}
-    }
-
-    return mode; //-Wall
+	    U -= pi[left];
+	    if (U <= 0) return left;
+	    else --left;
+	}
+    }
+
+    //Only reachable with bad arguments
+    return mode;
}

double DHyper::d(double z, PDFType type,
@@ -218,12 +268,12 @@
getParameters(n1, n2, m1, psi, parameters);

int x = static_cast<int>(z);
-    int ll = max(0U, m1 - n2);
-    int uu = min(n1, m1);
-
double den = 0;
-    if (x >= ll && x <= uu) {
-	den = density(n1, n2, m1, psi)[x - ll];
+    if (PDFType == PDF_PRIOR) {
+	den = density_unnormalized(x, n1, n2, m1, psi);
+    }
+    else {
+	den = density_normalized(x, n1, n2, m1, psi);
}

if (give_log) {
@@ -241,7 +291,7 @@
double psi;
getParameters(n1, n2, m1, psi, parameters);

-    int ll = max((int) 0, m1 - n2);
+    int ll = max(0U, m1 - n2);
int uu = min(n1, m1);

double sumpi = 0;
@@ -250,7 +300,7 @@
sumpi = 1;
}
else {
-	    vector<double> pi = density(n1, n2, m1, psi);
+	    vector<double> pi = density_full(n1, n2, m1, psi);
for (int i = ll; i <= x; ++i) {
sumpi += pi[i - ll];
}
@@ -278,7 +328,7 @@
int ll = max((int) 0, m1 - n2);
int uu = min(n1, m1);

-    vector<double> pi = density(n1, n2, m1, psi);
+    vector<double> pi = density_full(n1, n2, m1, psi);

if (log_p)
p = exp(p);
@@ -303,7 +353,7 @@

int ll = max(0U, m1 - n2);
int mode = modeCompute(n1, n2, m1, psi);
-    vector<double> pi = density(n1, n2, m1, psi);
+    vector<double> pi = density_full(n1, n2, m1, psi);
return sampleWithMode(mode - ll, pi, rng->uniform()) + ll;
}

@@ -350,8 +400,8 @@
if (lla < llb || uua > uub)
return JAGS_POSINF;

-    vector<double> da = density(n1a, n2a, m1a, psia);
-    vector<double> db = density(n1b, n2b, m1b, psib);
+    vector<double> da = density_full(n1a, n2a, m1a, psia);
+    vector<double> db = density_full(n1b, n2b, m1b, psib);

double y = 0;
for (int i = lla; i <= uua; ++i) {
```