From 9f8871e98ba05e052fe1fc036c529cc4b515282b Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Thu, 23 Oct 2025 14:05:48 -0700 Subject: [PATCH 1/6] src/radmeth/radmeth_model.hpp: fixing the Apple floating point parsing --- src/radmeth/radmeth_model.hpp | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/src/radmeth/radmeth_model.hpp b/src/radmeth/radmeth_model.hpp index 9be0271a..5fe721e6 100644 --- a/src/radmeth/radmeth_model.hpp +++ b/src/radmeth/radmeth_model.hpp @@ -24,6 +24,7 @@ #include #include #include +#include #include template struct mcounts { @@ -130,13 +131,13 @@ Regression::parse(const std::string &line) { field_e = std::find_if(field_s + 1, line_end, is_sep); { const auto [ptr, ec] = std::from_chars(field_s, field_e, mc1.n_reads); - failed = failed || (ptr != field_e); + failed = failed || (ec != std::errc{}); } field_s = std::find_if(field_e + 1, line_end, not_sep); field_e = std::find_if(field_s + 1, line_end, is_sep); - { + if constexpr (std::is_floating_point_v) { #ifdef __APPLE__ const int ret = std::sscanf(field_s, "%lf", &mc1.n_meth); failed = failed || (ret < 1); @@ -145,6 +146,11 @@ Regression::parse(const std::string &line) { failed = failed || (ec != std::errc{}); #endif } + else { + // Apple clang can do std::from_chars for int types + const auto [ptr, ec] = std::from_chars(field_s, field_e, mc1.n_meth); + failed = failed || (ec != std::errc{}); + } mc.push_back(mc1); } From 6fa6c0a7f0459c25bc416335f5d4a40336d64103 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Thu, 23 Oct 2025 19:19:11 -0700 Subject: [PATCH 2/6] src/radmeth/radmeth_nano.cpp and src/radmeth/radmeth_optimize_gamma.hpp: removing linear optimization functions and using only the fdf --- src/radmeth/radmeth_nano.cpp | 7 +++---- src/radmeth/radmeth_optimize_gamma.hpp | 10 ---------- 2 files changed, 3 insertions(+), 14 deletions(-) diff --git a/src/radmeth/radmeth_nano.cpp b/src/radmeth/radmeth_nano.cpp index 0932cb84..367b8e1e 100644 --- a/src/radmeth/radmeth_nano.cpp +++ b/src/radmeth/radmeth_nano.cpp @@ -162,14 +162,13 @@ that the design matrix and the proportion table are correctly formatted. if (has_extreme_counts(t_alt_model)) return std::tuple{1.0, row_status::na_extreme_cnt}; - fit_regression_model_gamma_fdf(t_alt_model, p_estim_alt, - phi_estim_alt); + fit_regression_model_gamma(t_alt_model, p_estim_alt, phi_estim_alt); t_null_model.mc = t_alt_model.mc; t_null_model.rowname = t_alt_model.rowname; - fit_regression_model_gamma_fdf(t_null_model, p_estim_null, - phi_estim_null); + fit_regression_model_gamma(t_null_model, p_estim_null, + phi_estim_null); const double p_value = llr_test(t_null_model.max_loglik, t_alt_model.max_loglik); diff --git a/src/radmeth/radmeth_optimize_gamma.hpp b/src/radmeth/radmeth_optimize_gamma.hpp index 1949b6fb..1f7ee9ad 100644 --- a/src/radmeth/radmeth_optimize_gamma.hpp +++ b/src/radmeth/radmeth_optimize_gamma.hpp @@ -31,14 +31,4 @@ fit_regression_model_gamma(Regression &r, std::vector &p_estimates, double &dispersion_estimate); -void -fit_regression_model_gamma_fdf(Regression &r, - std::vector &p_estimates, - double &dispersion_estimate); - -void -fit_regression_model_gamma_fdf(Regression &r, - std::vector &p_estimates, - double &dispersion_estimate); - #endif // RADMETH_OPTIMIZE_GAMMA_HPP From ec658845738fdeec10e969622c1b1e4feffaad41 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Thu, 23 Oct 2025 19:21:18 -0700 Subject: [PATCH 3/6] src/radmeth/radmeth_optimize_gamma.cpp: removed linear optimization functions and only using the fdf. Reordered some local functions. Added include for radmeth_lgamma and now using that instead of gamma --- src/radmeth/radmeth_optimize_gamma.cpp | 238 +++++++++---------------- 1 file changed, 82 insertions(+), 156 deletions(-) diff --git a/src/radmeth/radmeth_optimize_gamma.cpp b/src/radmeth/radmeth_optimize_gamma.cpp index d4efc75b..533ee81d 100644 --- a/src/radmeth/radmeth_optimize_gamma.cpp +++ b/src/radmeth/radmeth_optimize_gamma.cpp @@ -14,6 +14,7 @@ */ #include "radmeth_optimize_gamma.hpp" +#include "radmeth_lgamma.hpp" #include "radmeth_model.hpp" #include "radmeth_optimize_params.hpp" @@ -22,37 +23,98 @@ #include #include +#include #include #include #include +[[nodiscard]] static inline auto +logistic(const double x) -> double { + return 1.0 / (1.0 / std::exp(x) + 1.0); +} + +template +[[nodiscard]] static auto +get_p(const std::vector &v, const gsl_vector *params) -> double { + const auto a = v.data(); + return logistic(std::inner_product(a, a + std::size(v), params->data, 0.0)); +} + +static auto +get_cache_lgamma(const std::vector &group, + const gsl_vector *params, const double phi, + vars_cache &cache) { + const auto p = get_p(group, params); + const auto a = p * phi; + const auto b = (1.0 - p) * phi; + cache.p = p; + cache.a = a; + cache.b = b; + cache.lgamma_a = radmeth_lgamma::noneg_lgamma(a); + cache.lgamma_b = radmeth_lgamma::noneg_lgamma(b); + cache.lgamma_a_b = radmeth_lgamma::noneg_lgamma(a + b); +} + +template +[[nodiscard]] static auto +log_likelihood(const gsl_vector *params, Regression ®) -> double { + const auto phi = 1.0 / std::exp(gsl_vector_get(params, reg.n_factors())); + + const auto n_groups = reg.n_groups(); + const auto &groups = reg.design.groups; + auto &cache = reg.cache; + for (auto g_idx = 0u; g_idx < n_groups; ++g_idx) + get_cache_lgamma(groups[g_idx], params, phi, cache[g_idx]); + + const auto &group_id = reg.design.group_id; + const auto &mc = reg.mc; + double ll = 0.0; + const auto n_samples = reg.n_samples(); + for (auto i = 0u; i < n_samples; ++i) { + const auto y = mc[i].n_meth; + const auto n = mc[i].n_reads; + + const auto &c = cache[group_id[i]]; + const auto a = c.a; + const auto b = c.b; + + // clang-format off + ll += ((radmeth_lgamma::noneg_lgamma(y + a) - c.lgamma_a) + + (radmeth_lgamma::noneg_lgamma(n - y + b) - c.lgamma_b) + + (c.lgamma_a_b - radmeth_lgamma::noneg_lgamma(n + a + b))); + // clang-format on + } + + return ll; +} + /* Coefficients for the Chebyschev polynomial for the digamma function in the range 0-1 */ // clang-format off static constexpr std::array psi1_cs { -0.038057080835217922, // == -0.019028540417608961*2 - 0.491415393029387130, + +0.491415393029387130, -0.056815747821244730, - 0.008357821225914313, + +0.008357821225914313, -0.001333232857994342, - 0.000220313287069308, + +0.000220313287069308, -0.000037040238178456, - 0.000006283793654854, + +0.000006283793654854, -0.000001071263908506, - 0.000000183128394654, + +0.000000183128394654, -0.000000031353509361, - 0.000000005372808776, + +0.000000005372808776, -0.000000000921168141, - 0.000000000157981265, + +0.000000000157981265, -0.000000000027098646, - 0.000000000004648722, + +0.000000000004648722, -0.000000000000797527, - 0.000000000000136827, + +0.000000000000136827, -0.000000000000023475, - 0.000000000000004027, + +0.000000000000004027, -0.000000000000000691, - 0.000000000000000118, - -0.000000000000000020 + +0.000000000000000118, + -0.000000000000000020, }; // clang-format on @@ -62,20 +124,20 @@ static constexpr std::array psi1_cs { static constexpr std::array psi2_cs = { -0.0204749044678185, // == -0.01023745223390925*2 -0.0101801271534859, - 0.0000559718725387, + +0.0000559718725387, -0.0000012917176570, - 0.0000000572858606, + +0.0000000572858606, -0.0000000038213539, - 0.0000000003397434, + +0.0000000003397434, -0.0000000000374838, - 0.0000000000048990, + +0.0000000000048990, -0.0000000000007344, - 0.0000000000001233, + +0.0000000000001233, -0.0000000000000228, - 0.0000000000000045, + +0.0000000000000045, -0.0000000000000009, - 0.0000000000000002, - -0.0000000000000000 + +0.0000000000000002, + -0.0000000000000000, }; // clang-format on @@ -108,66 +170,6 @@ digamma(const double y) -> double { return chebyschev(psi1_cs, psi1_order, 2.0 * (y - 1.0) - 1.0); } -[[nodiscard]] static inline auto -logistic(const double x) -> double { - return 1.0 / (1.0 / std::exp(x) + 1.0); -} - -template -[[nodiscard]] static auto -get_p(const std::vector &v, const gsl_vector *params) -> double { - const auto a = v.data(); - return logistic(std::inner_product(a, a + std::size(v), params->data, 0.0)); -} - -static auto -get_cache_lgamma(const std::vector &group, - const gsl_vector *params, const double phi, - vars_cache &cache) { - const auto p = get_p(group, params); - const auto a = p * phi; - const auto b = (1.0 - p) * phi; - cache.p = p; - cache.a = a; - cache.b = b; - cache.lgamma_a = std::lgamma(a); - cache.lgamma_b = std::lgamma(b); - cache.lgamma_a_b = std::lgamma(a + b); -} - -template -[[nodiscard]] static auto -log_likelihood(const gsl_vector *params, Regression ®) -> double { - const auto phi = 1.0 / std::exp(gsl_vector_get(params, reg.n_factors())); - - const auto n_groups = reg.n_groups(); - const auto &groups = reg.design.groups; - auto &cache = reg.cache; - for (auto g_idx = 0u; g_idx < n_groups; ++g_idx) - get_cache_lgamma(groups[g_idx], params, phi, cache[g_idx]); - - const auto &group_id = reg.design.group_id; - const auto &mc = reg.mc; - double ll = 0.0; - const auto n_samples = reg.n_samples(); - for (auto i = 0u; i < n_samples; ++i) { - const auto y = mc[i].n_meth; - const auto n = mc[i].n_reads; - - const auto &c = cache[group_id[i]]; - const auto a = c.a; - const auto b = c.b; - - // clang-format off - ll += ((std::lgamma(y + a) - c.lgamma_a) + - (std::lgamma(n - y + b) - c.lgamma_b) + - (c.lgamma_a_b - std::lgamma(n + a + b))); - // clang-format on - } - - return ll; -} - static void get_cache_digamma(const std::vector &group, const gsl_vector *params, const double phi, @@ -274,68 +276,6 @@ fit_regression_model(Regression &r, std::vector &p_estimates, const auto n_groups = r.n_groups(); r.cache.resize(n_groups); // make sure scratch space is allocated - const std::size_t n_params = r.n_params(); - const auto tol = radmeth_optimize_params::tolerance; - - // set the parameters: zero for "p" parameters and the final one for - // dispersion using the constant - auto params = gsl_vector_alloc(n_params); - gsl_vector_set_all(params, 0.0); - gsl_vector_set(params, n_params - 1, init_dispersion_param); - - auto minex_func = gsl_multimin_function{ - &neg_loglik, - n_params, - static_cast(&r), - }; - - // set initial step size for all dims - auto step_sizes = gsl_vector_alloc(n_params); - gsl_vector_set_all(step_sizes, stepsize); - - const auto minimizer = gsl_multimin_fminimizer_nmsimplex2; - auto s = gsl_multimin_fminimizer_alloc(minimizer, n_params); - gsl_multimin_fminimizer_set(s, &minex_func, params, step_sizes); - - int status = 0; - std::size_t iter = 0; - - do { - status = gsl_multimin_fminimizer_iterate(s); // one iter and get status - if (status) - break; - - const auto size = gsl_multimin_fminimizer_size(s); - status = gsl_multimin_test_size(size, tol); - } while (status == GSL_CONTINUE && ++iter < max_iter); - // ADS: can't use (status != GSL_SUCCESS) - - const auto param_estimates = gsl_multimin_fminimizer_x(s); - - const auto &groups = r.design.groups; - p_estimates.clear(); - for (auto g_idx = 0u; g_idx < n_groups; ++g_idx) - p_estimates.push_back(get_p(groups[g_idx], param_estimates)); - const auto disp_param = gsl_vector_get(param_estimates, n_params - 1); - dispersion_estimate = 1.0 / std::exp(disp_param); - - r.max_loglik = log_likelihood(param_estimates, r); - - gsl_multimin_fminimizer_free(s); - gsl_vector_free(step_sizes); - gsl_vector_free(params); -} - -template -static void -fit_regression_model_fdf(Regression &r, std::vector &p_estimates, - double &dispersion_estimate) { - static constexpr auto init_dispersion_param = -2.5; - const auto stepsize = radmeth_optimize_params::stepsize; - const auto max_iter = radmeth_optimize_params::max_iter; - const auto n_groups = r.n_groups(); - r.cache.resize(n_groups); // make sure scratch space is allocated - const std::size_t n_params = r.n_params(); const auto tol = std::sqrt(n_params) * r.n_samples() * radmeth_optimize_params::tolerance; @@ -404,17 +344,3 @@ fit_regression_model_gamma(Regression &r, double &dispersion_estimate) { fit_regression_model(r, p_estimates, dispersion_estimate); } - -void -fit_regression_model_gamma_fdf(Regression &r, - std::vector &p_estimates, - double &dispersion_estimate) { - fit_regression_model_fdf(r, p_estimates, dispersion_estimate); -} - -void -fit_regression_model_gamma_fdf(Regression &r, - std::vector &p_estimates, - double &dispersion_estimate) { - fit_regression_model_fdf(r, p_estimates, dispersion_estimate); -} From 725c878eb5e5d1e222cef34b405bab8e5997856b Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Thu, 23 Oct 2025 19:57:37 -0700 Subject: [PATCH 4/6] src/radmeth/radmeth_lgamma.hpp: adding an lgamma implementation because the one in glibc is slow and I don't think we need all the branching, etc., which seemed to be pressuring icache --- src/radmeth/radmeth_lgamma.hpp | 266 +++++++++++++++++++++++++++++++++ 1 file changed, 266 insertions(+) create mode 100644 src/radmeth/radmeth_lgamma.hpp diff --git a/src/radmeth/radmeth_lgamma.hpp b/src/radmeth/radmeth_lgamma.hpp new file mode 100644 index 00000000..c719c104 --- /dev/null +++ b/src/radmeth/radmeth_lgamma.hpp @@ -0,0 +1,266 @@ +/* Copyright (C) 2025 Andrew D. + * + * Author: Andrew D. Smith + * + * This program is free software: you can redistribute it and/or modify it + * under the terms of the GNU General Public License as published by the Free + * Software Foundation, either version 3 of the License, or (at your option) + * any later version. + * + * This program 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 General Public License for + * more details. + */ + +#ifndef RADMETH_LGAMMA_HPP +#define RADMETH_LGAMMA_HPP + +/* ADS: the constants below, and some of the code, was taken from e_lgamma_r.c + in Google's bionic source and lgamma_r.c in the MUSL source. One or both + were modified from the original, maybe for style guidlines. Both files had + this header. + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunSoft, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== +*/ + +#include +#include +#include + +namespace radmeth_lgamma { + +// ADS: constant 'tc' below is the unique minimum of the gamma function on R+ +static constexpr std::double_t tc = +1.46163214496836224576e+00; +static constexpr std::double_t tf = -1.21486290535849611461e-01; +/* tt = -(tail of tf) */ +static constexpr std::double_t tt = -3.63867699703950536541e-18; + +// clang-format off + +static constexpr std::double_t A0[6] = { + +7.72156649015328655494e-02, + +6.73523010531292681824e-02, + +7.38555086081402883957e-03, + +1.19270763183362067845e-03, + +2.20862790713908385557e-04, + +2.52144565451257326939e-05, +}; + +static constexpr std::double_t A1[6] = { + +3.22467033424113591611e-01, + +2.05808084325167332806e-02, + +2.89051383673415629091e-03, + +5.10069792153511336608e-04, + +1.08011567247583939954e-04, + +4.48640949618915160150e-05, +}; + +static constexpr std::double_t T0[5] = { + +4.83836122723810047042e-01, + -3.27885410759859649565e-02, + +6.10053870246291332635e-03, + -1.40346469989232843813e-03, + +3.15632070903625950361e-04, +}; + +static constexpr std::double_t T1[5] = { + -1.47587722994593911752e-01, + +1.79706750811820387126e-02, + -3.68452016781138256760e-03, + +8.81081882437654011382e-04, + -3.12754168375120860518e-04, +}; + +static constexpr std::double_t T2[5] = { + +6.46249402391333854778e-02, + -1.03142241298341437450e-02, + +2.25964780900612472250e-03, + -5.38595305356740546715e-04, + +3.35529192635519073543e-04, +}; + +static constexpr std::double_t U[6] = { + -7.72156649015328655494e-02, + +6.32827064025093366517e-01, + +1.45492250137234768737e+00, + +9.77717527963372745603e-01, + +2.28963728064692451092e-01, + +1.33810918536787660377e-02, +}; + +static constexpr std::double_t V[5] = { + +2.45597793713041134822e+00, + +2.12848976379893395361e+00, + +7.69285150456672783825e-01, + +1.04222645593369134254e-01, + +3.21709242282423911810e-03, +}; + +static constexpr std::double_t S[7] = { + -7.72156649015328655494e-02, + +2.14982415960608852501e-01, + +3.25778796408930981787e-01, + +1.46350472652464452805e-01, + +2.66422703033638609560e-02, + +1.84028451407337715652e-03, + +3.19475326584100867617e-05, +}; + +static constexpr std::double_t R[6] = { + +1.39200533467621045958e+00, + +7.21935547567138069525e-01, + +1.71933865632803078993e-01, + +1.86459191715652901344e-02, + +7.77942496381893596434e-04, + +7.32668430744625636189e-06, +}; + +static constexpr std::double_t W[7] = { + +4.18938533204672725052e-01, + +8.33333333333329678849e-02, + -2.77777777728775536470e-03, + +7.93650558643019558500e-04, + -5.95187557450339963135e-04, + +8.36339918996282139126e-04, + -1.63092934096575273989e-03, +}; + +// clang-format on + +[[nodiscard]] inline auto +noneg_lgamma(double x) -> double { + // ADS: this function does not work for very small or very large values. + assert(std::isfinite(x) && x > 2.2e-16); + + union { + double f; // unused + std::uint64_t as_bits; + } u = {x}; // safe inspection of binary layout + + // get magnitude + const std::uint32_t ix = u.as_bits >> 32 & 0x7fffffff; + + std::double_t r{}; + // x exactly equal to 1 or 2 + if ((ix == 0x3ff00000 || ix == 0x40000000) && + static_cast(u.as_bits) == 0) + return 0; + + // ADS: for x > 2^58 the code below applies. Not using it because this is + // only needed for args that would be far from optimal in param estimates. + + // if (ix > 0x43900000) + // return x * (std::log(x) - 1.0); + + // x < 2.0 + if (ix < 0x40000000) { + int i{}; + std::double_t y{}; + if (ix <= 0x3feccccc) { /* lgamma(x) = lgamma(x+1)-log(x) */ + r = -std::log(x); + if (ix >= 0x3FE76944) { + y = 1.0 - x; + i = 0; + } + else if (ix >= 0x3FCDA661) { + y = x - (tc - 1.0); + i = 1; + } + else { + y = x; + i = 2; + } + } + else { + r = 0.0; + if (ix >= 0x3FFBB4C3) { /* [1.7316,2] */ + y = 2.0 - x; + i = 0; + } + else if (ix >= 0x3FF3B4C4) { /* [1.23,1.73] */ + y = x - tc; + i = 1; + } + else { + y = x - 1.0; + i = 2; + } + } + std::double_t p{}; + std::double_t w{}, z{}; + std::double_t p1{}, p2{}, p3{}; + switch (i) { + case 0: + z = y * y; + p1 = A0[0] + + z * (A0[1] + z * (A0[2] + z * (A0[3] + z * (A0[4] + z * A0[5])))); + p2 = + z * (A1[0] + + z * (A1[1] + z * (A1[2] + z * (A1[3] + z * (A1[4] + z * A1[5]))))); + p = y * p1 + p2; + r += -0.5 * y + p; + break; + case 1: + z = y * y; + w = z * y; + p1 = T0[0] + w * (T0[1] + w * (T0[2] + w * (T0[3] + w * T0[4]))); + p2 = T1[0] + w * (T1[1] + w * (T1[2] + w * (T1[3] + w * T1[4]))); + p3 = T2[0] + w * (T2[1] + w * (T2[2] + w * (T2[3] + w * T2[4]))); + p = z * p1 - (tt - w * (p2 + y * p3)); + r += tf + p; + break; + case 2: + p1 = y * (U[0] + + y * (U[1] + y * (U[2] + y * (U[3] + y * (U[4] + y * U[5]))))); + p2 = 1.0 + y * (V[0] + y * (V[1] + y * (V[2] + y * (V[3] + y * V[4])))); + r += -0.5 * y + p1 / p2; + } + return r; + } + + // x < 8.0 + if (ix < 0x40200000) { + const int i = static_cast(x); + const std::double_t y = x - static_cast(i); + const std::double_t p = + y * (S[0] + + y * (S[1] + + y * (S[2] + y * (S[3] + y * (S[4] + y * (S[5] + y * S[6])))))); + const std::double_t q = + 1.0 + + y * (R[0] + y * (R[1] + y * (R[2] + y * (R[3] + y * (R[4] + y * R[5]))))); + r = 0.5 * y + p / q; // RF P/Q + std::double_t z = 1.0; // lgamma(1+s) = log(s) + lgamma(s) + // clang-format off + switch (i) { + case 7: z *= y + 6.0; // FALLTHRU + case 6: z *= y + 5.0; // FALLTHRU + case 5: z *= y + 4.0; // FALLTHRU + case 4: z *= y + 3.0; // FALLTHRU + case 3: z *= y + 2.0; // FALLTHRU + r += std::log(z); + break; + } + // clang-format on + return r; + } + + // 8.0 <= x < 2^58 (ix < 0x43900000) + const std::double_t t = std::log(x); + const std::double_t z = 1.0 / x; + const std::double_t y = z * z; + const std::double_t w = + W[0] + + z * (W[1] + y * (W[2] + y * (W[3] + y * (W[4] + y * (W[5] + y * W[6]))))); + return (x - 0.5) * (t - 1.0) + w; +} +}; // namespace radmeth_lgamma + +#endif // RADMETH_LGAMMA_HPP From decb1f71530a194d3edc38aac28a56737cda9a3b Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Fri, 24 Oct 2025 09:17:40 -0700 Subject: [PATCH 5/6] src/radmeth/radmeth_lgamma.hpp: minor cleanup --- src/radmeth/radmeth_lgamma.hpp | 25 ++++++++++++++++--------- 1 file changed, 16 insertions(+), 9 deletions(-) diff --git a/src/radmeth/radmeth_lgamma.hpp b/src/radmeth/radmeth_lgamma.hpp index c719c104..d7332654 100644 --- a/src/radmeth/radmeth_lgamma.hpp +++ b/src/radmeth/radmeth_lgamma.hpp @@ -36,7 +36,7 @@ namespace radmeth_lgamma { -// ADS: constant 'tc' below is the unique minimum of the gamma function on R+ +// ADS: 'tc' is the unique minimum of the gamma function on R+ static constexpr std::double_t tc = +1.46163214496836224576e+00; static constexpr std::double_t tf = -1.21486290535849611461e-01; /* tt = -(tail of tf) */ @@ -137,7 +137,7 @@ static constexpr std::double_t W[7] = { [[nodiscard]] inline auto noneg_lgamma(double x) -> double { // ADS: this function does not work for very small or very large values. - assert(std::isfinite(x) && x > 2.2e-16); + assert(std::isfinite(x) && x > 2.22045e-16); // > 2^-52 union { double f; // unused @@ -153,17 +153,24 @@ noneg_lgamma(double x) -> double { static_cast(u.as_bits) == 0) return 0; - // ADS: for x > 2^58 the code below applies. Not using it because this is - // only needed for args that would be far from optimal in param estimates. + // ADS: for x < 2.22045e-16, just return the value of lgamma at + // that 2.22045e-16 + + // if (ix < 0x3CB00000) + // return 36.0437; + + // ADS: for x > 2^58 the code below applies. Maxing it out at 1.2994e+19, + // because that's large enough and if we see this our estimation is going in + // the wrong direction anyway. // if (ix > 0x43900000) - // return x * (std::log(x) - 1.0); + // return 1.12994e+19; // proper way: return x * (std::log(x) - 1.0); // x < 2.0 if (ix < 0x40000000) { int i{}; std::double_t y{}; - if (ix <= 0x3feccccc) { /* lgamma(x) = lgamma(x+1)-log(x) */ + if (ix <= 0x3feccccc) { // lgamma(x) = lgamma(x+1)-log(x) r = -std::log(x); if (ix >= 0x3FE76944) { y = 1.0 - x; @@ -180,11 +187,11 @@ noneg_lgamma(double x) -> double { } else { r = 0.0; - if (ix >= 0x3FFBB4C3) { /* [1.7316,2] */ + if (ix >= 0x3FFBB4C3) { // [1.7316,2] y = 2.0 - x; i = 0; } - else if (ix >= 0x3FF3B4C4) { /* [1.23,1.73] */ + else if (ix >= 0x3FF3B4C4) { // [1.23,1.73] y = x - tc; i = 1; } @@ -252,7 +259,7 @@ noneg_lgamma(double x) -> double { return r; } - // 8.0 <= x < 2^58 (ix < 0x43900000) + // 8.0 <= x (ix < 0x43900000) const std::double_t t = std::log(x); const std::double_t z = 1.0 / x; const std::double_t y = z * z; From 922d720515e254558a54b227855b59d76670ffb2 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Fri, 24 Oct 2025 09:20:09 -0700 Subject: [PATCH 6/6] Makefile.am: adding the radmeth_lgamma.hpp source --- Makefile.am | 1 + 1 file changed, 1 insertion(+) diff --git a/Makefile.am b/Makefile.am index 58a1c915..bb85391d 100644 --- a/Makefile.am +++ b/Makefile.am @@ -193,6 +193,7 @@ libdnmtools_a_SOURCES += \ src/radmeth/radmeth_optimize_gamma.hpp \ src/radmeth/radmeth_optimize_params.hpp \ src/radmeth/radmeth_model.hpp \ + src/radmeth/radmeth_lgamma.hpp \ src/radmeth/radmeth_utils.hpp \ src/radmeth/radmeth_design.hpp