// This library is part of PLINK 2.0, copyright (C) 2005-2026 Shaun Purcell,
// Christopher Chang.
//
// This library is free software: you can redistribute it and/or modify it
// under the terms of the GNU Lesser General Public License as published by the
// Free Software Foundation; either version 3 of the License, or (at your
// option) any later version.
//
// This library is distributed in the hope that it will be useful, but WITHOUT
// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
// FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public License
// for more details.
//
// You should have received a copy of the GNU Lesser General Public License
// along with this library.  If not, see <http://www.gnu.org/licenses/>.

#include "plink2_stats.h"

#include <assert.h>
#include <math.h>
#include <stdlib.h>  // exit()

#include "plink2_float.h"
#include "plink2_highprec.h"

#ifdef __cplusplus
namespace plink2 {
#endif

// Thread-unsafe portions of plink_stats.c have been replaced, mostly by code
// derived from boost/math/special_functions/gamma.hpp and
// boost/math/special_functions/detail/igamma_inverse.hpp in Boost 1.60 and
// 1.84 (Maddock et al.).  The derived portions are subject to the following
// license:
//
// *****
// Boost Software License - Version 1.0 - August 17th, 2003
//
// Permission is hereby granted, free of charge, to any person or organization
// obtaining a copy of the software and accompanying documentation covered by
// this license (the "Software") to use, reproduce, display, distribute,
// execute, and transmit the Software, and to prepare derivative works of the
// Software, and to permit third-parties to whom the Software is furnished to
// do so, all subject to the following:
//
// The copyright notices in the Software and this entire statement, including
// the above license grant, this restriction and the following disclaimer,
// must be included in all copies of the Software, in whole or in part, and
// all derivative works of the Software, unless such copies or derivative
// works are solely in the form of machine-executable object code generated by
// a source language processor.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE, TITLE AND NON-INFRINGEMENT. IN NO EVENT
// SHALL THE COPYRIGHT HOLDERS OR ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE
// FOR ANY DAMAGES OR OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE,
// ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
// DEALINGS IN THE SOFTWARE.
// *****

// ***** thread-safe ChisqToP *****
// port of Boost 1.60 implementation, float precision

static const double kLogMinValue = -708.0;
static const double kLogMaxValue = 709.0;

static const double kLentzFpmin = 1.0e-30;

// may move some of these tables to plink2_float
static const double kFactorialRecips[34] = {
  1.0,
  1.0,
  0.5,
  0.16666666666666666,
  0.041666666666666664,
  0.008333333333333333,
  0.001388888888888889,
  0.0001984126984126984,
  2.48015873015873e-05,
  2.7557319223985893e-06,
  2.755731922398589e-07,
  2.505210838544172e-08,
  2.08767569878681e-09,
  1.6059043836821613e-10,
  1.1470745597729725e-11,
  7.647163731819816e-13,
  4.779477332387385e-14,
  2.8114572543455206e-15,
  1.5619206968586225e-16,
  8.22063524662433e-18,
  4.110317623312165e-19,
  1.9572941063391263e-20,
  8.896791392450574e-22,
  3.8681701706306835e-23,
  1.6117375710961184e-24,
  6.446950284384474e-26,
  2.4795962632247972e-27,
  9.183689863795546e-29,
  3.279889237069838e-30,
  1.1309962886447718e-31,
  3.7699876288159054e-33,
  1.2161250415535181e-34,
  3.800390754854744e-36,
  1.151633562077195e-37
};

static const double kSmallRecips[30] = {
  0.0,  // could make this nan, though that's annoying for C++03
  1.0 / 1,
  1.0 / 2,
  1.0 / 3,
  1.0 / 4,
  1.0 / 5,
  1.0 / 6,
  1.0 / 7,
  1.0 / 8,
  1.0 / 9,
  1.0 / 10,
  1.0 / 11,
  1.0 / 12,
  1.0 / 13,
  1.0 / 14,
  1.0 / 15,
  1.0 / 16,
  1.0 / 17,
  1.0 / 18,
  1.0 / 19,
  1.0 / 20,
  1.0 / 21,
  1.0 / 22,
  1.0 / 23,
  1.0 / 24,
  1.0 / 25,
  1.0 / 26,
  1.0 / 27,
  1.0 / 28,
  1.0 / 29
};

double finite_gamma_q(uint32_t aa, double xx, double* p_derivative) {
  // a is a positive integer < 30; max(0.6, a-1) < x < kLogMaxValue
  // (e^{-x})(1 + x + x^2/2 + x^3/3! + x^4/4! + ... + x^{a-1}/(a-1)!)
  const double ee = exp(-xx);
  if (ee == 0.0) {
    return 0;
  }
  double sum = ee;
  double term = sum;
  for (uint32_t nn = 1; nn != aa; ++nn) {
    term *= xx * kSmallRecips[nn];
    sum += term;
  }
  if (p_derivative) {
    *p_derivative = ee * pow(xx, u31tod(aa)) * kFactorialRecips[aa - 1];
  }
  return sum;
}

static const double kSqrtPi = 1.7724538509055159;

double lower_gamma_series(double aa, double zz, double init_value) {
  // evaluate init_value + 1 + (z/(a+1)) + (z^2 / ((a+1)(a+2))) +
  //   (z^3 / ((a+1)(a+2)(a+3))) + ...
  // z shouldn't be larger than (a+170), otherwise overflow is possible
  double result = 1;
  double total = init_value + 1;
  do {
    aa += 1.0;
    result *= zz / aa;
    total += result;
  } while (fabs(result) > (kBigEpsilon * kBigEpsilon));
  return total;
}

// ~6 digits of precision is appropriate for p-value computations
static const double kContinuedFractionEpsilon = 3.0e-7;

double upper_gamma_fraction(double a1, double z1) {
  // No overflow issues, this tends to be close to 1.

  // evaluate a_1 / (b_1 + (a_2 / (b_2 + (a_3 / (b_3 + ...)))))
  // see Boost continued_fraction_a(), upper_incomplete_gamma_fract
  double cur_b = z1 - a1 + 3;

  double hh = cur_b;
  const double a0 = a1 - 1.0;
  if (fabs(hh) < kLentzFpmin) {
    hh = kLentzFpmin;
  }
  double cc = hh;
  double dd = 0.0;
  for (double kk = 2.0; kk <= 100.0; kk += 1.0) {
    const double cur_a = kk * (a1 - kk);
    cur_b += 2.0;
    dd = cur_a * dd + cur_b;
    if (fabs(dd) < kLentzFpmin) {
      dd = kLentzFpmin;
    }
    cc = cur_b + cur_a / cc;
    if (fabs(cc) < kLentzFpmin) {
      cc = kLentzFpmin;
    }
    dd = 1.0 / dd;
    const double delta = cc * dd;
    hh *= delta;
    if (fabs(delta - 1.0) < kContinuedFractionEpsilon) {
      break;
    }
  }
  // const double cont_frac = a0 / hh;
  // return 1 / (z1 - a1 + 1 + cont_frac);
  return hh / prefer_fma(hh, z1 - a1 + 1, a0);
}

// from Numerical Recipes in Fortran 77: The Art of Scientific Computing, via
// Wikipedia
// maximal error of 1.2e-7
double erfc_fast(double zz) {
  const double tt = 1.0 / (1.0 + 0.5 * zz);
  const double tt2 = tt * tt;
  const double tt3 = tt * tt2;
  const double tt4 = tt2 * tt2;
  // - zz * zz - 1.26551223
  // + tt * 1.00002368
  // + tt^2 * 0.37409196
  // + tt^3 * 0.09678418
  // + tt^4 * (-0.18628806)
  // + tt^5 * 0.27886807
  // + tt^6 * (-1.13520398)
  // + tt^7 * 1.48851587
  // + tt^8 * (-0.82215223)
  // + tt^9 * 0.17087277
  const double rem0 = prefer_fma(prefer_fma(-0.82215223, tt4, -0.18628806), tt4, -prefer_fma(zz, zz, 1.26551223));
  const double rem1 = prefer_fma(prefer_fma(0.17087277, tt4, 0.27886807), tt4, 1.00002368) * tt;
  const double rem2 = prefer_fma(-1.13520398, tt4, 0.37409196) * tt2;
  const double rem3 = prefer_fma(1.48851587, tt4, 0.09678418) * tt3;
  return tt * exp(rem0 + rem1 + (rem2 + rem3));
}

static const double kSmallHalfRecips[30] = {
  1.0 / 0.5,
  1.0 / 1.5,
  1.0 / 2.5,
  1.0 / 3.5,
  1.0 / 4.5,
  1.0 / 5.5,
  1.0 / 6.5,
  1.0 / 7.5,
  1.0 / 8.5,
  1.0 / 9.5,
  1.0 / 10.5,
  1.0 / 11.5,
  1.0 / 12.5,
  1.0 / 13.5,
  1.0 / 14.5,
  1.0 / 15.5,
  1.0 / 16.5,
  1.0 / 17.5,
  1.0 / 18.5,
  1.0 / 19.5,
  1.0 / 20.5,
  1.0 / 21.5,
  1.0 / 22.5,
  1.0 / 23.5,
  1.0 / 24.5,
  1.0 / 25.5,
  1.0 / 26.5,
  1.0 / 27.5,
  1.0 / 28.5,
  1.0 / 29.5
};

// This can underflow.  Use finite_half_gamma_q2_ln when that's a problem.
double finite_half_gamma_q2(uint32_t a_minus_half, double xx, double* p_derivative) {
  // a is in {0.5, 1.5, ..., 29.5}; max(0.2, a-1) < x < kLogMaxValue
  const double sqrt_x = sqrt(xx);
  double ee = erfc_fast(sqrt_x);
  if ((ee != 0) && a_minus_half) {
    double term = sqrt_x * exp(-xx) * (2.0 / kSqrtPi);
    double sum = term;
    for (uint32_t n_minus_half = 1; n_minus_half != a_minus_half; ++n_minus_half) {
      term *= xx * kSmallHalfRecips[n_minus_half];
      sum += term;
    }
    ee += sum;
    if (p_derivative) {
      *p_derivative = 0;
    }
  } else if (p_derivative) {
    *p_derivative = sqrt_x * exp(-xx) * (1.0 / kSqrtPi);
  }
  return ee;
}

static const double kLanczosFloatSumNumer[6] = {8706.3495925490091, 8523.650341121874, 3338.0292194764235, 653.64249942940087, 63.999518449381870, 2.5066282746310063};
static const double kLanczosFloatSumDenom[6] = {0, 24, 50, 35, 10, 1};
static const double kLanczosFloatSumExpgNumer[6] = {32.812445410297834, 32.123889414443320, 12.580347294552161, 2.4634444783532414, 0.2412010548258800, 0.0094469677045392};

// this depends on the polynomial coefficients above
static const double kLanczosFloatG = 5.581;

double lanczos_sum_f(double zz) {
  double s1;
  double s2;
  if (zz <= 1) {
    s1 = kLanczosFloatSumNumer[5];
    s2 = kLanczosFloatSumDenom[5];
    for (int32_t ii = 4; ii >= 0; --ii) {
      s1 = prefer_fma(s1, zz, kLanczosFloatSumNumer[S_CAST(uint32_t, ii)]);
      s2 = prefer_fma(s2, zz, kLanczosFloatSumDenom[S_CAST(uint32_t, ii)]);
    }
  } else {
    zz = 1 / zz;
    s1 = kLanczosFloatSumNumer[0];
    s2 = kLanczosFloatSumDenom[0];
    for (uint32_t uii = 1; uii != 6; ++uii) {
      s1 = prefer_fma(s1, zz, kLanczosFloatSumNumer[uii]);
      s2 = prefer_fma(s2, zz, kLanczosFloatSumDenom[uii]);
    }
  }
  return s1 / s2;
}

double lanczos_sum_f_expg_scaled_imp(double zz, double* s2_ptr) {
  double s1;
  double s2;
  if (zz <= 1) {
    s1 = kLanczosFloatSumExpgNumer[5];
    s2 = kLanczosFloatSumDenom[5];
    for (int32_t ii = 4; ii >= 0; --ii) {
      s1 = prefer_fma(s1, zz, kLanczosFloatSumExpgNumer[S_CAST(uint32_t, ii)]);
      s2 = prefer_fma(s2, zz, kLanczosFloatSumDenom[S_CAST(uint32_t, ii)]);
    }
  } else {
    zz = 1 / zz;
    s1 = kLanczosFloatSumExpgNumer[0];
    s2 = kLanczosFloatSumDenom[0];
    for (uint32_t uii = 1; uii != 6; ++uii) {
      s1 = prefer_fma(s1, zz, kLanczosFloatSumExpgNumer[uii]);
      s2 = prefer_fma(s2, zz, kLanczosFloatSumDenom[uii]);
    }
  }
  *s2_ptr = s2;
  return s1;
}

double lanczos_sum_f_expg_scaled(double zz) {
  double s2;
  double s1 = lanczos_sum_f_expg_scaled_imp(zz, &s2);
  return s1 / s2;
}

double lanczos_sum_f_expg_scaled_recip(double zz) {
  double s2;
  double s1 = lanczos_sum_f_expg_scaled_imp(zz, &s2);
  return s2 / s1;
}

double log1pmx(double xx) {
  // log(1+x) - x
  // assumes abs(xx) < 0.95
  const double aa = fabs(xx);
  if (aa < (kBigEpsilon / kSqrt2)) { // 2^{-21.5}
    return -xx * xx * 0.5;
  }
  double kk = 1.0;  // skip first term of usual log(1+x) series
  const double m_mult = -xx;
  double m_prod = xx;
  double total = 0.0;
  double rr;
  do {
    m_prod *= m_mult;
    kk += 1.0;
    rr = m_prod / kk;
    total += rr;
    // todo: tune these epsilons, but let's wait until we know all of the
    // callers of these functions
  } while (fabs(rr) > (kBigEpsilon * kBigEpsilon));
  return total;
}

// compute (z^a)(e^{-z})/tgamma(a)
double regularized_gamma_prefix(double aa, double zz) {
  // assumes a == 0.5 if a < 1.  assumes z > 0.
  // we are fine with float-level precision, so lanczos_n=6,
  // kLanczosFloatG=5.581
  if (aa < 1) {
    return sqrt(zz) * exp(-zz) * (1.0 / kSqrtPi);
  }
  const double agh = aa + kLanczosFloatG - 0.5;
  const double agh_recip = 1.0 / agh;
  const double dd = ((zz - aa) - (kLanczosFloatG - 0.5)) * agh_recip;
  double prefix;
  if ((fabs(dd * dd * aa) <= 100) && (aa > 150)) {
    // abs(dd) < sqrt(2/3) < 0.95
    prefix = zz * (0.5 - kLanczosFloatG) * agh_recip + aa * log1pmx(dd);
    prefix = exp(prefix);
  } else {
    const double alz = aa * log(zz * agh_recip);
    const double amz = aa - zz;
    const double cur_minv = MINV(alz, amz);
    if ((cur_minv <= kLogMinValue) || (MAXV(alz, amz) >= kLogMaxValue)) {
      const double amza = amz / aa;
      double sq;
      if ((cur_minv > 2 * kLogMinValue) && (MAXV(alz, amz) < 2 * kLogMaxValue)) {
        // need to structure this to avoid overflow
        sq = pow(zz * agh_recip, aa * 0.5) * exp(amz * 0.5);
        prefix = sq * sq;
      } else if ((cur_minv > 4 * kLogMinValue) && (MAXV(alz, amz) < 4 * kLogMaxValue) && (zz > aa)) {
        sq = pow(zz * agh_recip, aa * 0.25) * exp(amz * 0.25);
        prefix = sq * sq;
        prefix *= prefix;
      } else if ((amza > kLogMinValue) && (amza < kLogMaxValue)) {
        prefix = pow((zz * exp(amza)) * agh_recip, aa);
      } else {
        prefix = exp(alz + amz);
      }
    } else {
      prefix = pow(zz * agh_recip, aa) * exp(amz);
    }
  }
  prefix *= sqrt(agh * kRecipE) * lanczos_sum_f_expg_scaled_recip(aa);
  return prefix;
}

static const double kTemmeC0[7] = {-0.333333333, 0.0833333333, -0.0148148148, 0.00115740741, 0.000352733686, -0.000178755144, 0.391926318e-4};
static const double kTemmeC1[5] = {-0.00185185185, -0.00347222222, 0.00264550265, -0.000990226337, 0.000205761317};
static const double kTemmeC2[3] = {0.00413359788, -0.00268132716, 0.000771604938};

double igamma_temme_large(double aa, double xx) {
  // 24-bit precision is fine
  const double sigma = (xx - aa) / aa;
  // abs(sigma) < 0.4
  const double phi = -log1pmx(sigma);
  const double sqrt_a = sqrt(aa);
  const double sqrt_phi = sqrt(phi);
  const double yy = aa * phi;
  double zz = kSqrt2 * sqrt_phi;
  if (xx < aa) {
    zz = -zz;
  }
  double workspace[3];
  // todo: benchmark some other ways of organizing this computation.  there are
  // three independent calculation chains, which I'd expect to be decent?
  workspace[0] = prefer_fma(prefer_fma(prefer_fma(prefer_fma(prefer_fma(prefer_fma(kTemmeC0[6], zz, kTemmeC0[5]), zz, kTemmeC0[4]), zz, kTemmeC0[3]), zz, kTemmeC0[2]), zz, kTemmeC0[1]), zz, kTemmeC0[0]);
  workspace[1] = prefer_fma(prefer_fma(prefer_fma(prefer_fma(kTemmeC1[4], zz, kTemmeC1[3]), zz, kTemmeC1[2]), zz, kTemmeC1[1]), zz, kTemmeC1[0]);
  workspace[2] = prefer_fma(prefer_fma(kTemmeC2[2], zz, kTemmeC2[1]), zz, kTemmeC2[0]);
  const double a_recip = 1 / aa;
  double result = prefer_fma(prefer_fma(workspace[2], a_recip, workspace[1]), a_recip, workspace[0]);
  result *= exp(-yy) / ((kSqrt2 * kSqrtPi) * sqrt_a);
  if (xx < aa) {
    result = -result;
  }
  result += erfc_fast(sqrt_a * sqrt_phi) * 0.5;
  return result;
}

double gamma_incomplete_imp2(uint32_t df, double xx, uint32_t invert, double* p_derivative) {
  // normalized = true
  assert(df);
  assert(xx >= 0.0);
  const double aa = u31tod(df) * 0.5;
  const uint32_t is_small_a = (df < 60) && (aa <= xx + 1) && (xx < kLogMaxValue);
  uint32_t is_int = 0;
  uint32_t is_half_int = 0;
  if (is_small_a) {
    is_half_int = df % 2;
    is_int = !is_half_int;
  }
  uint32_t eval_method;
  if (is_int && (xx > 0.6)) {
    invert = !invert;
    eval_method = 0;
  } else if (is_half_int && (xx > 0.2)) {
    invert = !invert;
    eval_method = 1;
  } else if (xx < kSmallEpsilon) {
    // avoid computing log(0)
    // don't need more precision here, 6 digits is enough
    // invert always == 1
    assert(!p_derivative);
    return 1.0;
  } else if (xx < 1.1) {
    // (simplified from original Boost code)
    eval_method = 2;
  } else {
    const double x_minus_a = xx - aa;
    uint32_t use_temme = 0;
    if (aa > 20) {
      // sigma = abs((x - a) / a);
      // igamma_temme_large() assumes abs(sigma) < 0.95
      if (aa > 200) {
        // abs(sigma) < sqrt(20 / a) < 0.316...
        use_temme = (20 * aa > x_minus_a * x_minus_a);
      } else {
        // abs(sigma) < 0.4
        const double sigma_times_a = fabs(x_minus_a);
        use_temme = (sigma_times_a < 0.4 * aa);
      }
    }
    if (use_temme) {
      eval_method = 5;
    } else {
      // x - (1 / (3 * x)) < a
      // x * x - (1/3) < a * x
      // x * x - a * x < 1/3
      // x * (x - a) < 1/3
      if (xx * x_minus_a < (1.0 / 3.0)) {
        eval_method = 2;
      } else {
        eval_method = 4;
        invert = !invert;
      }
    }
  }
  double result;
  switch(eval_method) {
  case 0:
    result = finite_gamma_q(df / 2, xx, p_derivative);
    break;
  case 1:
    // previously used erfc, but that was up to ~3x as slow as dcdflib (e.g.
    // ChisqToP(2.706, 1) case).
    result = finite_half_gamma_q2(df / 2, xx, p_derivative);
    if (p_derivative && (*p_derivative == 0)) {
      *p_derivative = regularized_gamma_prefix(aa, xx);
    }
    break;
  case 2:
    result = regularized_gamma_prefix(aa, xx);
    if (p_derivative) {
      *p_derivative = result;
    }
    if (result != 0) {
      // uint32_t optimized_invert = 0;
      double init_value = 0;
      if (invert) {
        init_value = -aa / result;
        // optimized_invert = 1;
      }
      result *= lower_gamma_series(aa, xx, init_value) / aa;
      // if (optimized_invert) {
      if (invert) {
        invert = 0;
        result = -result;
      }
    }
    break;
    // case 3:
    // only when df=1 and 0.449 < x < 1.1, so no overflow/underflow issues.
    // ...wait a sec, we never reach here at all since eval_method is always
    // 1 there.
  case 4:
    result = regularized_gamma_prefix(aa, xx);
    if (p_derivative) {
      *p_derivative = result;
    }
    if (result != 0) {
      result *= upper_gamma_fraction(aa, xx);
    }
    break;
  case 5:
    result = igamma_temme_large(aa, xx);
    if (xx >= aa) {
      invert = !invert;
    }
    if (p_derivative) {
      *p_derivative = regularized_gamma_prefix(aa, xx);
    }
  }
  if (result > 1) {
    result = 1;
  }
  if (invert) {
    result = 1 - result;
  }
  if (p_derivative) {
    if ((xx < 1) && (DBL_MAX * xx < (*p_derivative))) {
      *p_derivative = DBL_MAX / 2;  // overflow; do we really need this?
    } else {
      *p_derivative /= xx;
    }
  }
  return result;
}

double ChisqToP(double chisq, uint32_t df) {
  // todo: figure out when we were depending on this to return -9, and decide
  // how to handle those situations now
  return gamma_incomplete_imp2(df, chisq * 0.5, 1, nullptr);
}

// ***** end thread-safe ChisqToP *****


// ***** ChisqToLnP *****

double finite_gamma_q_ln(uint32_t aa, double xx) {
  // a is a positive integer < 30; max(0.6, a-1) < x < 1e10
  // (e^{-x})(1 + x + x^2/2 + x^3/3! + x^4/4! + ... + x^{a-1}/(a-1)!)
  //
  // logarithm:
  // log(1 + x + ... + x^{a-1}/(a-1)!) - x
  // no overflow or underflow danger for main term thanks to bounds
  double sum = 1.0;
  double term = 1.0;
  for (uint32_t nn = 1; nn != aa; ++nn) {
    // division is slow enough that the lookup table speeds up this function by
    // >3x
    term *= xx * kSmallRecips[nn];
    sum += term;
  }
  return log(sum) - xx;
}

double erfc_fast2(double zz, double* tau_ln_plus_z2_ptr) {
  // probable todo: deduplicate with erfc_fast()
  const double tt = 1.0 / (1.0 + 0.5 * zz);
  const double tt2 = tt * tt;
  const double tt3 = tt * tt2;
  const double tt4 = tt2 * tt2;
  const double rem0 = prefer_fma(prefer_fma(-0.82215223, tt4, -0.18628806), tt4, -1.26551223);
  const double rem1 = prefer_fma(prefer_fma(0.17087277, tt4, 0.27886807), tt4, 1.00002368) * tt;
  const double rem2 = prefer_fma(-1.13520398, tt4, 0.37409196) * tt2;
  const double rem3 = prefer_fma(1.48851587, tt4, 0.09678418) * tt3;
  *tau_ln_plus_z2_ptr = rem0 + rem1 + (rem2 + rem3);
  return tt;
}

double finite_half_gamma_q2_ln(uint32_t a_minus_half, double xx) {
  // a is in {0.5, 1.5, ..., 29.5}; max(0.2, a-1) < x < 1e10
  const double sqrt_x = sqrt(xx);
  double tau_ln_plus_x;
  double tt = erfc_fast2(sqrt_x, &tau_ln_plus_x);
  if (!a_minus_half) {
    return log(tt) + tau_ln_plus_x - xx;
  }
  double ee = tt * exp(tau_ln_plus_x);
  // pre-multiply by e^x to avoid underflow
  double term = sqrt_x * (2.0 / kSqrtPi);
  double sum = term;
  for (uint32_t n_minus_half = 1; n_minus_half != a_minus_half; ++n_minus_half) {
    term *= xx * kSmallHalfRecips[n_minus_half];
    sum += term;
  }
  // tau_ln + xx guaranteed to be small
  ee += sum;
  return log(ee) - xx;
}

// compute -log((z^a)(e^{-z})/tgamma(a))
double regularized_gamma_prefix_ln(double aa, double zz) {
  // assumes a == 0.5 if a < 1.  assumes z > 0.
  // we are fine with float-level precision, so lanczos_n=6,
  // kLanczosFloatG=5.581
  if (aa < 1) {
    return -zz + 0.5 * log(zz) - kLnSqrtPi;
  }
  const double agh = aa + kLanczosFloatG - 0.5;
  const double agh_recip = 1.0 / agh;
  const double dd = ((zz - aa) - (kLanczosFloatG - 0.5)) * agh_recip;
  double prefix_ln;
  if ((fabs(dd * dd * aa) <= 100) && (aa > 150)) {
    // abs(dd) < sqrt(2/3) < 0.95
    prefix_ln = zz * (0.5 - kLanczosFloatG) * agh_recip + aa * log1pmx(dd);
  } else {
    prefix_ln = prefer_fma(aa, log(zz * agh_recip), aa - zz);
  }
  const double scaled_recip = lanczos_sum_f_expg_scaled_recip(aa);
  return prefix_ln + 0.5 * log(agh * kRecipE * scaled_recip * scaled_recip);
}

// does not guarantee return value <= 0 for now; caller must do that.
double gamma_incomplete_imp2_ln(uint32_t df, double xx) {
  // normalized = true, invert = false
  assert(df);
  assert(xx >= 0.0);
  const double aa = u31tod(df) * 0.5;
  const uint32_t is_small_a = (df < 60) && (aa <= xx + 1) && (xx < 1e10);
  uint32_t is_int = 0;
  uint32_t is_half_int = 0;
  if (is_small_a) {
    is_half_int = df % 2;
    is_int = !is_half_int;
  }
  uint32_t eval_method;
  if (is_int && (xx > 0.6)) {
    eval_method = 0;
  } else if (is_half_int && (xx > 0.2)) {
    eval_method = 1;
  } else if (xx < kSmallEpsilon) {
    // avoid computing log(0)
    // don't need more precision here, 6 digits is enough
    return 0.0;
  } else if (xx < 1.1) {
    eval_method = 2;
  } else {
    const double x_minus_a = xx - aa;
    uint32_t use_temme = 0;
    if (aa > 20) {
      // sigma = abs((x - a) / a);
      // igamma_temme_large() assumes abs(sigma) < 0.95
      if (aa > 200) {
        // abs(sigma) < sqrt(20 / a) < 0.316...
        use_temme = (20 * aa > x_minus_a * x_minus_a);
      } else {
        // abs(sigma) < 0.4
        const double sigma_times_a = fabs(x_minus_a);
        use_temme = (sigma_times_a < 0.4 * aa);
      }
    }
    if (use_temme) {
      eval_method = 5;
    } else {
      // x - (1 / (3 * x)) < a
      // x * x - (1/3) < a * x
      // x * x - a * x < 1/3
      // x * (x - a) < 1/3
      if (xx * x_minus_a < (1.0 / 3.0)) {
        eval_method = 2;
      } else {
        eval_method = 4;
      }
    }
  }
  switch(eval_method) {
  case 0:
    return finite_gamma_q_ln(df / 2, xx);
  case 1:
    return finite_half_gamma_q2_ln(df / 2, xx);
  case 2:
    {
      const double result_ln = regularized_gamma_prefix_ln(aa, xx);
      if (result_ln < kLogMinValue + 22) {
        // init_value overflows.  Not a big deal, this just ends up getting
        // inverted to pval=1.
        // (+22 since aa could theoretically be as large as 2^31.  Todo: find
        // the smallest result_ln value that could result in a nonzero value
        // being returned.)
        return 0.0;
      }
      const double init_value = -aa * exp(-result_ln);
      const double multiplier = -lower_gamma_series(aa, xx, init_value) / aa;
      return result_ln + log(multiplier);
    }
  case 4:
    {
      const double result1_ln = regularized_gamma_prefix_ln(aa, xx);
      const double result2_ln = log(upper_gamma_fraction(aa, xx));
      return result1_ln + result2_ln;
    }
  case 5:
  default:  // silence compiler warning
    {
      // aa large, fabs(xx - aa) relatively small, so no overflow/underflow
      // issues.
      double result = igamma_temme_large(aa, xx);
      if (xx < aa) {
        result = 1.0 - result;
      }
      return log(result);
    }
  }
}

double ChisqToLnP(double chisq, uint32_t df) {
  return MINV(gamma_incomplete_imp2_ln(df, chisq * 0.5), 0.0);
}
// ***** end ChisqToLnP *****


// ***** thread-safe PToChisq *****
// port of Boost 1.60 implementation

double find_inverse_gamma2(uint32_t df, double pp, double qq, uint32_t* has_10_digits_ptr) {
  // currently assumes *has_10_digits_ptr initialized to zero
  if (df == 2) {
    return -log(qq);
  }
  if (df == 1) {
    // g == tgamma(0.5) == sqrt(pi)
    const double bb = qq * kSqrtPi;
    if (bb >= 0.45) {
      // b * q > 1e-8, q > 1e-5 guaranteed
      // u = pow(p * g * a, 1/a)
      //   = pow(p * g * 0.5, 2)
      //   = p * p * g * g * 0.25
      //   = p * p * pi * 0.25
      const double uu = pp * pp * (0.25 * kPi);
      return uu / prefer_fma(uu, -1.0 / 1.5, 1);
    } else {
      const double yy = -log(bb);
      if (bb > 0.1) {
        const double uu = yy - 0.5 * log(yy);
        if (bb > 0.15) {
          return (yy - 0.5 * log(uu) - log(1 + 0.5 / uu));
        }
        return (yy - 0.5 * log(uu) - log(prefer_fma(uu + 5, uu, 3.75) / prefer_fma(uu + 4.5, uu, 2)));
      } else {
        const double c1 = -0.5 * log(yy);
        const double c1_2 = c1 * c1;
        const double c1_3 = c1_2 * c1;
        const double c1_4 = c1_2 * c1_2;
        // a_2 = 0.25
        // a_3 = 0.125

        const double c2 = -0.5 * (1 + c1);
        const double c3 = 0.25 * c1_2 + prefer_fma(0.75, c1, 0.875);
        const double c4 = prefer_fma(c1_3, -1.0 / 6.0, -26.75 / 12.0) - prefer_fma(0.875, c1_2, 1.875 * c1);
        const double c5 = prefer_fma(5.75 / 6.0, c1_3, 0.125 * c1_4) + prefer_fma(3.625, c1_2, prefer_fma(7.75, c1, 83.0625 / 12.0));

        const double y_recip = 1.0 / yy;
        const double y_recip_2 = y_recip * y_recip;
        const double y_recip_3 = y_recip_2 * y_recip;
        const double y_recip_4 = y_recip_2 * y_recip_2;
        if (bb < 1e-28) {
          *has_10_digits_ptr = 1;
        }
        // possible todo: benchmark other ways of organizing this calculation
        return prefer_fma(c2, y_recip, yy) + prefer_fma(c3, y_recip_2, c1) + prefer_fma(c5, y_recip_4, c4 * y_recip_3);
      }
    }
  }
  // not implemented yet
  assert(0);
  exit(1);
  return 0;
}

double gamma_p_inv_imp2(uint32_t df, double qq) {
  assert(df);
  assert(qq > 0.0);
  if (qq >= 1.0 - kSmallEpsilon) {
    return 0;
  }
  double pp = 1.0 - qq;
  uint32_t has_10_digits = 0;
  double guess = find_inverse_gamma2(df, pp, qq, &has_10_digits);
  if (has_10_digits) {
    return guess;
  }
  double min_guess = kSmallEpsilon;
  double max_guess = DBL_MAX;
  if (guess < kSmallEpsilon) {
    guess = kSmallEpsilon;
  }
  // halley iteration, digits == 24, lower == kSmallEpsilon
  // see second_order_root_finder in boost/math/tools/roots.hpp
  const uint32_t invert = (pp > 0.9);
  if (invert) {
    pp = qq;
  }
  const double a_minus_1 = 0.5 * S_CAST(double, S_CAST(int32_t, df) - 2);
  const double factor = 1.1920928955078125e-07;  // 2^{-23}
  double result = guess;
  double delta = 10000000;
  double delta1 = delta;
  uint32_t out_of_bounds_sentry = 0;
  do {
    double delta2 = delta1;
    delta1 = delta;

    // see gamma_p_inverse_func in
    // boost/math/special_functions/detail/igamma_inverse.hpp
    double f1;
    const double ff = gamma_incomplete_imp2(df, result, invert, &f1);
    const double div = (a_minus_1 - result) / result;
    double f2 = f1;
    if ((fabs(div) > 1) && (DBL_MAX / fabs(div) < f2)) {
      // overflow
      f2 = -DBL_MAX / 2;
    } else {
      f2 *= div;
    }
    if (invert) {
      f1 = -f1;
      f2 = -f2;
    }
    const double f0 = ff - pp;
    if (f0 == 0) {
      break;
    }
    assert(f1 != 0);  // shouldn't be possible, function is monotonic
    delta = f0 / f1;
    if (f2 != 0) {
      // delta = Stepper::step(result, f0, f1, f2);
      const double denom = 2 * f0;
      const double numer = prefer_fma(-f0, f2 / f1, 2 * f1);
      if ((fabs(numer) >= 1) || (fabs(denom) < fabs(numer) * DBL_MAX)) {
        const double halley_step = denom / numer;
        if (halley_step / delta < 0) {
          if (fabs(delta) > 2 * fabs(guess)) {
            delta = ((delta < 0)? -1 : 1) * 2 * fabs(guess);
          }
        } else {
          delta = halley_step;
        }
      }
    }
    double convergence = fabs(delta / delta2);
    if ((convergence > 0.8) && (convergence < 2)) {
      delta = (delta > 0)? (0.5 * (result - min_guess)) : (0.5 * (result - max_guess));
      if (fabs(delta) > result) {
        delta = ((delta > 0)? 1 : -1) * result;
      }
      // delta2 = delta * 3;
    }
    guess = result;
    result -= delta;
    // do we actually need this?
    if (result < min_guess) {
      double diff = ((fabs(min_guess) < 1) && (fabs(result) > 1) && ((DBL_MAX / fabs(result)) < fabs(min_guess)))? 1000 : (result / min_guess);
      if (fabs(diff) < 1) {
        diff = 1 / diff;
      }
      if ((!out_of_bounds_sentry) && (diff > 0) && (diff < 3)) {
        delta = 0.99 * (guess - min_guess);
        result = guess - delta;
        out_of_bounds_sentry = 1;
      } else {
        delta = (guess - min_guess) * 0.5;
        result = guess - delta;
        if ((result == min_guess) || (result == max_guess)) {
          break;
        }
      }
    } else if (result > max_guess) {
      double diff = ((fabs(max_guess) < 1) && (fabs(result) > 1) && ((DBL_MAX / fabs(result)) < fabs(max_guess)))? 1000 : (result / max_guess);
      if (fabs(diff) < 1) {
        diff = 1 / diff;
      }
      if ((!out_of_bounds_sentry) && (diff > 0) && (diff < 3)) {
        delta = 0.99 * (guess - max_guess);
        result = guess - delta;
        out_of_bounds_sentry = 1;
      } else {
        delta = (guess - max_guess) * 0.5;
        result = guess - delta;
        if ((result == min_guess) || (result == max_guess)) {
          break;
        }
      }
    }
    if (delta > 0) {
      max_guess = guess;
    } else {
      min_guess = guess;
    }
  } while (fabs(result * factor) < fabs(delta));
  return result;
}

// double PToChisq(double pval, uint32_t df) {
//   // only need this to handle df=1, 2, 4 for now
//   return gamma_p_inv_imp2(df, pval) * 2;
// }

// ***** end thread-safe PToChisq *****

double LnPToChisq(double ln_pval) {
  if (ln_pval > -65.04474754675797) {
    return gamma_p_inv_imp2(1, exp(ln_pval)) * 2;
  }
  // (bb < 1e-28) case in find_inverse_gamma2()
  // possible todo: deduplicate code
  const double yy = kLnSqrtPi - ln_pval;
  const double c1 = -0.5 * log(yy);
  const double c1_2 = c1 * c1;
  const double c1_3 = c1_2 * c1;
  const double c1_4 = c1_2 * c1_2;
  // a_2 = 0.25
  // a_3 = 0.125

  const double c2 = -0.5 * (1 + c1);
  const double c3 = 0.25 * c1_2 + prefer_fma(0.75, c1, 0.875);
  const double c4 = prefer_fma(c1_3, -1.0 / 6.0, -26.75 / 12.0) - prefer_fma(0.875, c1_2, 1.875 * c1);
  const double c5 = prefer_fma(5.75 / 6.0, c1_3, 0.125 * c1_4) + prefer_fma(3.625, c1_2, prefer_fma(7.75, c1, 83.0625 / 12.0));

  const double y_recip = 1.0 / yy;
  const double y_recip_2 = y_recip * y_recip;
  const double y_recip_3 = y_recip_2 * y_recip;
  const double y_recip_4 = y_recip_2 * y_recip_2;
  return 2 * (prefer_fma(c2, y_recip, yy) + prefer_fma(c3, y_recip_2, c1) + prefer_fma(c5, y_recip_4, c4 * y_recip_3));
}


// ***** thread-safe TstatToLnP *****

// Returns log(1 - x^y).
// Assumes x in (0, 1), y positive.
double neg_powm1_imp_ln(double xx, double yy) {
  const double ll = yy * log(xx);
  if (ll > -1.0) {
    // For tiny |l|, we may lose all precision when exp(l) evaluates to 1.
    return log(-expm1(ll));
  }
  // For larger |l|, it's the log step that may throw away precision.
  if (ll < kLogMinValue) {
    return 0.0;
  }
  return log1p(-exp(ll));
}

double binomial_ccdf_ln(uint32_t nn, uint32_t kk, double xx, double yy, uint32_t inv) {
  // x^n + (n choose 1)x^{n-1}y + ... + (n choose (n-k-1))x^{k+1}y^{n-k-1}
  // This is currently just designed to work with ibeta_imp2_ln().  So we
  // assume x and y positive, y = 1 - x, (n-k) < 40, (n-k) <= (k+1)(y/x).

  // Thanks to the (n-k) <= (k+1)(y/x) condition, it's always reasonable to
  // start from the rightmost term normalized to 1, and sum leftward; then we
  // multiply by the rightmost term at the end.
  const double n_d = u31tod(nn);
  const double k_plus1_d = u31tod(kk + 1);
  const double x_div_y = xx / yy;
  double multiplier = 1.0;
  double cur_term = 1.0;
  // need to shift to avoid overflow for n > 1 billion, (n-k) close to 40
  double shifted_inv_binom_coeff = k2p50 * (1LL << 33);
  double i_d = n_d - k_plus1_d;
  for (uint32_t uii = nn - kk - 1; uii; --uii) {
    //   (n choose (i - 1)) / (n choose i)
    //   i!(n-i)! / (i-1)!(n-i+1)!
    // = i / (n-i+1)
    if (cur_term < k2m53) {
      // Only need to finish computing shifted_inv_binom_coeff: multiply by
      // i!, divide by n!/(n-i)!.
      // Since we support n up to 2^31, n!/(n-i)! may overflow for i>33.
      for (; uii > 33; --uii) {
        shifted_inv_binom_coeff *= i_d;
        i_d -= 1.0;
        shifted_inv_binom_coeff /= n_d - i_d;
      }
      double denom = kFactorialRecips[uii];
      double n_minus_i = n_d - i_d;
      for (; uii; --uii) {
        n_minus_i += 1.0;
        denom *= n_minus_i;
      }
      shifted_inv_binom_coeff /= denom;
      break;
    }
    double cur_binom_ratio = i_d;
    i_d -= 1.0;
    cur_binom_ratio /= n_d - i_d;
    shifted_inv_binom_coeff *= cur_binom_ratio;
    cur_term *= cur_binom_ratio * x_div_y;
    multiplier += cur_term;
  }
  // tried taking log(multiplier / shifted_inv_binom_coeff), but that was
  // slower on my Mac?
  const double result_ln = log(multiplier) + k_plus1_d * log(xx) - log(shifted_inv_binom_coeff) + prefer_fma(n_d - k_plus1_d, log(yy), 83 * kLn2);
  if (!inv) {
    return result_ln;
  }
  // The nightmare scenario of tiny result_ln shouldn't happen thanks to the
  // (n-k) <= (k+1)(y/x) condition.
  return log1p(-exp(result_ln));
}

double ibeta_series_ln(double aa, double bb, double xx, uint32_t inv) {
  // BPSER in DiDonato and Morris.
  // ibeta(a, b, x) ~= gamma(a+b) / (gamma(a) * gamma(b))
  //                 * x^a / a
  //                 * (1 + a\sum_j=1^N [(1-b)(2-b)...(j-b)x^j / (j!(a+j))])
  // Currently always called with x <= 0.7, and either b*x <= 0.7 or b<40, so
  // guaranteed to converge at reasonable speed.

  // normalized always true

  const double cc = aa + bb;

  const double agh = aa + kLanczosFloatG - 0.5;
  const double bgh = bb + kLanczosFloatG - 0.5;
  const double cgh = cc + kLanczosFloatG - 0.5;
  double numer_a;
  double denom_a = lanczos_sum_f_expg_scaled_imp(aa, &numer_a);
  double numer_b;
  double denom_b = lanczos_sum_f_expg_scaled_imp(bb, &numer_b);
  double denom_c;
  double numer_c = lanczos_sum_f_expg_scaled_imp(cc, &denom_c);
  double result = (numer_a * numer_b * numer_c) / (denom_a * denom_b * denom_c);
  double l1 = log(cgh / bgh) * (bb - 0.5);
  double l2 = log(xx * cgh / agh) * aa;
  double result_ln = log(result * result * agh) * 0.5 + l1 + l2 - 0.5;

  double series_sum = 0.0;
  double term = 1.0;
  double apn = aa;
  double poch = 1.0 - bb;
  double nn = 1.0;
  while (1) {
    double rr = term / apn;
    series_sum += rr;
    if (fabs(rr) <= fabs(series_sum * kEpsilon)) {
      result_ln += log(series_sum);
      if (!inv) {
        return result_ln;
      }
      // assume for now that we wouldn't be inverting into the tiny result_ln
      // case.
      return log1p(-exp(result_ln));
    }
    apn += 1.0;
    term *= poch * xx / nn;
    nn += 1.0;
    poch += 1.0;
  }
}

double ln_sum(double aa, double bb) {
  if (aa > bb) {
    swap_f64(&aa, &bb);
  }
  const double diff = bb - aa;
  if ((aa == -DBL_MAX) || (diff >= kLn2 * 53)) {
    return bb;
  }
  return bb + log(1 + exp(-diff));
}

double tgamma_delta_ratio(double zz, double delta) {
  // gamma(z) / gamma(z + delta)
  // We only call this with delta = 0.5 or 1, so no overflow issues.
  // zz >= 15 for now.
  if (delta == 1.0) {
    // Trivial case.
    return 1.0 / zz;
  }
  assert(delta == 0.5);
  // Can skip z < epsilon and z + delta == z branches for now.
  double zgh = zz + kLanczosFloatG - 0.5;
  // Also skip fabs(delta) >= 10 branch for now.
  double result = exp((0.5 - zz) * log1p(delta / zgh));
  result *= lanczos_sum_f(zz) / lanczos_sum_f(zz + delta);
  // exploit forced delta == 0.5
  result *= sqrt(kE / (zgh + delta));
  return result;
}

double beta_small_b_large_a_series_ln(double aa, double bb, double xx, double yy, double s0_ln, uint32_t inv) {
  // BGRAT in DiDonato and Morris.

  // normalized always true, mult always 1
  // a >= 15, b == 0.5 or 1, and ibeta_imp was patched to ensure x > 0.7.
  double bm1 = bb - 1.0;
  double tt = aa + bm1 * 0.5;
  double lx;
  if (yy < 0.35) {
    lx = log1p(-yy);
  } else {
    lx = log(xx);
  }
  double uu = -tt * lx;
  double hh_ln = regularized_gamma_prefix_ln(bb, uu);
  double prefix_ln = prefer_fma(-bb, log(tt), hh_ln) - log(tgamma_delta_ratio(aa, bb));
  // validated up to this point

  double pp[15]; // ~8-15 digit accuracy
  pp[0] = 1.0;

  double jj;
  {
    double dxx;
    if (bb == 0.5) {
      // bugfix (17 Jun 2019): original expression could underflow
      // jj = finite_half_gamma_q2(0, uu, nullptr);
      dxx = finite_half_gamma_q2_ln(0, uu);
    } else {
      assert(bb == 1.0);
      dxx = -uu;
    }
    // Underflow is harmless here, jj becomes nonzero for the other terms in
    // the series in that case.
    jj = exp(dxx - hh_ln);
  }
  double sum = jj; // patch in s0 and prefix at the end
  uint32_t tnp1 = 1;
  double lx2 = lx * 0.5;
  lx2 *= lx2;
  double lxp = 1.0;
  double t4_recip = 0.25 / (tt * tt);
  double b2n = bb;
  for (uint32_t nn = 1; nn < 15; ++nn) {
    tnp1 += 2;
    double new_pn = 0.0;
    uint32_t tmp1 = 3;
    const double nn_d = u31tod(nn);
    // when this type of loop is a bottleneck, easiest win should be going up
    // to 2-8 accumulators
    for (uint32_t mm = 1; mm < nn; ++mm) {
      double mbn = prefer_fma(u31tod(mm), bb, -nn_d);
      new_pn += mbn * pp[nn - mm] * kFactorialRecips[tmp1];
      tmp1 += 2;
    }
    new_pn /= nn_d;
    new_pn += bm1 * kFactorialRecips[tnp1];
    pp[nn] = new_pn;

    jj = (b2n * (b2n + 1) * jj + (uu + b2n + 1) * lxp) * t4_recip;
    lxp *= lx2;
    b2n += 2;

    double rr = new_pn * jj;
    sum += rr;
    if (rr > 1) {
      if (fabs(rr) < fabs(kEpsilon * sum)) {
        break;
      }
    } else {
      if (fabs(rr * (1.0 / kEpsilon)) < fabs(sum)) {
        break;
      }
    }
  }
  double result_ln = ln_sum(s0_ln, prefix_ln + log(sum));
  if (!inv) {
    return result_ln;
  }
  return log1p(-exp(result_ln));
}

double ibeta_power_terms_ln(double aa, double bb, double xx, double yy) {
  // returns log((x^a)(y^b) / Beta(a,b))
  //
  // normalized always true
  // prefix always 1
  double cc = aa + bb;

  const double agh = aa + kLanczosFloatG - 0.5;
  const double bgh = bb + kLanczosFloatG - 0.5;
  const double cgh = cc + kLanczosFloatG - 0.5;
  double numer_a;
  double denom_a = lanczos_sum_f_expg_scaled_imp(aa, &numer_a);
  double numer_b;
  double denom_b = lanczos_sum_f_expg_scaled_imp(bb, &numer_b);
  double denom_c;
  double numer_c = lanczos_sum_f_expg_scaled_imp(cc, &denom_c);
  double result = (numer_a * numer_b * numer_c) / (denom_a * denom_b * denom_c);
  result *= sqrt(agh * bgh * kRecipE / cgh);
  double result_ln = log(result);
  double l1 = (xx * bb - yy * agh) / agh;
  double l2 = (yy * aa - xx * bgh) / bgh;
  // Since we're returning log(result) rather than the original result (thus,
  // no intermediate overflow/underflow problems), and we only need 6 digits of
  // precision, we shouldn't need any of the numerous cases in the Boost code.
  return prefer_fma(aa, log1p(l1), result_ln) + bb * log1p(l2);
}

double ibeta_a_step_ln(double aa, double bb, double xx, double yy, uint32_t kk) {
  double prefix_ln = ibeta_power_terms_ln(aa, bb, xx, yy) - log(aa);
  double sum = 1.0;
  double term = 1.0;
  const double a_plus_b = aa + bb;
  const double k_minus_1pt5 = u31tod(kk) - 1.5;
  for (double i_d = 0.0; i_d < k_minus_1pt5; ) {
    term *= (a_plus_b + i_d) * xx;
    i_d += 1.0;
    term /= aa + i_d;
    sum += term;
  }
  return prefix_ln + log(sum);
}

double ibeta_fraction2_ln(double aa, double bb, double xx, double yy, uint32_t inv) {
  // normalized always true

  // todo: original DiDonato and Morris paper notes that "x must also be a
  // sufficient distance from p when a > 100"; check if we have a problem there

  double result_ln = ibeta_power_terms_ln(aa, bb, xx, yy);

  // see Boost continued_fraction_b()
  const double ay_minus_bx_plus1 = aa * yy - prefer_fma(bb, xx, -1.0);
  double cc = (aa * ay_minus_bx_plus1) / (aa + 1.0);
  if (fabs(cc) < kLentzFpmin) {
    cc = kLentzFpmin;
  }
  const double a_plus_b = aa + bb;
  const double x2 = xx * xx;
  result_ln -= log(cc);
  double dd = 0.0;
  double mm = 0.0;
  while (1) {
    double cur_a = (aa + mm) * (a_plus_b + mm);
    mm += 1.0;
    cur_a *= mm * (bb - mm) * x2;
    double denom = aa + 2 * mm - 1.0;
    cur_a /= denom * denom;
    double cur_b = prefer_fma(mm * (bb - mm), xx / (aa + 2 * mm - 1.0), mm);
    cur_b += ((aa + mm) * prefer_fma(mm, 2.0 - xx, ay_minus_bx_plus1)) / (aa + 2 * mm + 1.0);
    dd = cur_b + cur_a * dd;
    if (fabs(dd) < kLentzFpmin) {
      dd = kLentzFpmin;
    }
    cc = cur_b + cur_a / cc;
    if (fabs(cc) < kLentzFpmin) {
      cc = kLentzFpmin;
    }
    dd = 1.0 / dd;
    const double delta = cc * dd;
    result_ln -= log(delta);
    if (fabs(delta - 1.0) <= kContinuedFractionEpsilon) {
      return inv? log1p(-exp(result_ln)) : result_ln;
    }
  }
}

// the '2' in imp2 refers to df1 and df2 arguments instead of a and b
double ibeta_imp2_ln(uint32_t df1, uint32_t df2, double xx, uint32_t inv) {
  // In addition to Boost beta.hpp and its dependencies, see DiDonato and
  // Morris's original paper at
  // https://apps.dtic.mil/dtic/tr/fulltext/u2/a210118.pdf .

  // normalized always true
  //
  // assume a >= 0, b >= 0, a+b > 0, x in [0, 1], a and b are multiples of 0.5
  // for now
  // in practice, x always <= 0.5 right now
  if (df1 == 0) {
    return inv? -DBL_MAX : 0.0;
  }
  if (df2 == 0) {
    return inv? 0.0 : -DBL_MAX;
  }
  if (xx == 0.0) {
    return inv? 0.0 : -DBL_MAX;
  }
  if (xx == 1.0) {
    return inv? -DBL_MAX : 0.0;
  }
  double yy = 1.0 - xx;
  if ((df1 == 1) && (df2 == 1)) {
    return log(asin(sqrt(inv? yy : xx)) * (2.0 / kPi));
  }
  if (df1 == 2) {
    df1 = df2;
    df2 = 2;

    swap_f64(&xx, &yy);

    inv = !inv;
  }
  double aa = u31tod(df1) * 0.5;
  if (df2 == 2) {
    if (df1 == 2) {
      return log(inv? yy : xx);
    }
    double ln_pp;
    if (yy < 0.5) {
      //       pp = inv? -expm1(aa * log1p(-yy)) : exp(aa * log1p(-yy))
      // -> ln_pp = inv? log(-exp(aa * log1p(-yy)) + 1.0) : (aa * log1p(-yy))
      //          = inv? log1p(-exp(aa * log1p(-yy))) : (aa * log1p(-yy))
      const double inner_term = aa * log1p(-yy);
      ln_pp = inv? log1p(-exp(inner_term)) : inner_term;
    } else {
      //       pp = inv? -powm1(xx, aa) : pow(xx, aa);
      // -> ln_pp = inv? log(-powm1(xx, aa)) : (aa * log(xx))
      ln_pp = inv? neg_powm1_imp_ln(xx, aa) : (aa * log(xx));
    }
    return ln_pp;
  }

  double bb = u31tod(df2) * 0.5;
  if ((df1 == 1) || (df2 == 1)) {
    if (xx > 0.5) {
      swap_u32(&df1, &df2);
      swap_f64(&aa, &bb);
      swap_f64(&xx, &yy);
      inv = !inv;
    }
    // can ignore max(a, b) <= 1 branch
    if ((df2 == 1) || ((xx < 0.1) && (bb * xx <= 0.49))) {
      return ibeta_series_ln(aa, bb, xx, inv);
    }
    // a/b swapped, x/y swapped, inv swapped
    if (xx >= 0.3) {
      return ibeta_series_ln(bb, aa, yy, !inv);
    }
    if (bb >= 15) {
      return beta_small_b_large_a_series_ln(bb, aa, yy, xx, -DBL_MAX, !inv);
    }
    double fract_ln = ibeta_a_step_ln(bb, aa, yy, xx, 20);
    return beta_small_b_large_a_series_ln(bb + 20.0, aa, yy, xx, fract_ln, !inv);
  }
  double lambda;
  if (aa < bb) {
    lambda = prefer_fma(aa + bb, -xx, aa);
  } else {
    lambda = prefer_fma(aa + bb, yy, -bb);
  }
  if (lambda < 0.0) {
    swap_u32(&df1, &df2);
    swap_f64(&aa, &bb);
    swap_f64(&xx, &yy);
    inv = !inv;
  }

  // a > 1, b > 1 guaranteed if we get here.
  if (df2 >= 80) {
    return ibeta_fraction2_ln(aa, bb, xx, yy, inv);
  }
  if ((!(df1 % 2)) && (!(df2 % 2)) && (yy != 1.0)) {
    const uint32_t a_int = df1 / 2;
    const uint32_t b_int = df2 / 2;
    const uint32_t kk = a_int - 1;
    const uint32_t nn = b_int + kk;
    return binomial_ccdf_ln(nn, kk, xx, yy, inv);
  }
  // Changed from b*x <= 0.7 since BGRAT has problems with small x, while BPSER
  // can handle this larger b*x case since b<40 is guaranteed.
  if (xx <= 0.7) {
    return ibeta_series_ln(aa, bb, xx, inv);
  }
  const uint32_t nn = (df2 - 1) / 2;
  const double bbar = bb - u31tod(nn);
  double fract_ln = ibeta_a_step_ln(bbar, aa, yy, xx, nn);
  if (df1 > 30) {
    return beta_small_b_large_a_series_ln(aa, bbar, xx, yy, fract_ln, inv);
  }
  fract_ln = ln_sum(fract_ln, ibeta_a_step_ln(aa, bbar, xx, yy, 20));
  return beta_small_b_large_a_series_ln(aa + 20.0, bbar, xx, yy, fract_ln, inv);
}

double TstatToLnP(double tt, uint32_t df) {
  const double df_d = u31tod(df);
  const double t2 = tt * tt;
  if (df_d > 2 * t2) {
    const double zz = t2 / (df_d + t2);
    return ibeta_imp2_ln(1, df, zz, 1);
  }
  const double zz = df_d / (df_d + t2);
  return ibeta_imp2_ln(df, 1, zz, 0);
}
// ***** end TstatToLnP *****

// ***** FstatToLnP *****

double FstatToLnP(double ff, uint32_t df1, uint32_t df2) {
  // See cdf() for complemented2_type in boost/math/distributions/fisher_f.hpp.
  const double df1_d = u31tod(df1);
  const double df2_d = u31tod(df2);
  const double v1x = df1_d * ff;
  if (v1x > df2_d) {
    return ibeta_imp2_ln(df2, df1, df2_d / (df2_d + v1x), 0);
  }
  return ibeta_imp2_ln(df1, df2, v1x / (df2_d + v1x), 1);
}

// ***** end FstatToLnP *****


// Inverse normal distribution
// (todo: check if boost implementation is better)

// Lower tail quantile for standard normal distribution function.
//
// This function returns an approximation of the inverse cumulative
// standard normal distribution function.  I.e., given P, it returns
// an approximation to the X satisfying P = Pr{Z <= X} where Z is a
// random variable from the standard normal distribution.
//
// The algorithm uses a minimax approximation by rational functions
// and the result has a relative error whose absolute value is less
// than 1.15e-9.
//
// Author:      Peter J. Acklam
// Time-stamp:  2002-06-09 18:45:44 +0200
// E-mail:      jacklam@math.uio.no
// WWW URL:     http://www.math.uio.no/~jacklam
//
// C implementation adapted from Peter's Perl version

// Coefficients in rational approximations.

static const double kIvnA[] =
  {
    -3.969683028665376e+01,
    2.209460984245205e+02,
    -2.759285104469687e+02,
    1.383577518672690e+02,
    -3.066479806614716e+01,
     2.506628277459239e+00
  };

static const double kIvnB[] =
  {
    -5.447609879822406e+01,
    1.615858368580409e+02,
    -1.556989798598866e+02,
    6.680131188771972e+01,
    -1.328068155288572e+01
  };

static const double kIvnC[] =
  {
    -7.784894002430293e-03,
    -3.223964580411365e-01,
    -2.400758277161838e+00,
    -2.549732539343734e+00,
    4.374664141464968e+00,
     2.938163982698783e+00
  };

static const double kIvnD[] =
  {
    7.784695709041462e-03,
    3.224671290700398e-01,
    2.445134137142996e+00,
    3.754408661907416e+00
  };

static const double kIvnLow = 0.02425;
static const double kIvnHigh = 0.97575;

double QuantileToZscore(double pval) {
  // This was named ltqnorm() in PLINK 1.x.
  // assumes 0 < pval < 1
  if ((pval < kIvnLow) || (pval > kIvnHigh)) {
    // Rational approximation for lower/upper region
    const double q2 = -2 * log((pval < kIvnLow)? pval : (1 - pval));
    const double q = sqrt(q2);
    // return (((((kIvnC[0]*q+kIvnC[1])*q+kIvnC[2])*q+kIvnC[3])*q+kIvnC[4])*q+kIvnC[5]) /
    //   ((((kIvnD[0]*q+kIvnD[1])*q+kIvnD[2])*q+kIvnD[3])*q+1);
    const double numer_rem0 = prefer_fma(prefer_fma(kIvnC[1], q2, kIvnC[3]), q2, kIvnC[5]);
    const double numer_rem1 = q * prefer_fma(prefer_fma(kIvnC[0], q2, kIvnC[2]), q2, kIvnC[4]);
    const double denom_rem0 = prefer_fma(prefer_fma(kIvnD[0], q2, kIvnD[2]), q2, 1);
    const double denom_rem1 = q * prefer_fma(kIvnD[1], q2, kIvnD[3]);
    const double frac = (numer_rem0 + numer_rem1) / (denom_rem0 + denom_rem1);
    return (pval < kIvnLow)? frac : -frac;
  }
  // Rational approximation for central region
  const double q = pval - 0.5;
  const double q2 = q*q;
  const double q4 = q2*q2;
  const double numer_rem1 = prefer_fma(prefer_fma(kIvnA[1], q4, kIvnA[3]), q4, kIvnA[5]);
  const double numer_rem3 = q2 * prefer_fma(prefer_fma(kIvnA[0], q4, kIvnA[2]), q4, kIvnA[4]);
  const double denom_rem0 = prefer_fma(prefer_fma(kIvnB[1], q4, kIvnB[3]), q4, 1);
  const double denom_rem2 = q2 * prefer_fma(prefer_fma(kIvnB[0], q4, kIvnB[2]), q4, kIvnB[4]);
  return q * (numer_rem1 + numer_rem3) / (denom_rem0 + denom_rem2);
}

double Lfact(double xx) {
  if (xx < 256) {
    return _ddr_ln_fact[S_CAST(int32_t, xx)].x[0];
  }
  // Since we're summing a fixed number of terms, start with the smaller terms
  // to reduce rounding error.
  const double invn = 1.0 / xx;
  const double invn2 = invn * invn;
  const double small_term_sum = prefer_fma(invn, prefer_fma(invn2, 1.0 / -360.0, 1.0 / 12.0), kLnSqrt2Pi - xx);
  return log(xx) * (xx + 0.5) + small_term_sum;
}

// - starting_lnprobv_ddr is expected to be initialized to
//     log(2^obs_hets / (obs_hets! obs_hom1! obs_hom2!))
//   or have x[0] initialized to DBL_MAX to indicate that the calculation
//   hasn't happened.  In the latter case, it may be set to the former value
//   if that is needed in the calculation.
//
// - *cmp_resultp is set to positive value if hets := obs_hets + 2*hom_decr has
//   higher likelihood than hets := obs_hets, 0 if identical likelihood, and
//   negative value if lower likelihood.
//
// - Error is returned iff malloc fails.
BoolErr HweCompare(uint32_t obs_hets, uint32_t obs_hom1, uint32_t obs_hom2, int32_t hom_decr, dd_real* starting_lnprobv_ddr_ptr, intptr_t* cmp_resultp, double* dbl_ptr) {
  // From e.g. the Wigginton paper, P(N_{AB}=n_{AB} | N, n_A) is
  //
  //      2^{n_{AB}} N! n_A! n_B!
  //   -----------------------------
  //   n_{AA}! n_{AB}! n_{BB}! (2N)!
  //
  // Thus, P(N_{AB}=obs_hets + 2*hom_decr) / P(N_{AB}=obs_hets) is
  //
  //       obs_hets!         obs_hom1! * obs_hom2!      2j
  //   ---------------- * -------------------------- * 2
  //   (obs_hets + 2j)!   (obs_hom1-j)!(obs_hom2-j)!
  //
  // where j=hom_decr.
  uint32_t numer_factorial_args[3];
  numer_factorial_args[0] = obs_hets;
  numer_factorial_args[1] = obs_hom1;
  numer_factorial_args[2] = obs_hom2;
  uint32_t denom_factorial_args[3];
  denom_factorial_args[0] = obs_hets + 2 * hom_decr;
  denom_factorial_args[1] = obs_hom1 - hom_decr;
  denom_factorial_args[2] = obs_hom2 - hom_decr;

  mp_limb_t* gmp_wkspace = nullptr;
  uintptr_t gmp_wkspace_limb_ct = 0;
  BoolErr reterr = CompareFactorialProducts(3, hom_decr * 2LL, obs_hets, numer_factorial_args, denom_factorial_args, starting_lnprobv_ddr_ptr, &gmp_wkspace, &gmp_wkspace_limb_ct, cmp_resultp, dbl_ptr);
  free_cond(gmp_wkspace);
  return reterr;
}

// obs_hets + obs_hom1 + obs_hom2 assumed to be <2^31.
BoolErr HweLnP(int32_t obs_hets, int32_t obs_hom1, int32_t obs_hom2, int32_t midp, double* resultp) {
  // This function implements an exact SNP test of Hardy-Weinberg
  // Equilibrium as described in Wigginton, JE, Cutler, DJ, and
  // Abecasis, GR (2005) A Note on Exact Tests of Hardy-Weinberg
  // Equilibrium. American Journal of Human Genetics. 76: 887 - 893.
  //
  // The original version was written by Jan Wigginton.
  //
  // This version was written by Christopher Chang.  It contains the following
  // improvements over the original SNPHWE():
  // 1. Proper handling of >64k genotypes.  Previously, there was a potential
  //    integer overflow.
  // 2. Detection and efficient handling of floating point overflow and
  //    underflow.  E.g. instead of summing a tail all the way down, the loop
  //    stops once the latest increment underflows the partial sum's 53-bit
  //    precision; this results in a large speedup when max heterozygote count
  //    >1k.
  // 3. No malloc() call in most cases: it's only necessary to keep track of a
  //    few partial sums.  (But see (6).)
  // 4. Support for the mid-p variant of this test.  See Graffelman J, Moreno V
  //    (2013) The mid p-value in exact tests for Hardy-Weinberg equilibrium.
  // 5. Log-p-value return (added Jan 2024).  A p-value of 1e-400 may be worth
  //    distinguishing from 1e-40000 in a biobank-scale dataset.
  // 6. Accurate handling of near-ties (added Mar 2026), using the QD and GMP
  //    libraries.  The QD library is used to efficiently perform computations
  //    on log-factorials with >60 bits of accuracy past the decimal point;
  //    this lets us jump from one tail to the other with negligible accuracy
  //    loss, and correctly resolve almost all near-ties without a full-blown
  //    rational expansion.  The GMP library supports the rational expansion
  //    when we need it; this can theoretically require large memory
  //    allocations, but that is very unlikely to come up in practice.
  //    (Note that we continue to allow the returned value to have a few bits
  //    of floating-point error.  But misclassification of a near-tie can
  //    result in a *large* relative error, so I've decided to go through the
  //    trouble of stamping that out despite its low analytical impact.)
  //
  // Note that the HweThreshLn() function is a lot more efficient for testing
  // against a p-value inclusion threshold.  HweLnP() should only be used if
  // you need the actual p-value.

  // Variables are mostly a mix of int32_ts and doubles, with a few
  // log-factorial-sum dd_reals (i.e. double-doubles).
  // int32 -> double casts are usually left implicit.
  // Naming conventions:
  // * obs_...: int32
  // * ..._ct: int32
  // * ..._ctd: double
  // * ..._ddr: dd_real
  // * Almost everything else is a double.
  int32_t obs_homc;
  int32_t obs_homr;
  if (obs_hom1 < obs_hom2) {
    obs_homc = obs_hom2;
    obs_homr = obs_hom1;
  } else {
    obs_homc = obs_hom1;
    obs_homr = obs_hom2;
  }
  const int32_t rare_ct = 2 * obs_homr + obs_hets;
  if (rare_ct < 2) {
    *resultp = midp * (-kLn2);
    return 0;
  }
  // MAF: rare_ct / allele_ct
  // modal #hets:
  //   sample_ct * 2 * MAF * (1 - MAF)
  // = rare_ct * (1 - MAF)
  const int32_t sample_ct = obs_hom1 + obs_hom2 + obs_hets;
  // allele_ct can be >= 2^31.
  const double allele_ctd = S_CAST(double, sample_ct * 2LL);
  const double maf = rare_ct / allele_ctd;
  // possible todo: check whether this type of expression is worth rewriting as
  // e.g. prefer_fma(rare_ct, -maf, rare_ct)
  // 'c' in cmodal_nhet is for 'continuous', i.e. the mode if we extend the
  // likelihood function to all reals in [0, sample_ct].
  const double cmodal_nhet = rare_ct * (1 - maf);
  double hets = obs_hets;
  double homr = obs_homr;
  double homc = obs_homc;
  double lastp = 1;
  double tailp = 1 - midp * 0.5;
  if (hets > cmodal_nhet) {
    const double het_delta = hets - cmodal_nhet;
    // From Feb 2024 - Feb 2026, (except for the p=1 fast path) we always
    // computed starting log-likelihood.  Then, if it was high enough, we
    // proceeded with computing 1 - [sum of center likelihoods]; otherwise we
    // summed both tails in a manner that could handle likelihoods < DBL_MIN.
    //
    // Numerical stability of both branches was investigated in Mar 2026.
    // * The slightly-unstable 1 - [sum of center likelihoods] branch was
    //   reverted to the old cancellation-avoiding relative-likelihood
    //   algorithm, and the entrance condition changed to het_delta < 344.
    //   This doesn't risk center-likelihood overflow, since the
    //   relative-likelihood of the mode is loosely bounded above by
    //     ((172^172) / 172!)^4 ~= 5.3e+292
    //   which leaves enough headroom to accumulate the rest of the center-sum
    //   and multiply by e.g. allele_ct without overflowing.
    // * Log-factorial computations in the tail-jumping branch are now
    //   performed with "double-double" precision.  (Possible todo: benchmark
    //   float128 on x86_64.  But web search results imply QD is better.)
    if ((!midp) && (het_delta < 2.0)) {
      // Fast path for p=1.
      if (obs_hets * (obs_hets - 1LL) <= 4 * (obs_homc + 1LL) * (obs_homr + 1)) {
        *resultp = 0;
        return 0;
      }
    }
    // Iterate outward to floating-point precision limit.
    // No need for homr > 0 check, tailp == preaddp will trigger when we hit 0.
    // (hets, homr, and homc never accumulate any floating-point error since
    // they start as small integers and are only changed by adding/subtracting
    // 1 or 2.)
    while (1) {
      hets += 2;
      lastp *= (4 * homr * homc) / (hets * (hets - 1));
      homr -= 1;
      homc -= 1;
      const double preaddp = tailp;
      tailp += lastp;
      if (tailp == preaddp) {
        break;
      }
    }
    if (het_delta < 344.0) {
      // Jump back to starting contingency table, and iterate inward.
      lastp = 1;
      hets = obs_hets;
      homr = obs_homr;
      homc = obs_homc;
      dd_real starting_lnprobv_ddr = {{DBL_MAX, 0.0}};
      double centerp = midp * 0.5;
      // No need for hets > 1 check, lastp checks do what we need.
      while (1) {
        homr += 1;
        homc += 1;
        lastp *= (hets * (hets - 1)) / (4 * homr * homc);
        hets -= 2;
        // Number of center tables is maximized with obs_hets - cmodal_nhet ~=
        // 344, obs_homr = 0, obs_homc and obs_hets both large.
        // Since 1 + 1/2 + ... + 1/172 < 1/173 + ... + 1/53000, we're limited
        // to ~53000 tables.  Each lastp update involves 4 operations which can
        // each introduce up to 0.5 ULP relative error under the default
        // rounding mode.  One ULP <= 2^{-52};
        //   (1 + 2^{-52})^k < 1 + (k+1)2^{-52}
        // for k < tens of millions, so in most cases we can safely add
        // error-bounds together as long as we add an extra 2^{-52} at the
        // beginning.
        if (lastp < 1 + 53000 * 2 * k2m52) {
          if (lastp <= 1 - 53000 * 2 * k2m52) {
            tailp += lastp;
            break;
          }
          // Near-tie.  True value of lastp can be greater than, equal to, or
          // less than 1.
          const int32_t hom_decr = obs_homr - S_CAST(int32_t, homr);
          intptr_t cmp_result;
          if (unlikely(HweCompare(obs_hets, obs_homr, obs_homc, hom_decr, &starting_lnprobv_ddr, &cmp_result, &lastp))) {
            return 1;
          }
          if (cmp_result <= 0) {
            tailp += lastp;
            if (midp && (cmp_result == 0)) {
              tailp -= 0.5;
              centerp += 0.5;
            }
            break;
          }
        }
        centerp += lastp;
      }
      // Continue down tail to floating-point precision limit.
      // No need for hets > 1 check, tailp == preaddp check does what we need
      // (even when hets is already -1 when entering the loop: in that case
      // lastp is 0).
      while (1) {
        homr += 1;
        homc += 1;
        lastp *= (hets * (hets - 1)) / (4 * homr * homc);
        hets -= 2;
        const double preaddp = tailp;
        tailp += lastp;
        if (tailp == preaddp) {
          break;
        }
      }
      *resultp = log(tailp / (tailp + centerp));
      return 0;
    }
    // starting_lnprobv_ddr is guaranteed to be negative for hets >= 4, and no
    // larger than ln(2) otherwise.
    const double c_minus_r = homc - homr;
    dd_real starting_lnprobv_ddr =
      ddr_sub(ddr_muld(_ddr_log2, obs_hets),
              ddr_add3_lfacts(obs_hets, obs_homr, obs_homc));
    // Now we want to jump near the other tail, without evaluating that many
    // contingency table log-likelihoods along the way.
    //
    // Each full log-likelihood evaluation requires 3 ddr_lfact() calls.  Since
    // they are now performed with extra precision, they require hundreds of
    // floating-point operations, so we want to limit ourselves to 1-2 full
    // evaluations most of the time.  (Possible todo: use lower-accuracy
    // Lfact() to jump around, followed by ddr_lfact() when exiting the loop.
    // Should be an easy performance win, but there's a complexity cost so I'll
    // wait until I see a scenario where this branch executes frequently...)
    //
    // The current heuristic starts by reflecting (obs_homr + homr) * 0.5
    // across the (continuous) mode, performing a full log-likelihood check at
    // an adjacent valid point.  Hopefully we find that we're in
    // (starting_lnprob - 62 * kLn2, starting_lnprob], so we're at or near a
    // table that actually contributes to the tail-sum.  (This window is chosen
    // to be wide enough to guarantee that at least one point falls inside when
    // sample_ct < 2^31.)
    //
    // If not, we jump again, using Newton's method.
    // If homr is too low (i.e. current log-likelihood is too high), when we
    // increase homr by 1, the likelihood gets multiplied by
    //   hets * (hets-1) / (4 * (homr+1) * (homc+1))
    // i.e. we're adding the logarithm of this value to the log-likelihood.
    // If homr is too high, when we decrease homr by 1, the likelihood gets
    // multiplied by
    //   4 * homr * homc / ((hets+2) * (hets+1))
    // We use the log of the first expression as the Newton's method f'(x) when
    // we're jumping to higher homr, and the negative-log of the second
    // expression when we're jumping to lower homr.
    // f''(x) is always negative, so we can aim for starting_lnprob instead of
    // the middle of the interval.

    // hets moves twice as fast as homr.  So if we add
    //   0.5 * (hets + obs_hets) - cmodal_nhet
    // to 0.5 * (homr + obs_homr), that reflects homr across the cmode.
    const double max_homr = S_CAST(double, rare_ct >> 1);
    {
      const double delta = 0.5 * (hets + obs_hets) - cmodal_nhet;
      homr = 0.5 * (homr + obs_homr) + delta;
      // Round up (to guarantee we've actually moved to the other side of the
      // cmode) and clamp.
      homr = ceil_smalleps_limit32(homr, max_homr);
    }
    // 'f' in lnprobf refers to fixed component of log-likelihood (unchanged
    // when rare_ct and sample_ct are held constant), 'v' refers to variable
    // component
    const dd_real lnprobf_ddr =
      ddr_sub(ddr_add3_lfacts(sample_ct, rare_ct, allele_ctd - rare_ct),
              ddr_lfact(allele_ctd));
    const double starting_lnprob = ddr_add(lnprobf_ddr, starting_lnprobv_ddr).x[0];
    while (1) {
      hets = rare_ct - homr * 2;
      homc = homr + c_minus_r;
      const dd_real lnprobv_ddr =
        ddr_sub(ddr_muld(_ddr_log2, hets),
                ddr_add3_lfacts(hets, homr, homc));
      const double lnprob_diff = ddr_sub(lnprobv_ddr, starting_lnprobv_ddr).x[0];
      // Could tighten this threshold further; I haven't performed a careful
      // error analysis yet but CompareFactorialProducts() includes a plausible
      // assumption that 2^{-60} is safe.  But code is correct as long as we're
      // guaranteed to enter the "lastp < 2 - one_minus_scaled_eps" branch for
      // positive lnprob_diff.
      if (lnprob_diff >= k2m53) {
        if (homr == max_homr) {
          // All tables on this tail are larger than the starting table.  Exit.
          // (This is possible when obs_hom1 == obs_hom2 == 0.)
          *resultp = starting_lnprob + log(tailp);
          return 0;
        }
        const double ll_deriv = log(hets * (hets - 1) / (4 * (homr + 1) * (homc + 1)));
        // Round absolute value up, to guarantee that we make progress.
        // (lnprob_diff is positive and ll_deriv is negative.)
        // This may overshoot.  But the function is guaranteed to terminate
        // because we never overshoot (and we do always make progress on each
        // step) once we're on the other side.
        homr += ceil_smalleps(-lnprob_diff / ll_deriv);
        if (homr > max_homr) {
          homr = max_homr;
        }
      } else if (lnprob_diff > -62 * kLn2) {
        lastp = exp(lnprob_diff);
        break;
      } else {
        const double ll_deriv = log((hets + 2) * (hets + 1) / (4 * homr * homc));
        // Round down, to guarantee we don't overshoot.
        // We're guaranteed to make progress, since lnprob_diff <= -62 * log(2)
        // and sample_ct < 2^31.
        homr -= S_CAST(int64_t, lnprob_diff / ll_deriv);
      }
    }
    // Sum toward center, until lastp >= 1.  (No more risk of double-counting
    // the starting table, since we don't enter this branch at all unless the
    // starting table is >= 172 steps from the cmode.)
    //
    // lastp should be accurate to 3 ULP as we enter this loop (max 1.5 ULP
    // observed error from exp, tiny bit over 0.5 from lnprob_diff, we round up
    // all the way to 3 so we don't have to worry about "2 -
    // one_minus_scaled_eps" rounding behavior), so near-tie detection can use
    // a tight epsilon here.
    double one_minus_scaled_eps = 1 - 3 * k2m52;
    const double lastp_tail = lastp;
    const double homr_tail = homr;
    const double homc_tail = homc;
    const double hets_tail = hets;
    while (lastp <= one_minus_scaled_eps) {
      tailp += lastp;
      hets += 2;
      lastp *= (4 * homr * homc) / (hets * (hets - 1));
      homr -= 1;
      homc -= 1;
      one_minus_scaled_eps -= 2 * k2m52;
    }
    if (lastp < 2 - one_minus_scaled_eps) {
      const int32_t hom_decr = obs_homr - S_CAST(int32_t, homr);
      intptr_t cmp_result;
      if (unlikely(HweCompare(obs_hets, obs_homr, obs_homc, hom_decr, &starting_lnprobv_ddr, &cmp_result, &lastp))) {
        return 1;
      }
      if (cmp_result <= 0) {
        tailp += lastp;
        if (midp && (cmp_result == 0)) {
          tailp -= 0.5;
        }
      }
    }
    // Sum away from center, until sums stop changing.
    lastp = lastp_tail;
    homr = homr_tail;
    homc = homc_tail;
    hets = hets_tail;
    while (1) {
      homr += 1;
      homc += 1;
      lastp *= (hets * (hets - 1)) / (4 * homr * homc);
      const double preaddp = tailp;
      tailp += lastp;
      if (tailp == preaddp) {
        break;
      }
      hets -= 2;
    }
    *resultp = starting_lnprob + log(tailp);
    return 0;
  }
  // Same as above, just with directions flipped.
  const double het_delta = cmodal_nhet - hets;
  if ((!midp) && (het_delta < 2.0)) {
    // Fast path for p=1.
    if ((4LL * obs_homr) * obs_homc <= (obs_hets + 2LL) * (obs_hets + 1LL)) {
      *resultp = 0;
      return 0;
    }
  }
  // Iterate outward to floating-point precision limit.
  while (1) {
    homr += 1;
    homc += 1;
    lastp *= (hets * (hets - 1)) / (4 * homr * homc);
    hets -= 2;
    const double preaddp = tailp;
    tailp += lastp;
    if (tailp == preaddp) {
      break;
    }
  }
  if (het_delta < 344.0) {
    // Jump back to starting table, and iterate inward.
    lastp = 1;
    hets = obs_hets;
    homr = obs_homr;
    homc = obs_homc;
    dd_real starting_lnprobv_ddr = {{DBL_MAX, 0.0}};
    double centerp = midp * 0.5;
    // No need for hets > 1 check, lastp checks do what we need.
    while (1) {
      hets += 2;
      lastp *= (4 * homr * homc) / (hets * (hets - 1));
      homr -= 1;
      homc -= 1;
      // If we're 172 steps from the center, number of center tables is limited
      // to ~2*172 = 344, when obs_homr ~= obs_homc.
      if (lastp < 1 + (344 * 2 + 1) * k2m52) {
        if (lastp <= 1 - (344 * 2 + 1) * k2m52) {
          tailp += lastp;
          break;
        }
        // Near-tie.  True value of lastp can be greater than, equal to, or
        // less than 1.
        const int32_t hom_decr = obs_homr - S_CAST(int32_t, homr);
        intptr_t cmp_result;
        if (unlikely(HweCompare(obs_hets, obs_homr, obs_homc, hom_decr, &starting_lnprobv_ddr, &cmp_result, &lastp))) {
          return 1;
        }
        if (cmp_result <= 0) {
          tailp += lastp;
          if (midp && (cmp_result == 0)) {
            tailp -= 0.5;
            centerp += 0.5;
          }
          break;
        }
      }
      centerp += lastp;
    }
    // Continue down tail to floating-point precision limit.
    while (1) {
      hets += 2;
      lastp *= (4 * homr * homc) / (hets * (hets - 1));
      homr -= 1;
      homc -= 1;
      const double preaddp = tailp;
      tailp += lastp;
      if (tailp == preaddp) {
        break;
      }
    }
    *resultp = log(tailp / (tailp + centerp));
    return 0;
  }
  const double c_minus_r = homc - homr;
  dd_real starting_lnprobv_ddr =
    ddr_sub(ddr_muld(_ddr_log2, obs_hets),
            ddr_add3_lfacts(obs_hets, obs_homr, obs_homc));
  // Jump to other tail.
  {
    const double delta = cmodal_nhet - 0.5 * (hets + obs_hets);
    homr = 0.5 * (homr + obs_homr) - delta;
    // Round down (to guarantee we've actually moved to the other side of the
    // cmode) and clamp.
    homr = S_CAST(int32_t, homr);
    if (homr < 0) {
      homr = 0;
    }
  }
#ifndef NDEBUG
  const double max_homr = S_CAST(double, rare_ct >> 1);
#endif
  const dd_real lnprobf_ddr =
    ddr_sub(ddr_add3_lfacts(sample_ct, rare_ct, allele_ctd - rare_ct),
            ddr_lfact(allele_ctd));
  const double starting_lnprob = ddr_add(lnprobf_ddr, starting_lnprobv_ddr).x[0];
  while (1) {
    hets = rare_ct - homr * 2;
    homc = homr + c_minus_r;
    const dd_real lnprobv_ddr =
      ddr_sub(ddr_muld(_ddr_log2, hets),
              ddr_add3_lfacts(hets, homr, homc));
    const double lnprob_diff = ddr_sub(lnprobv_ddr, starting_lnprobv_ddr).x[0];
    if (lnprob_diff >= k2m53) {
      if (homr == 0) {
        // All tables on this tail have higher likelihood than the starting
        // table.  Exit.
        *resultp = starting_lnprob + log(tailp);
        return 0;
      }
      const double ll_deriv = log(4 * homr * homc / ((hets + 2) * (hets + 1)));
      homr -= ceil_smalleps(-lnprob_diff / ll_deriv);
      if (homr < 0) {
        homr = 0;
      }
    } else if (lnprob_diff > -62 * kLn2) {
      lastp = exp(lnprob_diff);
      break;
    } else {
      const double ll_deriv = log(4 * (homr + 1) * (homc + 1) / (hets * (hets - 1)));
      homr += S_CAST(int64_t, lnprob_diff / ll_deriv);
      assert(homr <= max_homr);
    }
  }
  // Sum toward center, until lastp >= 1.
  double one_minus_scaled_eps = 1 - 3 * k2m52;
  const double lastp_tail = lastp;
  const double homr_tail = homr;
  const double homc_tail = homc;
  const double hets_tail = hets;
  while (lastp <= one_minus_scaled_eps) {
    tailp += lastp;
    homr += 1;
    homc += 1;
    lastp *= hets * (hets - 1) / (4 * homr * homc);
    hets -= 2;
    one_minus_scaled_eps -= 2 * k2m52;
  }
  if (lastp < 2 - one_minus_scaled_eps) {
    const int32_t hom_decr = obs_homr - S_CAST(int32_t, homr);
    intptr_t cmp_result;
    if (unlikely(HweCompare(obs_hets, obs_homr, obs_homc, hom_decr, &starting_lnprobv_ddr, &cmp_result, &lastp))) {
      return 1;
    }
    if (cmp_result <= 0) {
      tailp += lastp;
      if (midp && (cmp_result == 0)) {
        tailp -= 0.5;
      }
    }
  }
  // Sum away from center, until sums stop changing.
  lastp = lastp_tail;
  homr = homr_tail;
  homc = homc_tail;
  hets = hets_tail;
  while (1) {
    hets += 2;
    lastp *= 4 * homr * homc / (hets * (hets - 1));
    const double preaddp = tailp;
    tailp += lastp;
    if (tailp == preaddp) {
      break;
    }
    homr -= 1;
    homc -= 1;
  }
  *resultp = starting_lnprob + log(tailp);
  return 0;
}

// 2^{-83} bias to give plink 1.9-style exact tests maximum ability to
// determine tiny p-values.  (~2^{-53} is necessary to take advantage of
// denormalized small numbers, then allow tail sum to be up to 2^30.  ...okay,
// HweThresh[Midp]() is not responsible for denormal values of thresh, and
// plink2 now just flushes denormals to zero.  But configuring this
// constant to be compatible with them doesn't cost us anything.)
static const double kExactTestBias = k2m50 / (1LL << 33);

BoolErr HweThresh(int32_t obs_hets, int32_t obs_hom1, int32_t obs_hom2, double pval_thresh, uint32_t* out_of_eqp) {
  // Threshold-test-only version of HweLnP() which is usually able to exit
  // from the calculation earlier.  Assumes DBL_MIN <= pval_thresh <= 1 (note
  // that some older versions of this function didn't handle pval_thresh=1
  // correctly, and plink2 still avoids calling this with pval_thresh=1).
  // Returns *out_of_eqp=0 if these counts are close enough to Hardy-Weinberg
  // equilibrium, 1 otherwise.
  //
  // BoolErr corresponds to malloc failure when resolving a near-tie.
  //
  // Suppose, for definiteness, that the number of observed hets is no less
  // than expectation.  (Same ideas apply for the other case.)  We proceed as
  // follows:
  // - Sum the *relative* likelihoods of more likely smaller het counts.
  // - Determine the minimum tail mass to pass the threshold.
  // - The majority of the time, the tail boundary elements are enough to pass
  //   the threshold; we never need to sum the remainder of the tails.
  // - And in the case of disequilibrium, we will often be able to immediately
  //   determine that the tail sum cannot possibly pass the threshold, just by
  //   looking at the tail boundary elements and using a geometric series to
  //   upper-bound the tail sums.
  // - Only when neither of these conditions hold do we start traveling down
  //   the tails.
  int32_t obs_homc;
  int32_t obs_homr;
  if (obs_hom1 < obs_hom2) {
    obs_homc = obs_hom2;
    obs_homr = obs_hom1;
  } else {
    obs_homc = obs_hom1;
    obs_homr = obs_hom2;
  }
  const int64_t genotypes2 = (obs_hets + obs_homc + obs_homr) * 2LL;
  if (!genotypes2) {
    *out_of_eqp = 0;
    return 0;
  }
  const int32_t rare_copies = 2 * obs_homr + obs_hets;
  double hets_t2 = obs_hets;  // tail 2
  double homr_t2 = obs_homr;
  double homc_t2 = obs_homc;

  double tailp1 = kExactTestBias;
  double centerp = 0;
  double lastp2 = tailp1;
  double tailp2 = 0;

  // const double center_div_tail_thresh = (1 - pval_thresh) / pval_thresh;
  // An initial upper bound on the tail sum is useful, since it lets us
  // report test failure before summing the entire center.  We use the
  // trivial bound of 1 + floor(rare_copies / 2): that's the total number
  // of possible het counts, and the relative probability for each count must
  // be <= 1 if it's in the tail.
  const double center_div_tail_thresh = (1 - pval_thresh) / pval_thresh;
  const double centerp_exit_thresh = u31tod(1 + (rare_copies >> 1)) * (center_div_tail_thresh * kExactTestBias);
  double scaled_one_plus_eps = kExactTestBias * (1 + k2m52);

  // Expected het count:
  //   2 * rarefreq * (1 - rarefreq) * genotypes
  // = 2 * (rare_copies / (2 * genotypes)) * (1 - rarefreq) * genotypes
  // = rare_copies * (1 - (rare_copies / (2 * genotypes)))
  // = (rare_copies * (2 * genotypes - rare_copies)) / (2 * genotypes)
  //
  // The computational identity is
  //   P(nhets == n) := P(nhets == n+2) * (n+2) * (n+1) /
  //                    (4 * homr(n) * homc(n))
  // where homr() and homc() are the number of homozygous rares/commons needed
  // to maintain the same allele frequencies.
  // This probability is always decreasing when proceeding away from the
  // expected het count.
  if (obs_hets * genotypes2 > rare_copies * (genotypes2 - rare_copies)) {
    // tail 1 = upper hets
    if (obs_hets < 2) {
      *out_of_eqp = 0;
      return 0;
    }

    // het_probs[hets] = 1
    // het_probs[hets - 2] = het_probs[hets] * hets * (hets - 1) / (4 * (homr + 1) * (homc + 1))
    do {
      homr_t2 += 1;
      homc_t2 += 1;
      lastp2 *= (hets_t2 * (hets_t2 - 1)) / (4 * homr_t2 * homc_t2);
      hets_t2 -= 2;
      scaled_one_plus_eps += kExactTestBias * 2 * k2m52;
      if (lastp2 < scaled_one_plus_eps) {
        if (lastp2 <= 2 * kExactTestBias - scaled_one_plus_eps) {
          tailp2 = lastp2;
          break;
        }
        double lastp = lastp2 * (1.0 / kExactTestBias);
        const int32_t hom_decr = obs_homr - S_CAST(int32_t, homr_t2);
        dd_real starting_lnprobv_ddr = {{DBL_MAX, 0.0}};
        intptr_t cmp_result;
        if (unlikely(HweCompare(obs_hets, obs_homr, obs_homc, hom_decr, &starting_lnprobv_ddr, &cmp_result, &lastp))) {
          return 1;
        }
        lastp2 = lastp * kExactTestBias;
        if (cmp_result <= 0) {
          tailp2 = lastp2;
          break;
        }
        // lastp can be of the form exp(lnprob_diff).  We are conservatively
        // assuming that exp() introduces up to 1.5 ULP error, and lnprob_diff
        // may be off by slightly over 0.5 ULP.
        scaled_one_plus_eps = kExactTestBias * (1 + 3 * k2m52);
      }
      centerp += lastp2;
      if (centerp >= centerp_exit_thresh) {
        *out_of_eqp = 1;
        return 0;
      }
    } while (hets_t2 > 1);
    // hets_t2 guaranteed to be nonnegative on loop exit, so ratio
    // calculation works

    // This is NaN when pval_thresh=1, so we write the next if-condition to be
    // true on NaN.
    const double tailp_exit_thresh = centerp / center_div_tail_thresh;
    if (!(tailp1 + tailp2 < tailp_exit_thresh)) {
      *out_of_eqp = 0;
      return 0;
    }
    // c + cr + cr^2 + ... = c/(1-r), which is an upper bound for the tail sum
    const double ratio = (hets_t2 * (hets_t2 - 1)) / (4 * (homr_t2 + 1) * (homc_t2 + 1));
    const double tail2_ceil = tailp2 / (1 - ratio);
    double hets_t1 = obs_hets + 2;
    double homr_t1 = obs_homr;
    double homc_t1 = obs_homc;
    // ratio for the other tail
    double lastp1 = (4 * homr_t1 * homc_t1) / (hets_t1 * (hets_t1 - 1));
    const double tail1_ceil = tailp1 / (1 - lastp1);
    if (tail1_ceil + tail2_ceil < tailp_exit_thresh) {
      *out_of_eqp = 1;
      return 0;
    }
    lastp1 *= tailp1;
    tailp1 += lastp1;

    if (obs_homr > 1) {
      // het_probs[hets + 2] = het_probs[hets] * 4 * homr * homc / ((hets + 2) * (hets + 1))
      const double tailp1_exit_thresh = tailp_exit_thresh - tailp2;
      while (1) {
        hets_t1 += 2;
        homr_t1 -= 1;
        homc_t1 -= 1;
        lastp1 *= (4 * homr_t1 * homc_t1) / (hets_t1 * (hets_t1 - 1));
        const double preaddp = tailp1;
        tailp1 += lastp1;
        if (tailp1 >= tailp1_exit_thresh) {
          *out_of_eqp = 0;
          return 0;
        }
        // homr_t1 == 1 check isn't necessary, but it provides a
        // noticeable speedup in my testing on real data.
        if ((tailp1 == preaddp) || (homr_t1 == 1)) {
          break;
        }
      }
    }
    if (tailp1 + tail2_ceil < tailp_exit_thresh) {
      *out_of_eqp = 1;
      return 0;
    }
    const double tailp2_exit_thresh = tailp_exit_thresh - tailp1;
    while (1) {
      homr_t2 += 1;
      homc_t2 += 1;
      lastp2 *= (hets_t2 * (hets_t2 - 1)) / (4 * homr_t2 * homc_t2);
      const double preaddp = tailp2;
      tailp2 += lastp2;
      if (tailp2 >= tailp2_exit_thresh) {
        *out_of_eqp = 0;
        return 0;
      }
      if (tailp2 == preaddp) {
        *out_of_eqp = 1;
        return 0;
      }
      hets_t2 -= 2;
    }
  }
  // tail 1 = lower hets
  if (!obs_homr) {
    *out_of_eqp = 0;
    return 0;
  }
  do {
    hets_t2 += 2;
    lastp2 *= (4 * homr_t2 * homc_t2) / (hets_t2 * (hets_t2 - 1));
    homr_t2 -= 1;
    homc_t2 -= 1;
    scaled_one_plus_eps += kExactTestBias * 2 * k2m52;
    if (lastp2 < scaled_one_plus_eps) {
      if (lastp2 <= 2 * kExactTestBias - scaled_one_plus_eps) {
        tailp2 = lastp2;
        break;
      }
      double lastp = lastp2 * (1.0 / kExactTestBias);
      const int32_t hom_decr = obs_homr - S_CAST(int32_t, homr_t2);
      dd_real starting_lnprobv_ddr = {{DBL_MAX, 0.0}};
      intptr_t cmp_result;
      if (unlikely(HweCompare(obs_hets, obs_homr, obs_homc, hom_decr, &starting_lnprobv_ddr, &cmp_result, &lastp))) {
        return 1;
      }
      lastp2 = lastp * kExactTestBias;
      if (cmp_result <= 0) {
        tailp2 = lastp2;
        break;
      }
      scaled_one_plus_eps = kExactTestBias * (1 + 3 * k2m52);
    }
    centerp += lastp2;
    if (centerp >= centerp_exit_thresh) {
      *out_of_eqp = 1;
      return 0;
    }
  } while (homr_t2 > 0);
  // homr_t2 guaranteed to be nonnegative on loop exit, so ratio
  // calculation works

  const double tailp_exit_thresh = centerp / center_div_tail_thresh;
  if (!(tailp1 + tailp2 < tailp_exit_thresh)) {
    *out_of_eqp = 0;
    return 0;
  }
  const double ratio = (4 * homr_t2 * homc_t2) / ((hets_t2 + 2) * (hets_t2 + 1));
  const double tail2_ceil = tailp2 / (1 - ratio);
  double hets_t1 = obs_hets;
  double homr_t1 = obs_homr + 1;
  double homc_t1 = obs_homc + 1;
  double lastp1 = (hets_t1 * (hets_t1 - 1)) / (4 * homr_t1 * homc_t1);
  const double tail1_ceil = tailp1 / (1 - lastp1);
  if (tail1_ceil + tail2_ceil < tailp_exit_thresh) {
    *out_of_eqp = 1;
    return 0;
  }
  lastp1 *= tailp1;
  tailp1 += lastp1;

  if (obs_hets >= 4) {
    const double tailp1_exit_thresh = tailp_exit_thresh - tailp2;
    while (1) {
      hets_t1 -= 2;
      homr_t1 += 1;
      homc_t1 += 1;
      lastp1 *= (hets_t1 * (hets_t1 - 1)) / (4 * homr_t1 * homc_t1);
      const double preaddp = tailp1;
      tailp1 += lastp1;
      if (tailp1 >= tailp1_exit_thresh) {
        *out_of_eqp = 0;
        return 0;
      }
      if ((tailp1 == preaddp) || (hets_t1 < 4)) {
        break;
      }
    }
  }
  if (tailp1 + tail2_ceil < tailp_exit_thresh) {
    *out_of_eqp = 1;
    return 0;
  }
  const double tailp2_exit_thresh = tailp_exit_thresh - tailp1;
  while (1) {
    hets_t2 += 2;
    lastp2 *= (4 * homr_t2 * homc_t2) / (hets_t2 * (hets_t2 - 1));
    homr_t2 -= 1;
    homc_t2 -= 1;
    const double preaddp = tailp2;
    tailp2 += lastp2;
    if (tailp2 >= tailp2_exit_thresh) {
      *out_of_eqp = 0;
      return 0;
    }
    if (tailp2 == preaddp) {
      *out_of_eqp = 1;
      return 0;
    }
  }
}

BoolErr HweThreshMidp(int32_t obs_hets, int32_t obs_hom1, int32_t obs_hom2, double pval_thresh, uint32_t* out_of_eqp) {
  // Mid-p version of HweThresh().  (There are enough fiddly differences that I
  // think it's better for this to be a separate function.)  Assumes
  // DBL_MIN <= pval_thresh < 0.5.
  int32_t obs_homc;
  int32_t obs_homr;
  if (obs_hom1 < obs_hom2) {
    obs_homc = obs_hom2;
    obs_homr = obs_hom1;
  } else {
    obs_homc = obs_hom1;
    obs_homr = obs_hom2;
  }
  const int64_t genotypes2 = (obs_hets + obs_homc + obs_homr) * 2LL;
  if (!genotypes2) {
    *out_of_eqp = 0;
    return 0;
  }
  int32_t rare_copies = 2 * obs_homr + obs_hets;
  double hets_t2 = obs_hets;  // tail 2
  double homr_t2 = obs_homr;
  double homc_t2 = obs_homc;
  double tailp1 = kExactTestBias * 0.5;
  double centerp = tailp1;
  double lastp2 = kExactTestBias;
  double tailp2 = 0;
  const double center_div_tail_thresh = (1 - pval_thresh) / pval_thresh;
  const double centerp_exit_thresh = u31tod(1 + (rare_copies >> 1)) * (center_div_tail_thresh * kExactTestBias);
  double scaled_one_plus_eps = kExactTestBias * (1 + k2m52);
  if (obs_hets * genotypes2 > rare_copies * (genotypes2 - rare_copies)) {
    if (obs_hets < 2) {
      *out_of_eqp = 0;
      return 0;
    }
    do {
      homr_t2 += 1;
      homc_t2 += 1;
      lastp2 *= (hets_t2 * (hets_t2 - 1)) / (4 * homr_t2 * homc_t2);
      hets_t2 -= 2;
      scaled_one_plus_eps += kExactTestBias * 2 * k2m52;
      if (lastp2 < scaled_one_plus_eps) {
        if (lastp2 <= 2 * kExactTestBias - scaled_one_plus_eps) {
          tailp2 = lastp2;
          break;
        }
        double lastp = lastp2 * (1.0 / kExactTestBias);
        const int32_t hom_decr = obs_homr - S_CAST(int32_t, homr_t2);
        dd_real starting_lnprobv_ddr = {{DBL_MAX, 0.0}};
        intptr_t cmp_result;
        if (unlikely(HweCompare(obs_hets, obs_homr, obs_homc, hom_decr, &starting_lnprobv_ddr, &cmp_result, &lastp))) {
          return 1;
        }
        lastp2 = lastp * kExactTestBias;
        if (cmp_result <= 0) {
          if (cmp_result == 0) {
            tailp2 = tailp1;
            centerp += tailp1;
          } else {
            tailp2 = lastp2;
          }
          break;
        }
        scaled_one_plus_eps = kExactTestBias * (1 + 3 * k2m52);
      }
      centerp += lastp2;
      if (centerp >= centerp_exit_thresh) {
        *out_of_eqp = 1;
        return 0;
      }
    } while (hets_t2 > 1);

    const double tailp_exit_thresh = centerp / center_div_tail_thresh;
    if (!(tailp1 + tailp2 < tailp_exit_thresh)) {
      *out_of_eqp = 0;
      return 0;
    }
    const double ratio = (hets_t2 * (hets_t2 - 1)) / (4 * (homr_t2 + 1) * (homc_t2 + 1));
    // this needs to work in both the tie and no-tie cases
    const double tail2_ceil = prefer_fma(lastp2, ratio / (1 - ratio), tailp2);
    double hets_t1 = obs_hets + 2;
    double homr_t1 = obs_homr;
    double homc_t1 = obs_homc;
    double lastp1 = (4 * homr_t1 * homc_t1) / (hets_t1 * (hets_t1 - 1));
    // always a tie here
    const double tail1_ceil = (tailp1 * 2) / (1 - lastp1) - tailp1;
    if (tail1_ceil + tail2_ceil < tailp_exit_thresh) {
      *out_of_eqp = 1;
      return 0;
    }
    lastp1 *= tailp1 * 2;
    tailp1 += lastp1;

    if (obs_homr > 1) {
      const double tailp1_exit_thresh = tailp_exit_thresh - tailp2;
      while (1) {
        hets_t1 += 2;
        homr_t1 -= 1;
        homc_t1 -= 1;
        lastp1 *= (4 * homr_t1 * homc_t1) / (hets_t1 * (hets_t1 - 1));
        const double preaddp = tailp1;
        tailp1 += lastp1;
        if (tailp1 >= tailp1_exit_thresh) {
          *out_of_eqp = 0;
          return 0;
        }
        if ((tailp1 == preaddp) || (homr_t1 == 1)) {
          break;
        }
      }
    }
    if (tailp1 + tail2_ceil < tailp_exit_thresh) {
      *out_of_eqp = 1;
      return 0;
    }
    const double tailp2_exit_thresh = tailp_exit_thresh - tailp1;
    while (1) {
      homr_t2 += 1;
      homc_t2 += 1;
      lastp2 *= (hets_t2 * (hets_t2 - 1)) / (4 * homr_t2 * homc_t2);
      const double preaddp = tailp2;
      tailp2 += lastp2;
      if (tailp2 >= tailp2_exit_thresh) {
        *out_of_eqp = 0;
        return 0;
      }
      if (tailp2 == preaddp) {
        *out_of_eqp = 1;
        return 0;
      }
      hets_t2 -= 2;
    }
  }
  if (!obs_homr) {
    *out_of_eqp = 0;
    return 0;
  }
  do {
    hets_t2 += 2;
    lastp2 *= (4 * homr_t2 * homc_t2) / (hets_t2 * (hets_t2 - 1));
    homr_t2 -= 1;
    homc_t2 -= 1;
    scaled_one_plus_eps += kExactTestBias * 2 * k2m52;
    if (lastp2 < scaled_one_plus_eps) {
      if (lastp2 <= 2 - scaled_one_plus_eps) {
        tailp2 = lastp2;
        break;
      }
      double lastp = lastp2 * (1.0 / kExactTestBias);
      const int32_t hom_decr = obs_homr - S_CAST(int32_t, homr_t2);
      dd_real starting_lnprobv_ddr = {{DBL_MAX, 0.0}};
      intptr_t cmp_result;
      if (unlikely(HweCompare(obs_hets, obs_homr, obs_homc, hom_decr, &starting_lnprobv_ddr, &cmp_result, &lastp))) {
        return 1;
      }
      lastp2 = lastp * kExactTestBias;
      if (cmp_result <= 0) {
        if (cmp_result == 0) {
          tailp2 = tailp1;
          centerp += tailp1;
        } else {
          tailp2 = lastp2;
        }
        break;
      }
      scaled_one_plus_eps = kExactTestBias * (1 + 3 * k2m52);
    }
    centerp += lastp2;
    if (centerp >= centerp_exit_thresh) {
      *out_of_eqp = 1;
      return 0;
    }
  } while (homr_t2 > 0);

  const double tailp_exit_thresh = centerp / center_div_tail_thresh;
  if (!(tailp1 + tailp2 < tailp_exit_thresh)) {
    *out_of_eqp = 0;
    return 0;
  }
  const double ratio = (4 * homr_t2 * homc_t2) / ((hets_t2 + 2) * (hets_t2 + 1));
  const double tail2_ceil = prefer_fma(lastp2, ratio / (1 - ratio), tailp2);
  double hets_t1 = obs_hets;
  double homr_t1 = obs_homr + 1;
  double homc_t1 = obs_homc + 1;
  double lastp1 = (hets_t1 * (hets_t1 - 1)) / (4 * homr_t1 * homc_t1);
  const double tail1_ceil = (2 * tailp1) / (1 - lastp1) - tailp1;
  lastp1 *= 2 * tailp1;
  tailp1 += lastp1;

  if (tail1_ceil + tail2_ceil < tailp_exit_thresh) {
    *out_of_eqp = 1;
    return 0;
  }
  if (obs_hets >= 4) {
    const double tailp1_exit_thresh = tailp_exit_thresh - tailp2;
    while (1) {
      hets_t1 -= 2;
      homr_t1 += 1;
      homc_t1 += 1;
      lastp1 *= (hets_t1 * (hets_t1 - 1)) / (4 * homr_t1 * homc_t1);
      const double preaddp = tailp1;
      tailp1 += lastp1;
      if (tailp1 >= tailp1_exit_thresh) {
        *out_of_eqp = 0;
        return 0;
      }
      if ((tailp1 == preaddp) || (hets_t1 < 4)) {
        break;
      }
    }
  }
  if (tailp1 + tail2_ceil < tailp_exit_thresh) {
    *out_of_eqp = 1;
    return 0;
  }
  const double tailp2_exit_thresh = tailp_exit_thresh - tailp1;
  while (1) {
    hets_t2 += 2;
    lastp2 *= (4 * homr_t2 * homc_t2) / (hets_t2 * (hets_t2 - 1));
    homr_t2 -= 1;
    homc_t2 -= 1;
    const double preaddp = tailp2;
    tailp2 += lastp2;
    if (tailp2 >= tailp2_exit_thresh) {
      *out_of_eqp = 0;
      return 0;
    }
    if (tailp2 == preaddp) {
      *out_of_eqp = 1;
      return 0;
    }
  }
}

BoolErr HweThreshLnMain(int32_t obs_hets, int32_t obs_hom1, int32_t obs_hom2, int32_t midp, double ln_thresh, uint32_t* out_of_eqp) {
  assert(ln_thresh < -708);
  // Threshold-test-only version of HweLnP() which is usually able to exit
  // from the calculation earlier.  Returns *out_of_eqp=0 if these counts are
  // close enough to Hardy-Weinberg equilibrium, 1 otherwise.
  //
  // Assumes ln_thresh < -708, otherwise the |curr_hets - cmodal_nhet| < 344
  // early-exit doesn't work.  (possible todo: make the 344 constant a function
  // of ln_thresh, and set it closer to optimally.)
  //
  // Caller is responsible for including a tolerance in ln_thresh when
  // appropriate.
  int32_t obs_homc;
  int32_t obs_homr;
  if (obs_hom1 < obs_hom2) {
    obs_homc = obs_hom2;
    obs_homr = obs_hom1;
  } else {
    obs_homc = obs_hom1;
    obs_homr = obs_hom2;
  }
  const int32_t rare_ct = 2 * obs_homr + obs_hets;
  const int32_t sample_ct = obs_hom1 + obs_hom2 + obs_hets;
  const double allele_ctd = S_CAST(double, sample_ct * 2LL);
  const double maf = rare_ct / allele_ctd;
  const double cmodal_nhet = rare_ct * (1 - maf);
  double hets = obs_hets;
  // Change this to "rare_ct < 2" if ln_thresh restriction is being loosened
  // (to e.g. compare results against HweThresh()).
  if (fabs(hets - cmodal_nhet) < 344.0) {
    *out_of_eqp = 0;
    return 0;
  }

  // 1. Compute log-likelihood of starting table.  This may be high enough on
  //    its own to immediately return 0.
  //    If likelihood is lower than threshold / <total # of tables>, we can
  //    immediately return 1.
  // 2. Determine tailsum we must hit to return 0.
  // 3. The rest follows HweLnP(), except with an extra geometric-series-based
  //    early-exit attempt near the end.
  const dd_real lnprobf_ddr =
    ddr_sub(ddr_add3_lfacts(sample_ct, rare_ct, allele_ctd - rare_ct),
            ddr_lfact(allele_ctd));
  double homr = obs_homr;
  double homc = obs_homc;
  dd_real starting_lnprobv_ddr =
    ddr_sub(ddr_muld(_ddr_log2, hets),
            ddr_add3_lfacts(hets, homr, homc));
  const double starting_lnprob = ddr_add(lnprobf_ddr, starting_lnprobv_ddr).x[0];
  if (ln_thresh <= starting_lnprob - midp * kLn2) {
    *out_of_eqp = 0;
    return 0;
  }
  const double max_homr = S_CAST(double, rare_ct >> 1);
  if (ln_thresh > starting_lnprob + log(max_homr + 1 - midp * 0.5)) {
    *out_of_eqp = 1;
    return 0;
  }

  const double c_minus_r = homc - homr;
  // This should be in (0.5, 2^30].
  const double tail_thresh = exp(ln_thresh - starting_lnprob);
  double tailp = 1 - midp * 0.5;
  double lastp = 1;
  if (hets > cmodal_nhet) {
    // No center-sum (or p=1) code path, since it doesn't make sense to choose
    // a HWE threshold that makes these relevant; we should have already
    // returned 0.
    // (So we can't assume obs_hets >= 2.)
    while (1) {
      hets += 2;
      lastp *= (4 * homr * homc) / (hets * (hets - 1));
      const double preaddp = tailp;
      tailp += lastp;
      if (tailp == preaddp) {
        break;
      }
      if (tailp >= tail_thresh) {
        *out_of_eqp = 0;
        return 0;
      }
      homr -= 1;
      homc -= 1;
    }
    {
      const double delta = 0.5 * (hets + obs_hets) - cmodal_nhet;
      homr = 0.5 * (homr + obs_homr) + delta;
      homr = ceil_smalleps_limit32(homr, max_homr);
    }
    while (1) {
      hets = rare_ct - homr * 2;
      homc = homr + c_minus_r;
      const dd_real lnprobv_ddr =
        ddr_sub(ddr_muld(_ddr_log2, hets),
                ddr_add3_lfacts(hets, homr, homc));
      const double lnprob_diff = ddr_sub(lnprobv_ddr, starting_lnprobv_ddr).x[0];
      if (lnprob_diff >= k2m53) {
        if (homr == max_homr) {
          *out_of_eqp = 1;
          return 0;
        }
        const double ll_deriv = log(hets * (hets - 1) / (4 * (homr + 1) * (homc + 1)));
        homr += ceil_smalleps(-lnprob_diff / ll_deriv);
        if (homr > max_homr) {
          homr = max_homr;
        }
      } else if (lnprob_diff > -62 * kLn2) {
        lastp = exp(lnprob_diff);
        break;
      } else {
        const double ll_deriv = log((hets + 2) * (hets + 1) / (4 * homr * homc));
        homr -= S_CAST(int64_t, lnprob_diff / ll_deriv);
        assert(homr >= 0);
      }
    }
    double one_minus_scaled_eps = 1 - 3 * k2m52;
    const double lastp_tail = lastp;
    const double homr_tail = homr;
    const double homc_tail = homc;
    const double hets_tail = hets;
    while (lastp <= one_minus_scaled_eps) {
      tailp += lastp;
      if (tailp >= tail_thresh) {
        *out_of_eqp = 0;
        return 0;
      }
      hets += 2;
      lastp *= (4 * homr * homc) / (hets * (hets - 1));
      homr -= 1;
      homc -= 1;
      one_minus_scaled_eps -= 2 * k2m52;
    }
    if (lastp < 2 - one_minus_scaled_eps) {
      const int32_t hom_decr = obs_homr - S_CAST(int32_t, homr);
      intptr_t cmp_result;
      if (unlikely(HweCompare(obs_hets, obs_homr, obs_homc, hom_decr, &starting_lnprobv_ddr, &cmp_result, &lastp))) {
        return 1;
      }
      if (cmp_result <= 0) {
        if (cmp_result == 0) {
          tailp += 1 - midp * 0.5;
        } else {
          tailp += lastp;
        }
        if (tailp >= tail_thresh) {
          *out_of_eqp = 0;
          return 0;
        }
      }
    }
    // We're down to one tail that can be tightly bounded by a geometric
    // series.  (ratio is always decreasing)
    // c + cr + cr^2 + ... = c/(1-r)
    lastp = lastp_tail;
    homr = homr_tail;
    homc = homc_tail;
    hets = hets_tail;

    homr += 1;
    homc += 1;
    const double cur_ratio = (hets * (hets - 1)) / (4 * homr * homc);
    lastp *= cur_ratio;
    const double remaining_ceil = lastp_tail / (1 - cur_ratio);
    if (tailp + remaining_ceil < tail_thresh) {
      *out_of_eqp = 1;
      return 0;
    }
    while (1) {
      const double preaddp = tailp;
      tailp += lastp;
      if (tailp == preaddp) {
        *out_of_eqp = 1;
        return 0;
      }
      if (tailp >= tail_thresh) {
        *out_of_eqp = 0;
        return 0;
      }
      hets -= 2;
      homr += 1;
      homc += 1;
      lastp *= (hets * (hets - 1)) / (4 * homr * homc);
    }
    // unreachable
  }
  while (1) {
    homr += 1;
    homc += 1;
    lastp *= hets * (hets - 1) / (4 * homr * homc);
    hets -= 2;
    const double preaddp = tailp;
    tailp += lastp;
    if (tailp == preaddp) {
      break;
    }
    if (tailp >= tail_thresh) {
      *out_of_eqp = 0;
      return 0;
    }
  }
  {
    const double delta = cmodal_nhet - 0.5 * (hets + obs_hets);
    homr = 0.5 * (homr + obs_homr) - delta;
    homr = S_CAST(int32_t, homr);
    if (homr < 0) {
      homr = 0;
    }
  }
  while (1) {
    hets = rare_ct - homr * 2;
    homc = homr + c_minus_r;
    const dd_real lnprobv_ddr =
      ddr_sub(ddr_muld(_ddr_log2, hets),
              ddr_add3_lfacts(hets, homr, homc));
    const double lnprob_diff = ddr_sub(lnprobv_ddr, starting_lnprobv_ddr).x[0];
    if (lnprob_diff >= k2m53) {
      if (homr <= 0) {
        *out_of_eqp = 1;
        return 0;
      }
      const double ll_deriv = log(4 * homr * homc / ((hets + 2) * (hets + 1)));
      homr -= ceil_smalleps(-lnprob_diff / ll_deriv);
      if (homr < 0) {
        homr = 0;
      }
    } else if (lnprob_diff > -62 * kLn2) {
      lastp = exp(lnprob_diff);
      break;
    } else {
      const double ll_deriv = log(4 * (homr + 1) * (homc + 1) / (hets * (hets - 1)));
      homr += S_CAST(int64_t, lnprob_diff / ll_deriv);
      assert(homr <= max_homr);
    }
  }
  double one_minus_scaled_eps = 1 - 3 * k2m52;
  const double lastp_tail = lastp;
  const double homr_tail = homr;
  const double homc_tail = homc;
  const double hets_tail = hets;
  while (lastp <= one_minus_scaled_eps) {
    tailp += lastp;
    if (tailp >= tail_thresh) {
      *out_of_eqp = 0;
      return 0;
    }
    homr += 1;
    homc += 1;
    lastp *= hets * (hets - 1) / (4 * homr * homc);
    hets -= 2;
    one_minus_scaled_eps -= 2 * k2m52;
  }
  if (lastp < 2 - one_minus_scaled_eps) {
    const int32_t hom_decr = obs_homr - S_CAST(int32_t, homr);
    intptr_t cmp_result;
    if (unlikely(HweCompare(obs_hets, obs_homr, obs_homc, hom_decr, &starting_lnprobv_ddr, &cmp_result, &lastp))) {
      return 1;
    }
    if (cmp_result <= 0) {
      if (cmp_result == 0) {
        tailp += 1 - midp * 0.5;
      } else {
        tailp += lastp;
      }
      if (tailp >= tail_thresh) {
        *out_of_eqp = 0;
        return 0;
      }
    }
  }
  lastp = lastp_tail;
  homr = homr_tail;
  homc = homc_tail;
  hets = hets_tail;

  hets += 2;
  const double cur_ratio = 4 * homr * homc / (hets * (hets - 1));
  lastp *= cur_ratio;
  const double remaining_ceil = lastp / (1 - cur_ratio);
  if (tailp + remaining_ceil < tail_thresh) {
    *out_of_eqp = 1;
    return 0;
  }
  while (1) {
    const double preaddp = tailp;
    tailp += lastp;
    if (tailp == preaddp) {
      *out_of_eqp = 1;
      return 0;
    }
    if (tailp >= tail_thresh) {
      *out_of_eqp = 0;
      return 0;
    }
    homr -= 1;
    homc -= 1;
    hets += 2;
    lastp *= 4 * homr * homc / (hets * (hets - 1));
  }
}

// Switch between log- and regular representations at kSwitchThresh.
// 2^890 leaves enough headroom for at least 4 more multiplies by allele_ct.
// (Useful to reduce this to e.g. 2^150 when testing correctness, though it
// isn't currently safe to go all the way down to 1 due to algorithmic
// assumptions.)
static const double kSwitchThresh = k2p800 * k2p50 * (1LL << 40);
static const double kLnSwitchThresh = 890.0 * kLn2;
static const double kJumpThresh = 314.0; // chosen to guarantee base_prob < kSwitchThresh in non-jumping case
// static const double kSwitchThresh = k2p100 * k2p50;
// static const double kLnSwitchThresh = 150.0 * kLn2;

BoolErr HweLnFirstRow(int32_t obs_fhets, int32_t obs_fhom1, int32_t obs_fhom2, double* tailp_ptr, dd_real* starting_lnprobv_ddr_ptr, int32_t* tie_ct_ptr, double* orig_base_probl_ptr, double* orig_base_lnprobl_ptr, double* orig_base_epsl_ptr, double* orig_base_probr_ptr, double* orig_base_lnprobr_ptr, double* orig_base_epsr_ptr, double* orig_saved_lhets_ptr, double* orig_saved_lhom1_ptr, double* orig_saved_lhom2_ptr, double* orig_saved_rhets_ptr, double* orig_saved_rhom1_ptr, double* orig_saved_rhom2_ptr) {
  // possible todo: have this and HweLnP() call a shared function
  const int32_t female_ct = obs_fhets + obs_fhom1 + obs_fhom2;
  const double f1 = S_CAST(double, obs_fhets + 2LL * obs_fhom1);
  const double f2 = S_CAST(double, obs_fhets + 2LL * obs_fhom2);
  const double female_allele_ctd = S_CAST(double, female_ct * 2LL);
  *starting_lnprobv_ddr_ptr =
    ddr_sub(ddr_muld(_ddr_log2, obs_fhets),
            ddr_add3_lfacts(obs_fhets, obs_fhom1, obs_fhom2));
  const double cmodal_fhets = f1 * f2 / female_allele_ctd;
  double lastp = 1;
  int32_t tie_ct = 1;
  double hets = obs_fhets;
  double homr = obs_fhom1;
  double homc = obs_fhom2;
  double tailp = 1;
  double female_rare_ctd = f1;
  double c_minus_r = homc - homr;
  if (hets > cmodal_fhets) {
    *orig_base_probr_ptr = 1;
    *orig_base_epsr_ptr = k2m52;
    *orig_saved_rhets_ptr = obs_fhets;
    *orig_saved_rhom1_ptr = obs_fhom1;
    *orig_saved_rhom2_ptr = obs_fhom2;
    double* orig_saved_lhomr_ptr = orig_saved_lhom1_ptr;
    double* orig_saved_lhomc_ptr = orig_saved_lhom2_ptr;
    if (c_minus_r < 0.0) {
      orig_saved_lhomr_ptr = orig_saved_lhom2_ptr;
      orig_saved_lhomc_ptr = orig_saved_lhom1_ptr;
      female_rare_ctd = f2;
      c_minus_r = -c_minus_r;
      homr = obs_fhom2;
      homc = obs_fhom1;
    }
    const double obs_fhomr_d = homr;
    const double het_delta = hets - cmodal_fhets;
    while (1) {
      hets += 2;
      lastp *= (4 * homr * homc) / (hets * (hets - 1));
      homr -= 1;
      homc -= 1;
      const double preaddp = tailp;
      tailp += lastp;
      if (tailp == preaddp) {
        break;
      }
    }
    if (het_delta < kJumpThresh) {
      // Jump back to starting table, and iterate inward until we find the
      // start of the other tail.
      lastp = 1;
      hets = obs_fhets;
      homr = obs_fhomr_d;
      homc = obs_fhomr_d + c_minus_r;
      double one_plus_scaled_eps = 1 + k2m52;
      while (hets > 1) {
        homr += 1;
        homc += 1;
        lastp *= (hets * (hets - 1)) / (4 * homr * homc);
        hets -= 2;
        one_plus_scaled_eps += 2 * k2m52;
        if (lastp < one_plus_scaled_eps) {
          if (lastp <= 2 - one_plus_scaled_eps) {
            tailp += lastp;
            break;
          }
          // Near-tie.  True value of lastp can be greater than, equal to, or
          // less than 1.
          const int32_t hom_decr = S_CAST(int32_t, obs_fhomr_d - homr);
          intptr_t cmp_result;
          if (unlikely(HweCompare(obs_fhets, obs_fhom1, obs_fhom2, hom_decr, starting_lnprobv_ddr_ptr, &cmp_result, &lastp))) {
            return 1;
          }
          one_plus_scaled_eps = 1 + 3 * k2m52;
          if (cmp_result <= 0) {
            tailp += lastp;
            tie_ct += (cmp_result == 0);
            break;
          }
        }
      }
      *orig_saved_lhets_ptr = hets;
      *orig_saved_lhomr_ptr = homr;
      *orig_saved_lhomc_ptr = homc;
      *orig_base_probl_ptr = lastp;
      *orig_base_epsl_ptr = one_plus_scaled_eps - 1;
    } else {
      const double max_homr = S_CAST(int32_t, female_rare_ctd * 0.5);
      {
        const double delta = 0.5 * (hets + obs_fhets) - cmodal_fhets;
        // Round up and clamp.
        homr = 0.5 * (homr + obs_fhomr_d) + delta;
        homr = ceil_smalleps_limit32(homr, max_homr);
      }
      while (1) {
        hets = female_rare_ctd - homr * 2;
        homc = homr + c_minus_r;

        const dd_real lnprobv_ddr =
          ddr_sub(ddr_muld(_ddr_log2, hets),
                  ddr_add3_lfacts(hets, homr, homc));
        const dd_real lnprob_diff_ddr = ddr_sub(lnprobv_ddr, *starting_lnprobv_ddr_ptr);
        const double lnprob_diff = lnprob_diff_ddr.x[0];
        if (lnprob_diff >= k2m53) {
          if (homr == max_homr) {
            // All tables on this tail have higher likelihood than the starting
            // table.  Exit.
            assert(homr == max_homr);
            *tie_ct_ptr = 1;
            *orig_saved_lhets_ptr = hets;
            *orig_saved_lhomr_ptr = homr;
            *orig_saved_lhomc_ptr = homc;
            *tailp_ptr = tailp;
            if (lnprob_diff < kLnSwitchThresh) {
              *orig_base_probl_ptr = ddr_exp(lnprob_diff_ddr).x[0];
              *orig_base_epsl_ptr = 2 * k2m52;
            } else {
              *orig_base_probl_ptr = 0;
              *orig_base_lnprobl_ptr = lnprob_diff;
              *orig_base_epsl_ptr = (1 + ceil_smalleps(lnprob_diff)) * k2m52;
            }
            return 0;
          }
          const double ll_deriv = log(hets * (hets - 1) / (4 * (homr + 1) * (homc + 1)));
          // Round absolute value up, to guarantee that we make progress.
          // (lnprob_diff is positive and ll_deriv is negative.)
          // This may overshoot.  But the function is guaranteed to terminate
          // because we never overshoot (and we do always make progress on each
          // step) once we're on the other side.
          homr += ceil_smalleps(-lnprob_diff / ll_deriv);
          if (homr > max_homr) {
            homr = max_homr;
          }
        } else if (lnprob_diff > -62 * kLn2) {
          lastp = exp(lnprob_diff);
          break;
        } else {
          const double ll_deriv = log((hets + 2) * (hets + 1) / (4 * homr * homc));
          // Round down, to guarantee we don't overshoot.
          // We're guaranteed to make progress, since lnprob_diff <=
          // -62 * log(2) and female_ct < 2^31.
          homr -= S_CAST(int64_t, lnprob_diff / ll_deriv);
          assert(homr >= 0);
        }
      }
      // Sum toward center, until lastp >= 1.  (No more risk of double-counting
      // the starting table, since we separated out the started-near-mode
      // case.)
      double one_minus_scaled_eps = 1 - 3 * k2m52;
      const double lastp_tail = lastp;
      const double homr_tail = homr;
      const double homc_tail = homc;
      const double hets_tail = hets;
      while (lastp <= one_minus_scaled_eps) {
        tailp += lastp;
        hets += 2;
        lastp *= (4 * homr * homc) / (hets * (hets - 1));
        homr -= 1;
        homc -= 1;
        one_minus_scaled_eps -= 2 * k2m52;
      }
      if (lastp < 2 - one_minus_scaled_eps) {
        const int32_t hom_decr = S_CAST(int32_t, obs_fhomr_d - homr);
        intptr_t cmp_result;
        if (unlikely(HweCompare(obs_fhets, obs_fhom1, obs_fhom2, hom_decr, starting_lnprobv_ddr_ptr, &cmp_result, &lastp))) {
          return 1;
        }
        one_minus_scaled_eps = 1 - 3 * k2m52;
        if (cmp_result <= 0) {
          tailp += lastp;
          tie_ct += (cmp_result == 0);
        }
      }
      *orig_saved_lhets_ptr = hets;
      *orig_saved_lhomr_ptr = homr;
      *orig_saved_lhomc_ptr = homc;
      *orig_base_probl_ptr = lastp;
      *orig_base_epsl_ptr = 1 - one_minus_scaled_eps;
      lastp = lastp_tail;
      homr = homr_tail;
      homc = homc_tail;
      hets = hets_tail;
    }
    *tie_ct_ptr = tie_ct;
    // Sum away from center, until sums stop changing.
    while (1) {
      homr += 1;
      homc += 1;
      lastp *= (hets * (hets - 1)) / (4 * homr * homc);
      const double preaddp = tailp;
      tailp += lastp;
      if (tailp == preaddp) {
        break;
      }
      hets -= 2;
    }
    *tailp_ptr = tailp;
    return 0;
  }
  *orig_base_probl_ptr = 1;
  *orig_base_epsl_ptr = k2m52;
  *orig_saved_lhets_ptr = obs_fhets;
  *orig_saved_lhom1_ptr = obs_fhom1;
  *orig_saved_lhom2_ptr = obs_fhom2;
  double* orig_saved_rhomr_ptr = orig_saved_rhom1_ptr;
  double* orig_saved_rhomc_ptr = orig_saved_rhom2_ptr;
  if (c_minus_r < 0.0) {
    orig_saved_rhomr_ptr = orig_saved_rhom2_ptr;
    orig_saved_rhomc_ptr = orig_saved_rhom1_ptr;
    female_rare_ctd = f2;
    c_minus_r = -c_minus_r;
    homr = obs_fhom2;
    homc = obs_fhom1;
  }
  const double obs_fhomr_d = homr;
  const double het_delta = cmodal_fhets - hets;
  while (1) {
    homr += 1;
    homc += 1;
    lastp *= hets * (hets - 1) / (4 * homr * homc);
    hets -= 2;
    const double preaddp = tailp;
    tailp += lastp;
    if (tailp == preaddp) {
      break;
    }
  }
  if (het_delta < kJumpThresh) {
    lastp = 1;
    hets = obs_fhets;
    homr = obs_fhomr_d;
    homc = obs_fhomr_d + c_minus_r;
    double one_plus_scaled_eps = 1 + k2m52;
    while (homr > 0) {
      hets += 2;
      lastp *= (4 * homr * homc) / (hets * (hets - 1));
      homr -= 1;
      homc -= 1;
      one_plus_scaled_eps += 2 * k2m52;
      if (lastp < one_plus_scaled_eps) {
        if (lastp <= 2 - one_plus_scaled_eps) {
          tailp += lastp;
          break;
        }
        const int32_t hom_decr = S_CAST(int32_t, obs_fhomr_d - homr);
        intptr_t cmp_result;
        if (unlikely(HweCompare(obs_fhets, obs_fhom1, obs_fhom2, hom_decr, starting_lnprobv_ddr_ptr, &cmp_result, &lastp))) {
          return 1;
        }
        one_plus_scaled_eps = 1 + 3 * k2m52;
        if (cmp_result <= 0) {
          tailp += lastp;
          tie_ct += (cmp_result == 0);
          break;
        }
      }
    }
    *orig_saved_rhets_ptr = hets;
    *orig_saved_rhomr_ptr = homr;
    *orig_saved_rhomc_ptr = homc;
    *orig_base_probr_ptr = lastp;
    *orig_base_epsr_ptr = one_plus_scaled_eps - 1;
  } else {
    // Jump to other tail.
    {
      const double delta = cmodal_fhets - 0.5 * (hets + obs_fhets);
      homr = 0.5 * (homr + obs_fhomr_d) - delta;
      // Round down and clamp.
      homr = S_CAST(int32_t, homr);
      if (homr < 0) {
        homr = 0;
      }
    }
    while (1) {
      hets = female_rare_ctd - homr * 2;
      homc = homr + c_minus_r;
      const dd_real lnprobv_ddr =
        ddr_sub(ddr_muld(_ddr_log2, hets),
                ddr_add3_lfacts(hets, homr, homc));
      const dd_real lnprob_diff_ddr = ddr_sub(lnprobv_ddr, *starting_lnprobv_ddr_ptr);
      const double lnprob_diff = lnprob_diff_ddr.x[0];
      if (lnprob_diff >= k2m53) {
        if (homr <= 0) {
          // All tables on this tail have higher likelihood than the starting
          // table.  Exit.
          assert(homr == 0.0);
          *tie_ct_ptr = 1;
          *orig_saved_rhets_ptr = hets;
          *orig_saved_rhomr_ptr = homr;
          *orig_saved_rhomc_ptr = homc;
          *tailp_ptr = tailp;
          if (lnprob_diff < kLnSwitchThresh) {
            *orig_base_probr_ptr = ddr_exp(lnprob_diff_ddr).x[0];
            *orig_base_epsr_ptr = 2 * k2m52;
          } else {
            *orig_base_probr_ptr = 0;
            *orig_base_lnprobr_ptr = lnprob_diff;
            *orig_base_epsr_ptr = (1 + ceil_smalleps(lnprob_diff)) * k2m52;
          }
          return 0;
        }
        const double ll_deriv = log(4 * homr * homc / ((hets + 2) * (hets + 1)));
        homr -= ceil_smalleps(-lnprob_diff / ll_deriv);
        if (homr < 0) {
          homr = 0;
        }
      } else if (lnprob_diff > -62 * kLn2) {
        lastp = exp(lnprob_diff);
        break;
      } else {
        const double ll_deriv = log(4 * (homr + 1) * (homc + 1) / (hets * (hets - 1)));
        homr += S_CAST(int64_t, lnprob_diff / ll_deriv);
        assert(homr >= 0);
      }
    }
    // Sum toward center, until lastp >= 1.
    double one_minus_scaled_eps = 1 - 3 * k2m52;
    const double lastp_tail = lastp;
    const double homr_tail = homr;
    const double homc_tail = homc;
    const double hets_tail = hets;
    while (lastp <= one_minus_scaled_eps) {
      tailp += lastp;
      homr += 1;
      homc += 1;
      lastp *= hets * (hets - 1) / (4 * homr * homc);
      hets -= 2;
      one_minus_scaled_eps -= 2 * k2m52;
    }
    if (lastp < 2 - one_minus_scaled_eps) {
      const int32_t hom_decr = S_CAST(int32_t, obs_fhomr_d - homr);
      intptr_t cmp_result;
      if (unlikely(HweCompare(obs_fhets, obs_fhom1, obs_fhom2, hom_decr, starting_lnprobv_ddr_ptr, &cmp_result, &lastp))) {
        return 1;
      }
      one_minus_scaled_eps = 1 - 3 * k2m52;
      if (cmp_result <= 0) {
        tailp += lastp;
        tie_ct += (cmp_result == 0);
      }
    }
    *orig_saved_rhets_ptr = hets;
    *orig_saved_rhomr_ptr = homr;
    *orig_saved_rhomc_ptr = homc;
    *orig_base_probr_ptr = lastp;
    *orig_base_epsr_ptr = 1 - one_minus_scaled_eps;
    lastp = lastp_tail;
    homr = homr_tail;
    homc = homc_tail;
    hets = hets_tail;
  }
  *tie_ct_ptr = tie_ct;
  // Sum away from center, until sums stop changing.
  while (1) {
    hets += 2;
    lastp *= 4 * homr * homc / (hets * (hets - 1));
    const double preaddp = tailp;
    tailp += lastp;
    if (tailp == preaddp) {
      break;
    }
    homr -= 1;
    homc -= 1;
  }
  *tailp_ptr = tailp;
  return 0;
}

// - starting_lnprobv_ddr_ptr is expected to be initialized to
//     log(2^obs_fhets / (obs_fhets! obs_fhom1! obs_fhom2! obs_m1! obs_m2!))
//   or have x[0] initialized to DBL_MAX to indicate that the calculation
//   hasn't happened.  In the latter case, it may be set to the former value
//   if that is needed in the calculation.
//
// - *cmp_resultp is set to positive value if the (cur_fhom1, cur_fhom2)
//   contingency table has higher likelihood than (obs_fhom1, obs_fhom2), 0 if
//   identical likelihood, and negative value if lower likelihood.
//
// - Error is returned iff malloc fails.
BoolErr HweXchrCompare(uint32_t obs_fhets, uint32_t obs_fhom1, uint32_t obs_fhom2, uint32_t obs_m1, uint32_t obs_m2, uint32_t cur_fhom1, uint32_t cur_fhom2, int32_t m1_decr, dd_real* starting_lnprobv_ddr_ptr, intptr_t* cmp_resultp, double* dbl_ptr) {
  // From the Graffelman-Weir 2016 paper, likelihood is
  //
  //   (m1+f1)! (m2+f2)! (m1+m2)! (fhets+fhom1+fhom2)!    fhets
  //   ----------------------------------------------- * 2
  //     (m1+m2+f1+f2)! m1! m2! fhets! fhom1! fhom2!
  //
  // where f1 := fhets + 2*fhom1, and f2 := fhets + 2*fhom2.
  //
  // Likelihood ratio with all numerator terms (not counting 2^{fhets}) held
  // constant is then
  //
  //   obs_fhets! obs_fhom1! obs_fhom2! obs_m1! obs_m2!    cur_fhets-obs_fhets
  //   ------------------------------------------------ * 2
  //   cur_fhets! cur_fhom1! cur_fhom2! cur_m1! cur_m2!
  uint32_t numer_factorial_args[5];
  numer_factorial_args[0] = obs_fhets;
  numer_factorial_args[1] = obs_fhom1;
  numer_factorial_args[2] = obs_fhom2;
  numer_factorial_args[3] = obs_m1;
  numer_factorial_args[4] = obs_m2;
  const int32_t fhets_incr = obs_fhom1 + obs_fhom2 - cur_fhom1 - cur_fhom2;
  uint32_t denom_factorial_args[5];
  denom_factorial_args[0] = obs_fhets + fhets_incr;
  denom_factorial_args[1] = cur_fhom1;
  denom_factorial_args[2] = cur_fhom2;
  denom_factorial_args[3] = obs_m1 - m1_decr;
  denom_factorial_args[4] = obs_m2 + m1_decr;

  mp_limb_t* gmp_wkspace = nullptr;
  uintptr_t gmp_wkspace_limb_ct = 0;
  BoolErr reterr = CompareFactorialProducts(5, fhets_incr, obs_fhets, numer_factorial_args, denom_factorial_args, starting_lnprobv_ddr_ptr, &gmp_wkspace, &gmp_wkspace_limb_ct, cmp_resultp, dbl_ptr);
  free_cond(gmp_wkspace);

  return reterr;
}

// center_status set to 2 if all tables in current row have relative-likelihood
// <= 1, ties in current row have been counted, and we can prove all tables in
// subsequent rows in have relative-likelihood < 1.
// center_status set to 1 if only the first is true.
BoolErr HweXchrLnPLowHetTailsum(dd_real starting_lnprobv_ddr, uint32_t obs_fhets, uint32_t obs_fhom1, uint32_t obs_fhom2, uint32_t obs_m1, uint32_t obs_m2, int32_t m1_decr, double* base_probp, double* base_lnprobp, double* base_epsp, double* saved_fhetsp, double* saved_fhom1p, double* saved_fhom2p, int32_t* tie_ctp, double *totalp, uint32_t* center_statusp) {
  // similar to plink 1.9 fisher23_tailsum()
  double total = 0;
  double lastp = *base_probp;
  double cur_eps = *base_epsp;
  double hets = *saved_fhetsp;
  double hom1 = *saved_fhom1p;
  double hom2 = *saved_fhom2p;
  double hets_tail;
  double hom1_tail;
  double hom2_tail;
  *center_statusp = 0;
  if (lastp == 0.0) {
    double last_lnp = *base_lnprobp;
    if (last_lnp >= kLnSwitchThresh) {
      while (1) {
        if (hets <= 1) {
          *base_lnprobp = last_lnp;
          *base_epsp = cur_eps;
          *saved_fhetsp = hets;
          *saved_fhom1p = hom1;
          *saved_fhom2p = hom2;
          return 0;
        }
        hom1 += 1;
        hom2 += 1;
        const double lnprob_incr = log((hets * (hets - 1)) / (4 * hom1 * hom2));
        last_lnp += lnprob_incr;
        hets -= 2;
        cur_eps += 3 * k2m52;
        if (lnprob_incr <= -2) {
          cur_eps += (trunc(-lnprob_incr) - 1) * k2m52;
        }
      }
    }
    lastp = exp(last_lnp);
  }
  if (lastp >= 1 + cur_eps) {
    while (1) {
      if (hets <= 1) {
        *saved_fhetsp = hets;
        *saved_fhom1p = hom1;
        *saved_fhom2p = hom2;
        if (lastp < kSwitchThresh) {
          *base_probp = lastp;
          *base_epsp = cur_eps;
        } else {
          *base_probp = 0;
          const double last_lnp = log(lastp);
          *base_lnprobp = last_lnp;
          *base_epsp = cur_eps + ceil_smalleps(last_lnp) * k2m52;
        }
        return 0;
      }
      hom1 += 1;
      hom2 += 1;
      lastp *= (hets * (hets - 1)) / (4 * hom1 * hom2);
      hets -= 2;
      cur_eps += 2 * k2m52;
      if (lastp < 1 + cur_eps) {
        if (lastp <= 1 - cur_eps) {
          break;
        }
        intptr_t cmp_result;
        if (unlikely(HweXchrCompare(obs_fhets, obs_fhom1, obs_fhom2, obs_m1, obs_m2, S_CAST(int32_t, hom1), S_CAST(int32_t, hom2), m1_decr, &starting_lnprobv_ddr, &cmp_result, &lastp))) {
          return 1;
        }
        cur_eps = 3 * k2m52;
        if (cmp_result <= 0) {
          *tie_ctp += (cmp_result == 0);
          break;
        }
      }
    }
    *base_probp = lastp;
    total = lastp;
    hets_tail = hets;
    hom1_tail = hom1;
    hom2_tail = hom2;
  } else {
    hets_tail = hets;
    hom1_tail = hom1;
    hom2_tail = hom2;
    cur_eps += 2 * k2m52;
    const double lastp_tail = lastp;
    uint32_t tie_ct_incr = 0;
    while (1) {
      if (lastp > 1 - cur_eps) {
        if (lastp >= 1 + cur_eps) {
          break;
        }
        intptr_t cmp_result;
        if (unlikely(HweXchrCompare(obs_fhets, obs_fhom1, obs_fhom2, obs_m1, obs_m2, S_CAST(int32_t, hom1), S_CAST(int32_t, hom2), m1_decr, &starting_lnprobv_ddr, &cmp_result, &lastp))) {
          return 1;
        }
        cur_eps = 3 * k2m52;
        if (cmp_result > 0) {
          break;
        }
        tie_ct_incr += (cmp_result == 0);
      }
      total += lastp;
      hets += 2;
      double prob_mult = (4 * hom1 * hom2) / (hets * (hets - 1));
      if (prob_mult <= 1 + 2 * k2m52) {
        // If prob_mult <= 1, we know we're already at a mode for this row,
        // and thus there are no tables in this row with relative-likelihood
        // > 1.
        //
        // There are a few huge cases where 53-bit accuracy isn't enough to
        // distinguish == 1 from > 1; fortunately the function contract
        // ensures int64 is still good enough.
        const int64_t hom1_i = S_CAST(int32_t, hom1);
        const int64_t hom2_i = S_CAST(int32_t, hom2);
        const int64_t hets_i = S_CAST(int32_t, hets);
        const int64_t prob_mult_numer = 4 * hom1_i * hom2_i;
        const int64_t prob_mult_denom = hets_i * (hets_i - 1);
        if (prob_mult_numer <= prob_mult_denom) {
          // Step back to guarantee we stay in range (otherwise hom1 or hom2
          // could become negative).
          hets -= 2;

          // bugfix (25 Apr 2026): When male_ct is more than the ~square of
          // female_ct, the maximum likelihood in row k+1 is no longer
          // guaranteed to be between that of row k and row (k+2).
          // Example: fhets=63, fhom1=27, fhom2=36, m1=8000, m2=9000.
          // - If (from the starting table on row m2=9000) you move to the
          //   fhets+=1, fhom2-=1 adjacent table on row m2=9001, the likelihood
          //   is multiplied by
          //     2 * old_m1 * old_fhom2 / (new_m2 * new_fhets)
          //   = 2 * 8000 * 36 / (9001 * 64) < 1
          // - If you move to the fhets-=1, fhom1+=1 adjacent table,
          //     old_m1 * old_fhets / (2 * new_m2 * new_fhom1)
          //     8000 * 63 / (2 * 9001 * 28) < 1
          // - But if you jump two rows to the fhom1+=1, fhom2-=1 table,
          //     old_m1 * (old_m1-1) * old_fhom2 / ((new_m2-1) * new_m2 * new_fhom1)
          //     8000 * 7999 * 36 / (9001 * 9002 * 28) > 1
          // What's happening here?  If the fhom1 continuous-mode is very close
          // to an integer on one row and far from the nearest integer on the
          // next, it's possible for the likelihood at the continuous-mode to
          // increase from the first row to the second while the integer-mode
          // likelihood decreases.
          //
          // However, if the following additional conditions are true, we can
          // establish center_status=2:
          // 1. If m1_decreasing is false, for each table (fhets, fhom1,
          //    fhom2, m1, m2) in the current or next row, consider the table
          //      (fhets, fhom1-1, fhom2+1, m1+2, m2-2).
          //    That table's likelihood divided by the current table's
          //    likelihood is
          //      m2 * (m2-1) * fhom1 / ((m1+1) * (m1+2) * (fhom2+1))
          //    and if all these ratios <= 1, we know all further-out ratios
          //    are even smaller.
          //    If m1_decreasing true, swap allele 1 and 2 in the expression
          //    above.
          // 2. The next row's mode has relative-likelihood < 1.
          // Alternatively, if the total row relative-likelihood < 1, we can
          // also prove it; the caller is expected to check for this.
          //
          // Both of these conditions are loose.  I'd rather use the
          // observation that the continuous-mode cannot change by more than
          // one fhet per row: so we can compute an upper-bound for the
          // current-row continuous-mode likelihood based on the
          // adjacent-table likelihoods, and then look at gamma-function
          // derivatives.  If digamma(x+1) was equal to log(x+0.5),
          //   m2 * ((fhets/2) + 0.25) <= (m1+1) * (fhom2 + 0.5) and
          //   m2 * (fhom1 + 0.5) <= (m1+1) * ((fhets/2) + 0.25)
          // (along with current continuous-mode relative-likelihood <= 1)
          // would be enough to establish what we need for m1_decreasing false.
          // But that's just an approximation; I haven't yet worked out a
          // provably-correct approach of this type so I will stick to the two
          // loose exit conditions for now.
          const double cur_m1 = u31tod(obs_m1 - m1_decr);
          const double cur_m2 = u31tod(obs_m2 + m1_decr);
          uint32_t potential_exit = 0;
          // hets in between the two tables we just evaluated is guaranteed to
          // be in-range on the next row.
          double next_hets = hets + 1;
          double next_eps = cur_eps + 2 * k2m52;
          double next_hom1;
          double next_hom2;
          double next_m1;
          double next_m2;
          double next_prob;
          // (1) Check extremal likelihood ratios between current-row and
          // row+2, and between row+1 and row+3.
          if (m1_decr < 0) {
            // m1_decreasing is false.
            next_m1 = cur_m1 + 1;
            next_m2 = cur_m2 - 1;
            next_hom1 = hom1 - 1;
            next_hom2 = hom2;
            next_prob = lastp * 2 * cur_m2 * hom1 / (next_m1 * next_hets);
            if (next_prob < 1 - next_eps) {
              // If hom1 > hom2, fhom1/(fhom2+1) is maximized with fhom2=0;
              // otherwise, the limiting case is fhets<2.
              double numer = cur_m2 * (cur_m2 - 1);
              double denom = (cur_m1 + 1) * (cur_m1 + 2);
              if (hom1 > hom2) {
                numer *= hom1 - hom2;
              } else {
                const double hom_incr = trunc(0.5 * hets);
                numer *= hom1 + hom_incr;
                denom *= hom2 + hom_incr;
              }
              if (numer <= denom * (1 - 3 * k2m52)) {
                numer = next_m2 * (next_m2 - 1);
                denom = (next_m1 + 1) * (next_m1 + 2);
                if (next_hom1 > next_hom2) {
                  numer *= next_hom1 - hom2 - 1;
                } else {
                  const double hom_incr = trunc(0.5 * next_hets);
                  numer *= next_hom1 + hom_incr;
                  denom *= next_hom2 + hom_incr;
                }
                potential_exit = (numer <= denom * (1 - 3 * k2m52));
              }
            }
          } else {
            next_m1 = cur_m1 - 1;
            next_m2 = cur_m2 + 1;
            next_hom1 = hom1;
            next_hom2 = hom2 - 1;
            next_prob = lastp * 2 * cur_m1 * hom2 / (next_m2 * next_hets);
            if (next_prob < 1 - next_eps) {
              double numer = cur_m1 * (cur_m1 - 1);
              double denom = (cur_m2 + 1) * (cur_m2 + 2);
              if (hom2 > hom1) {
                numer *= hom2 - hom1;
              } else {
                const double hom_incr = trunc(0.5 * hets);
                numer *= hom2 + hom_incr;
                denom *= hom1 + hom_incr;
              }
              if (numer <= denom * (1 - 3 * k2m52)) {
                numer = next_m1 * (next_m1 - 1);
                denom = (next_m2 + 1) * (next_m2 + 2);
                if (next_hom2 > next_hom1) {
                  numer *= next_hom2 - next_hom1;
                } else {
                  const double hom_incr = trunc(0.5 * next_hets);
                  numer *= next_hom2 + hom_incr;
                  denom *= next_hom1 + hom_incr;
                }
                potential_exit = (numer <= denom * (1 - 3 * k2m52));
              }
            }
          }
          if (potential_exit) {
            // (2) Now check next row's mode.
            const double next_a1_ctd = 2 * next_hom1 + next_hets;
            const double next_a2_ctd = 2 * next_hom2 + next_hets;
            const double next_cmodal_nhet = next_a1_ctd * next_a2_ctd / (next_a1_ctd + next_a2_ctd);
            if (next_hets > next_cmodal_nhet) {
              // Only need to inspect lower het counts.
              while (1) {
                next_hom1 += 1;
                next_hom2 += 1;
                prob_mult = next_hets * (next_hets - 1) / (4 * next_hom1 * next_hom2);
                if (prob_mult < 1) {
                  break;
                }
                next_hets -= 2;
                next_prob *= prob_mult;
                next_eps += 2 * k2m52;
              }
            } else {
              // Only need to inspect higher het counts.
              while (1) {
                next_hets += 2;
                prob_mult = 4 * next_hom1 * next_hom2 / (next_hets * (next_hets - 1));
                if (prob_mult < 1) {
                  break;
                }
                next_hom1 += 1;
                next_hom2 += 1;
                next_prob *= prob_mult;
                next_eps += 2 * k2m52;
              }
            }
            if (next_prob < 1 - next_eps) {
              *center_statusp = 2;
              if (tie_ct_incr) {
                // previous table must have relative-likelihood exactly 1.
                // bugfix (25 Apr 2026): since we won't call this function
                // again for the high-het tail, we must check for its
                // potential tie here.
                *tie_ctp += tie_ct_incr + (prob_mult_numer == prob_mult_denom);
              }
              return 0;
            }
          }
          *center_statusp = 1;
          break;
        }
      }
      lastp *= prob_mult;
      hom1 -= 1;
      hom2 -= 1;
      cur_eps += 2 * k2m52;
    }
    *base_probp = lastp;
    *tie_ctp += tie_ct_incr;
    lastp = lastp_tail;
  }
  *base_epsp = cur_eps;
  *saved_fhetsp = hets;
  *saved_fhom1p = hom1;
  *saved_fhom2p = hom2;
  // {hets,hom1,tmp2}_tail is now at a table with relative-likelihood <= 1, and
  // 'total' is the sum of all relative-likelihoods for tables at least as
  // close to the center.

  // sum tail to floating point precision limit
  hets = hets_tail;
  hom1 = hom1_tail;
  hom2 = hom2_tail;
  while (1) {
    hom1 += 1;
    hom2 += 1;
    lastp *= (hets * (hets - 1)) / (4 * hom1 * hom2);
    hets -= 2;
    const double preaddp = total;
    total += lastp;
    if (total == preaddp) {
      break;
    }
  }
  *totalp = total;
  return 0;
}

// It is assumed that HweXchrLnPLowHetTailsum() was called first for the row.
// In the event that there are no 'center' tables at all, fhets_stop is the
// first table we don't count.
BoolErr HweXchrLnPHighHetTailsum(double fhets_stop, dd_real starting_lnprobv_ddr, uint32_t obs_fhets, uint32_t obs_fhom1, uint32_t obs_fhom2, uint32_t obs_m1, uint32_t obs_m2, int32_t m1_decr, double* base_probp, double* base_lnprobp, double* base_epsp, double* saved_fhetsp, double* saved_fhom1p, double* saved_fhom2p, int32_t* tie_ctp, double *totalp) {
  // similar to plink 1.9 fisher23_tailsum()
  double total = 0;
  double lastp = *base_probp;
  double cur_eps = *base_epsp;
  double hets = *saved_fhetsp;
  double hom1 = *saved_fhom1p;
  double hom2 = *saved_fhom2p;
  double hets_tail;
  double hom1_tail;
  double hom2_tail;
  // identify beginning of tail
  if (lastp == 0.0) {
    double last_lnp = *base_lnprobp;
    if (last_lnp >= kLnSwitchThresh) {
      while (1) {
        const double prev_numer = hom1 * hom2;
        if (prev_numer == 0) {
          // lowest-likelihood table on this side is still too probable
          *base_lnprobp = last_lnp;
          *base_epsp = cur_eps;
          *saved_fhetsp = hets;
          *saved_fhom1p = hom1;
          *saved_fhom2p = hom2;
          return 0;
        }
        hets += 2;
        const double lnprob_incr = log((4 * prev_numer) / (hets * (hets - 1)));
        last_lnp += lnprob_incr;
        hom1 -= 1;
        hom2 -= 1;
        cur_eps += 3 * k2m52;
        if (lnprob_incr <= -2) {
          cur_eps += (trunc(-lnprob_incr) - 1) * k2m52;
        }
        // no risk of last_lnp dropping below cur_eps
      }
    }
    lastp = exp(last_lnp);
  }
  if (lastp >= 1 + cur_eps) {
    while (1) {
      const double prev_numer = hom1 * hom2;
      if (prev_numer == 0) {
        // lowest-likelihood table on this side is still too probable
        *saved_fhetsp = hets;
        *saved_fhom1p = hom1;
        *saved_fhom2p = hom2;
        if (lastp < kSwitchThresh) {
          *base_probp = lastp;
          *base_epsp = cur_eps;
        } else {
          *base_probp = 0;
          const double last_lnp = log(lastp);
          *base_lnprobp = last_lnp;
          *base_epsp = cur_eps + ceil_smalleps(last_lnp) * k2m52;
        }
        return 0;
      }
      hets += 2;
      lastp *= (4 * prev_numer) / (hets * (hets - 1));
      hom1 -= 1;
      hom2 -= 1;
      cur_eps += 2 * k2m52;
      if (lastp < 1 + cur_eps) {
        if (lastp <= 1 - cur_eps) {
          break;
        }
        intptr_t cmp_result;
        if (unlikely(HweXchrCompare(obs_fhets, obs_fhom1, obs_fhom2, obs_m1, obs_m2, S_CAST(int32_t, hom1), S_CAST(int32_t, hom2), m1_decr, &starting_lnprobv_ddr, &cmp_result, &lastp))) {
          return 1;
        }
        cur_eps = 3 * k2m52;
        if (cmp_result <= 0) {
          *tie_ctp += (cmp_result == 0);
          break;
        }
      }
    }
    *base_probp = lastp;
    total = lastp;
    hets_tail = hets;
    hom1_tail = hom1;
    hom2_tail = hom2;
  } else {
    hets_tail = hets;
    hom1_tail = hom1;
    hom2_tail = hom2;
    const double lastp_tail = lastp;
    while (hets > fhets_stop) {
      if (lastp > 1 - cur_eps) {
        if (lastp >= 1 + cur_eps) {
          break;
        }
        intptr_t cmp_result;
        if (unlikely(HweXchrCompare(obs_fhets, obs_fhom1, obs_fhom2, obs_m1, obs_m2, S_CAST(int32_t, hom1), S_CAST(int32_t, hom2), m1_decr, &starting_lnprobv_ddr, &cmp_result, &lastp))) {
          return 1;
        }
        cur_eps = 3 * k2m52;
        if (cmp_result > 0) {
          break;
        }
        *tie_ctp += (cmp_result == 0);
      }
      total += lastp;
      hom1 += 1;
      hom2 += 1;
      lastp *= (hets * (hets - 1)) / (4 * hom1 * hom2);
      hets -= 2;
      cur_eps += 2 * k2m52;
    }
    *base_probp = lastp;
    lastp = lastp_tail;
  }
  *base_epsp = cur_eps;
  *saved_fhetsp = hets;
  *saved_fhom1p = hom1;
  *saved_fhom2p = hom2;
  // {hets,hom1,tmp2}_tail is now at a table with relative-likelihood <= 1, and
  // 'total' is the sum of all relative-likelihoods for tables at least as
  // close to the center.

  // sum tail to floating point precision limit
  hets = hets_tail;
  hom1 = hom1_tail;
  hom2 = hom2_tail;
  while (1) {
    hets += 2;
    lastp *= (4 * hom1 * hom2) / (hets * (hets - 1));
    hom1 -= 1;
    hom2 -= 1;
    const double preaddp = total;
    total += lastp;
    if (total == preaddp) {
      break;
    }
  }
  *totalp = total;
  return 0;
}

// This could optionally take a longer precomputed ddr_lfact table as an
// additional pair of parameters, but it isn't a big deal: large problem
// instances are dominated by ordinary arithmetic operations rather than
// ddr_lfact calls, and the fixed overhead amounts to less than a second of
// wall-clock time when I test millions of chrX variants on an ordinary Mac.
BoolErr HweXchrLnP(int32_t obs_fhets, int32_t obs_fhom1, int32_t obs_fhom2, int32_t obs_m1, int32_t obs_m2, uint32_t midp, double* resultp) {
  // See Graffelman J, Weir BS (2016) Testing for Hardy-Weinberg equilibrium at
  // biallelic genetic markers on the X chromosome.
  // Evaluation strategy is similar to plink 1.9 fisher23().
  if ((!obs_m1) && (!obs_m2)) {
    return HweLnP(obs_fhets, obs_fhom1, obs_fhom2, midp, resultp);
  }
  // 1. Determine relative tail vs. center masses for the starting row.
  //    (Within a row, m1/m2 are constant, so the calculation reduces to the
  //    non-chrX case.  'm' stands for male here.)
  const double a1_ctd = S_CAST(double, obs_fhets + 2LL * obs_fhom1 + obs_m1);
  const double a2_ctd = S_CAST(double, obs_fhets + 2LL * obs_fhom2 + obs_m2);
  // "left" = low hets side, "right" = high hets side
  // orig_base_lnprob{l,r} is valid iff orig_base_prob{l,r} is zero.
  double orig_base_lnprobl = 0;
  double orig_base_lnprobr = 0;
  double orig_base_probl;
  double orig_base_epsl;
  double orig_base_probr;
  double orig_base_epsr;
  double tailp;
  dd_real starting_lnprobv_ddr;
  int32_t tie_ct;
  double orig_saved_lhets;
  double orig_saved_lhom1;
  double orig_saved_lhom2;
  double orig_saved_rhets;
  double orig_saved_rhom1;
  double orig_saved_rhom2;
  if (unlikely(HweLnFirstRow(obs_fhets, obs_fhom1, obs_fhom2, &tailp, &starting_lnprobv_ddr, &tie_ct, &orig_base_probl, &orig_base_lnprobl, &orig_base_epsl, &orig_base_probr, &orig_base_lnprobr, &orig_base_epsr, &orig_saved_lhets, &orig_saved_lhom1, &orig_saved_lhom2, &orig_saved_rhets, &orig_saved_rhom1, &orig_saved_rhom2))) {
    return 1;
  }

  // Returned starting_lnprobv_ddr is for a single row, corresponding to
  // 2^obs_fhets / (obs_fhets! obs_fhom1! obs_fhom2!).
  // Modify this to include the m1! m2! denominator terms that vary between
  // rows.
  starting_lnprobv_ddr = ddr_sub(starting_lnprobv_ddr, ddr_add_lfacts(obs_m1, obs_m2));

  // Other log-factorial expressions we want:
  //
  // * lnprobf, so we can add it to starting_lnprobv at the end and convert
  //   tailp (which is a multiple of starting_prob) into the final
  //   log-[mid]p-value.
  //
  //     (m1+f1)! (m2+f2)! (m1+m2)! (fhets+fhom1+fhom2)!
  //     -----------------------------------------------
  //                     (m1+m2+f1+f2)!
  //
  // * Likelihood of a table is
  //
  //     (m1+f1)! (m2+f2)! (m1+m2)!   (fhets+fhom1+fhom2)!    fhets
  //     -------------------------- * -------------------- * 2
  //       (m1+m2+f1+f2)! m1! m2!     fhets! fhom1! fhom2!
  //
  //   where the first fraction is constant, and the remainder of the
  //   is (f1+f2)! / (f1!f2!) times the likelihood for a female-only HWE test.
  //   Thus, the likelihood sum over an entire row is
  //
  //     (m1+f1)! (m2+f2)! (m1+m2)!   (f1+f2)!
  //     -------------------------- * --------
  //       (m1+m2+f1+f2)! m1! m2!     f1! f2!
  //
  //   We actually want this row-likelihood in starting_prob units:
  //
  //     (m1+f1) (m2+f2)! (m1+m2)! (f1+f2)!    (m1+m2+f1+f2)! / starting_probv
  //     ---------------------------------- * ---------------------------------
  //       (m1+m2+f1+f2)! m1! m2! f1! f2!     (m1+f1)! ... (fhets+fhom1+fhom2)!
  //
  //                          (f1+f2)!
  //   = ---------------------------------------------------
  //     m1! m2! f1! f2! (fhets+fhom1+fhom2)! starting_probv
  const double female_ctd = obs_fhets + obs_fhom1 + obs_fhom2;
  const dd_real female_ct_lfact_ddr = ddr_lfact(female_ctd);
  const dd_real lnprobf_ddr =
    ddr_sub(ddr_add(ddr_add3_lfacts(a1_ctd, a2_ctd, obs_m1 + obs_m2), female_ct_lfact_ddr),
            ddr_lfact(a1_ctd + a2_ctd));
  const dd_real row_relative_lnprobf_ddr =
    ddr_sub(ddr_lfact(2 * female_ctd),
            ddr_add(female_ct_lfact_ddr, starting_lnprobv_ddr));

  // tailp, orig_base_probl, and orig_base_probr are in starting_prob units.
  for (uint32_t m1_decreasing = 0; m1_decreasing != 2; ++m1_decreasing) {
    double m1 = obs_m1;
    double m2 = obs_m2;
    double lhets = orig_saved_lhets;
    double lhom1 = orig_saved_lhom1;
    double lhom2 = orig_saved_lhom2;
    double rhets = orig_saved_rhets;
    double rhom1 = orig_saved_rhom1;
    double rhom2 = orig_saved_rhom2;
    double base_probl = orig_base_probl;
    double base_lnprobl = orig_base_lnprobl;
    double base_epsl = orig_base_epsl;
    double base_probr = orig_base_probr;
    double base_lnprobr = orig_base_lnprobr;
    double base_epsr = orig_base_epsr;
    int32_t m1_decr_last;
    if (m1_decreasing) {
      m1_decr_last = obs_m1;
      const uint32_t uii = 2 * S_CAST(uint32_t, obs_fhom2) + S_CAST(uint32_t, obs_fhets);
      if (uii < S_CAST(uint32_t, obs_m1)) {
        m1_decr_last = uii;
      }
    } else {
      m1_decr_last = -obs_m2;
      const uint32_t uii = 2 * S_CAST(uint32_t, obs_fhom1) + S_CAST(uint32_t, obs_fhets);
      if (uii < S_CAST(uint32_t, obs_m2)) {
        m1_decr_last = -S_CAST(int32_t, uii);
      }
    }
    for (int32_t m1_decr = 0; m1_decr != m1_decr_last; ) {
      m1_decr += m1_decreasing * 2 - 1;
      double prob_mult;
      if (m1_decreasing) {
        const double prev_m1 = m1;
        m2 += 1;
        m1 -= 1;
        // bugfix (19 Apr 2017): We cannot move to the right of the mode here.
        // Otherwise, if the mode itself is more probable than our initial
        // table, but the table to the immediate right of the mode is not,
        // we'll fail to count the mode.
        // ("right" = high het count, "left" = low het count.)
        if (lhets != 0.0) {
          lhom1 += 1;
          prob_mult = (prev_m1 * lhets) / (2 * m2 * lhom1);
          lhets -= 1;
        } else {
          lhets += 1;
          prob_mult = (2 * prev_m1 * lhom2) / (m2 * lhets);
          lhom2 -= 1;
        }
      } else {
        const double prev_m2 = m2;
        m1 += 1;
        m2 -= 1;
        if (lhets != 0.0) {
          lhom2 += 1;
          prob_mult = (prev_m2 * lhets) / (2 * m1 * lhom2);
          lhets -= 1;
        } else {
          lhets += 1;
          prob_mult = (2 * prev_m2 * lhom1) / (m1 * lhets);
          lhom1 -= 1;
        }
      }
      if (base_probl == 0.0) {
        const double lnprob_incr = log(prob_mult);
        base_lnprobl += lnprob_incr;
        base_epsl += 3 * k2m52;
        if (fabs(lnprob_incr) >= 2) {
          base_epsl += (trunc(fabs(lnprob_incr)) - 1) * k2m52;
        }
      } else {
        base_probl *= prob_mult;
        base_epsl += 2 * k2m52;
      }
      double tail_incr = 0.0;
      uint32_t center_status;
      if (unlikely(HweXchrLnPLowHetTailsum(starting_lnprobv_ddr, obs_fhets, obs_fhom1, obs_fhom2, obs_m1, obs_m2, m1_decr, &base_probl, &base_lnprobl, &base_epsl, &lhets, &lhom1, &lhom2, &tie_ct, &tail_incr, &center_status))) {
        return 1;
      }
      if (center_status < 2) {
        if (m1_decreasing) {
          const double prev_m1 = m1 + 1;
          if (rhom2 != 0.0) {
            rhets += 1;
            prob_mult = 2 * prev_m1 * rhom2 / (m2 * rhets);
            rhom2 -= 1;
          } else {
            rhom1 += 1;
            prob_mult = prev_m1 * rhets / (2 * m2 * rhom1);
            rhets -= 1;
          }
        } else {
          const double prev_m2 = m2 + 1;
          if (rhom1 != 0.0) {
            rhets += 1;
            prob_mult = 2 * prev_m2 * rhom1 / (m1 * rhets);
            rhom1 -= 1;
          } else {
            rhom2 += 1;
            prob_mult = prev_m2 * rhets / (2 * m1 * rhom2);
            rhets -= 1;
          }
        }
        if (base_probr == 0.0) {
          const double lnprob_incr = log(prob_mult);
          base_lnprobr += lnprob_incr;
          base_epsr += 3 * k2m52;
          if (fabs(lnprob_incr) >= 2) {
            base_epsr += (trunc(fabs(lnprob_incr)) - 1) * k2m52;
          }
        } else {
          base_probr *= prob_mult;
          base_epsr += 2 * k2m52;
        }
        if (rhets < lhets) {
          // Make sure we don't double-count any tables when center is empty.
          // (may be able to prove this condition never triggers?)
          assert(base_probl != 0.0);
          base_probr = base_probl;
          base_epsr = base_epsl;
          rhets = lhets;
          rhom1 = lhom1;
          rhom2 = lhom2;
        }
        double tail_incr2 = 0.0;
        if (unlikely(HweXchrLnPHighHetTailsum(lhets, starting_lnprobv_ddr, obs_fhets, obs_fhom1, obs_fhom2, obs_m1, obs_m2, m1_decr, &base_probr, &base_lnprobr, &base_epsr, &rhets, &rhom1, &rhom2, &tie_ct, &tail_incr2))) {
          return 1;
        }
        tail_incr += tail_incr2;
        if ((center_status == 1) && (tail_incr < 1.0)) {
          center_status = 2;
        }
      }
      if (center_status == 2) {
        // All tables in subsequent rows have relative-likelihood < 1, all
        // tables in this row have relative-likelihood <= 1 with ties counted.
        double f1 = a1_ctd - m1;
        double f2 = a2_ctd - m2;
        // On the 1000 Genomes dataset, this approach is more expensive than
        // incrementally updating ln_row_prob.  However, it should scale better
        // to larger female_ct.
        const dd_real row_adj_ddr = ddr_add4_lfacts(m1, f1, m2, f2);
        double row_prob = exp(ddr_sub(row_relative_lnprobf_ddr, row_adj_ddr).x[0]);
        if (m1_decreasing) {
          while (1) {
            const double preaddp = tailp;
            tailp += row_prob;
            if (tailp == preaddp) {
              break;
            }
            m2 += 1;
            f1 += 1;
            row_prob *= m1 * f2 / (m2 * f1);
            m1 -= 1;
            f2 -= 1;
          }
        } else {
          while (1) {
            const double preaddp = tailp;
            tailp += row_prob;
            if (tailp == preaddp) {
              break;
            }
            m1 += 1;
            f2 += 1;
            row_prob *= m2 * f1 / (m1 * f2);
            m2 -= 1;
            f1 -= 1;
          }
        }
        break;
      }
      tailp += tail_incr;
    }
  }
  const double starting_lnprob =
    ddr_add(lnprobf_ddr, starting_lnprobv_ddr).x[0];
  if (midp) {
    tailp -= tie_ct * 0.5;
  }
  const double result = log(tailp) + starting_lnprob;
  *resultp = result;
  if (result > -k2m35) {
    // true p-value should always be 1 here
    // (possible todo: check boundary cases with female_ct near 2^31)
    *resultp = 0;
  }
  return 0;
}

#ifdef __cplusplus
}
#endif
