From 7ee1619f786a316170228b9bfe23a804f712a6bb Mon Sep 17 00:00:00 2001
From: JereKoskela <jere.koskela@newcastle.ac.uk>
Date: Tue, 6 Feb 2024 20:37:19 +0000
Subject: [PATCH] Numerical stability of beta model

---
 CHANGELOG.md  | 6 ++++++
 lib/msprime.c | 4 ++--
 2 files changed, 8 insertions(+), 2 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index a216b748c..c3fd114c0 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,5 +1,11 @@
 # Changelog
 
+## [1.3.1] - 2024-XX-XX
+
+**Bug fixes**:
+
+- Change tolerance of polynomial approximation in Beta-coalescent acceptance probabilities ({issue}`2256`, {pr}`2257`, {user}`JereKoskela`)
+
 ## [1.3.0] - 2022-12-13
 
 **New features**
diff --git a/lib/msprime.c b/lib/msprime.c
index b3c8cbc3e..3ad619a89 100644
--- a/lib/msprime.c
+++ b/lib/msprime.c
@@ -7056,7 +7056,7 @@ msp_beta_common_ancestor_event(msp_t *self, population_id_t pop_id, label_id_t l
     beta_x = ran_inc_beta(self->rng, 2.0 - alpha, alpha, truncation_point);
 
     /* We calculate the probability of accepting the event */
-    if (beta_x > 1e-9) {
+    if (beta_x > 1e-5) {
         u = (n - 1) * log(1 - beta_x) + log(1 + (n - 1) * beta_x);
         u = exp(log(1 - exp(u)) - 2 * log(beta_x) - gsl_sf_lnchoose(n, 2));
     } else {
@@ -7065,7 +7065,7 @@ msp_beta_common_ancestor_event(msp_t *self, population_id_t pop_id, label_id_t l
         u = 0;
         for (j = 2; j <= n; j += 2) {
             increment = (j - 1) * exp(gsl_sf_lnchoose(n, j) + (j - 2) * log(beta_x));
-            if (increment / u < 1e-12) {
+            if (increment / (u + increment) < 1e-12) {
                 /* We truncate the expansion adaptively once the increment
                  * becomes negligible. */
                 break;