From 18067c163ca461763e9011934f90f0dff1d0a946 Mon Sep 17 00:00:00 2001
From: "Viral B. Shah" <ViralBShah@users.noreply.github.com>
Date: Sat, 18 May 2024 23:33:01 -0400
Subject: [PATCH 1/3] Update to R 4.4.0

---
 RVERSION                  |   2 +-
 include/R_ext/Error.h     |  32 ++++++++---
 include/R_ext/Print.h     |  32 +++++++++--
 include/R_ext/RS.h        |   2 +-
 include/R_ext/libextern.h |   7 ++-
 include/Rmath.h           |  33 ++++++------
 src/choose.c              |  11 ++--
 src/dbinom.c              |  55 ++++++++++++++++---
 src/dnorm.c               |   2 +-
 src/gamma.c               |  22 ++++----
 src/lgammacor.c           |   7 ++-
 src/mlutils.c             |   7 +--
 src/nmath.h               |  10 ++--
 src/pbeta.c               |  17 +++---
 src/pnorm.c               |  14 ++---
 src/qnorm.c               |  63 ++++++++++++++--------
 src/rnchisq.c             |   2 +-
 src/signrank.c            |   4 +-
 src/stirlerr.c            |  70 +++++++++++++++++++-----
 src/toms708.c             | 108 +++++++++++++++++++++-----------------
 src/wilcox.c              |   4 +-
 21 files changed, 330 insertions(+), 174 deletions(-)

diff --git a/RVERSION b/RVERSION
index af8c8ec..fdc6698 100644
--- a/RVERSION
+++ b/RVERSION
@@ -1 +1 @@
-4.2.2
+4.4.0
diff --git a/include/R_ext/Error.h b/include/R_ext/Error.h
index 8d17d37..91a5353 100644
--- a/include/R_ext/Error.h
+++ b/include/R_ext/Error.h
@@ -1,6 +1,6 @@
 /*
  *  R : A Computer Language for Statistical Data Analysis
- *  Copyright (C) 1998-2005   The R Core Team
+ *  Copyright (C) 1998-2023   The R Core Team
  *
  *  This header file is free software; you can redistribute it and/or modify
  *  it under the terms of the GNU Lesser General Public License as published by
@@ -26,21 +26,37 @@
 #ifndef R_ERROR_H_
 #define R_ERROR_H_
 
+#include <R_ext/Print.h>
+
 #ifdef  __cplusplus
 extern "C" {
 #endif
 
-#if defined(__GNUC__) && __GNUC__ >= 3
-#define NORET __attribute__((noreturn))
+/* C23 has a [[noreturn]] attribute supported in GCC 13 and LLVM clang
+ * 15 with -std=c2x but not Apple clang 14.  All have version 202000L.
+ * In C11 there is _Noreturn * (or noreturn in header <stdnoreturn.h>).
+ */
+#if defined NORET
+#elif (defined(__STDC_VERSION__) && __STDC_VERSION__ >= 202301L)
+# define NORET [[noreturn]]
+#elif defined(__STDC_VERSION__) && __STDC_VERSION__ >= 201102L
+# define NORET _Noreturn
+#elif defined(__GNUC__) && __GNUC__ >= 3
+// LLVM and Apple clang identify themselves as 4.
+// But Mandriva (or OpenMandriva) is said to patch clang to 11.
+// Boost also uses this for __SUNPRO_CC >= 0x590
+# define NORET __attribute__((noreturn))
 #else
-#define NORET
+# define NORET
 #endif
 
-void NORET Rf_error(const char *, ...);
-void NORET UNIMPLEMENTED(const char *);
-void NORET WrongArgCount(const char *);
+NORET void Rf_error(const char *, ...) R_PRINTF_FORMAT(1, 2);
+
+NORET void UNIMPLEMENTED(const char *);
+NORET void WrongArgCount(const char *);
+
+void	Rf_warning(const char *, ...) R_PRINTF_FORMAT(1,2);
 
-void	Rf_warning(const char *, ...);
 void 	R_ShowMessage(const char *s);
     
 
diff --git a/include/R_ext/Print.h b/include/R_ext/Print.h
index 2e4cdd3..ff6a3d9 100644
--- a/include/R_ext/Print.h
+++ b/include/R_ext/Print.h
@@ -1,6 +1,6 @@
 /*
  *  R : A Computer Language for Statistical Data Analysis
- *  Copyright (C) 1998-2019    The R Core Team
+ *  Copyright (C) 1998-2024    The R Core Team
  *
  *  This header file is free software; you can redistribute it and/or modify
  *  it under the terms of the GNU Lesser General Public License as published by
@@ -42,11 +42,33 @@ extern "C" {
 # define R_VA_LIST va_list
 #endif
 
-void Rprintf(const char *, ...);
-void REprintf(const char *, ...);
+#ifdef __GNUC__
+# ifdef _WIN32
+#  if defined(_UCRT) || ((__MSVCRT_VERSION__ >= 0x1400) || \
+                        (__MSVCRT_VERSION__ >= 0xE00 && __MSVCRT_VERSION__ < 0x1000))
+#   if defined(__clang__)
+#    define R_PRINTF_FORMAT(M,N) __attribute__ ((format (printf, M, N)))    
+#   else
+#    define R_PRINTF_FORMAT(M,N) __attribute__ ((format (gnu_printf, M, N)))    
+#   endif
+#  else
+#   define R_PRINTF_FORMAT(M,N)
+#  endif
+# else
+#  define R_PRINTF_FORMAT(M,N) __attribute__ ((format (printf, M, N)))
+# endif
+#else
+# define R_PRINTF_FORMAT(M,N)
+#endif
+
+void Rprintf(const char *, ...) R_PRINTF_FORMAT(1, 2);
+void REprintf(const char *, ...) R_PRINTF_FORMAT(1, 2);
+
 #if !defined(__cplusplus) || defined R_USE_C99_IN_CXX
-void Rvprintf(const char *, R_VA_LIST);
-void REvprintf(const char *, R_VA_LIST);
+
+void Rvprintf(const char *, R_VA_LIST) R_PRINTF_FORMAT(1, 0);
+void REvprintf(const char *, R_VA_LIST) R_PRINTF_FORMAT(1, 0);
+
 #endif
 
 #ifdef  __cplusplus
diff --git a/include/R_ext/RS.h b/include/R_ext/RS.h
index 50f3c8b..e993c81 100644
--- a/include/R_ext/RS.h
+++ b/include/R_ext/RS.h
@@ -83,7 +83,7 @@ extern void R_chk_free(void *);
 #define F77_COM(x)     F77_CALL(x)
 #define F77_COMDECL(x) F77_CALL(x)
 
-/* Depreacated in R 2.15.0, non-API
+/* Deprecated in R 2.15.0, non-API
 #if !defined(NO_CALL_R) && defined(DECLARE_LEGACY_CALL_R)
 void	call_R(char*, long, void**, char**, long*, char**, long, char**);
 #endif
diff --git a/include/R_ext/libextern.h b/include/R_ext/libextern.h
index fc3716e..ad16113 100644
--- a/include/R_ext/libextern.h
+++ b/include/R_ext/libextern.h
@@ -1,6 +1,6 @@
 /*
  *  R : A Computer Language for Statistical Data Analysis
- *  Copyright (C) 2001, 2004  The R Core Team.
+ *  Copyright (C) 2001, 2022  The R Core Team.
  *
  *  This header file is free software; you can redistribute it and/or modify
  *  it under the terms of the GNU Lesser General Public License as published by
@@ -32,11 +32,10 @@
 #undef LibImport
 #undef LibExport
 
-/* Don't try to include CYGWIN here: decorating some symbols breaks
-   the auto-export that it relies on, even if R_DLL_BUILD were set. */
 #ifdef _WIN32 /* _WIN32 as does not depend on config.h */
 #define LibImport __declspec(dllimport)
-#define LibExport __declspec(dllexport)
+/* exporting is now done via .def file in R */
+#define LibExport /* __declspec(dllexport) */
 #else
 #define LibImport
 #define LibExport
diff --git a/include/Rmath.h b/include/Rmath.h
index 609eb9f..834c046 100644
--- a/include/Rmath.h
+++ b/include/Rmath.h
@@ -1,6 +1,6 @@
 /* -*- C -*-
  *  Mathlib : A C Library of Special Functions
- *  Copyright (C) 1998-2022  The R Core Team
+ *  Copyright (C) 1998-2024  The R Core Team
  *  Copyright (C) 2004       The R Foundation
  *
  *  This program is free software; you can redistribute it and/or modify
@@ -252,6 +252,7 @@ double  Rlog1p(double);
 #define lgammafn	Rf_lgammafn
 #define lgammafn_sign	Rf_lgammafn_sign
 #define lgamma1p	Rf_lgamma1p
+#define pow1p		Rf_pow1p
 #define log1mexp       	Rf_log1mexp
 #define log1pexp       	Rf_log1pexp
 #define log1pmx		Rf_log1pmx
@@ -394,8 +395,17 @@ double  log1pexp(double); // <-- ../nmath/plogis.c
 double  log1mexp(double);
 double  lgamma1p(double);/* accurate log(gamma(x+1)), small x (0 < x < 0.5) */
 
-double  logspace_add(double, double);
-double  logspace_sub(double, double);
+double  pow1p(double, double); /* pow1p(x, y) := (1+x)^y  accurately also for |x| << 1 */
+
+/* Compute the log of a sum or difference from logs of terms, i.e.,
+ *
+ *     log (exp (logx) + exp (logy))
+ * or  log (exp (logx) - exp (logy))
+ *
+ * without causing overflows or throwing away too much accuracy:
+ */
+double  logspace_add(double logx, double logy);
+double  logspace_sub(double logx, double logy);
 double  logspace_sum(const double *, int);
 
 	/* Beta Distribution */
@@ -426,14 +436,14 @@ double	pnchisq(double, double, double, int, int);
 double	qnchisq(double, double, double, int, int);
 double	rnchisq(double, double);
 
-	/* F Distibution */
+	/* F Distribution */
 
 double	df(double, double, double, int);
 double	pf(double, double, double, int, int);
 double	qf(double, double, double, int, int);
 double	rf(double, double);
 
-	/* Student t Distibution */
+	/* Student t Distribution */
 
 double	dt(double, double, int);
 double	pt(double, double, int, int);
@@ -473,7 +483,7 @@ double	pgeom(double, double, int, int);
 double	qgeom(double, double, int, int);
 double	rgeom(double);
 
-	/* Hypergeometric Distibution */
+	/* Hypergeometric Distribution */
 
 double	dhyper(double, double, double, double, int);
 double	phyper(double, double, double, double, int, int);
@@ -608,17 +618,6 @@ double tanpi(double);
 #endif
 double Rtanpi(double); /* our own in any case */
 
-/* Compute the log of a sum or difference from logs of terms, i.e.,
- *
- *     log (exp (logx) + exp (logy))
- * or  log (exp (logx) - exp (logy))
- *
- * without causing overflows or throwing away too much accuracy:
- */
-double  logspace_add(double logx, double logy);
-double  logspace_sub(double logx, double logy);
-
-
 /* ----------------- Private part of the header file ------------------- */
 
 #if defined(MATHLIB_STANDALONE) && !defined(MATHLIB_PRIVATE_H)
diff --git a/src/choose.c b/src/choose.c
index 3ec9f4f..542b9de 100644
--- a/src/choose.c
+++ b/src/choose.c
@@ -1,7 +1,7 @@
 /*
  *  Mathlib : A C Library of Special Functions
+ *  Copyright (C) 2004-2024 The R Foundation
  *  Copyright (C) 1998      Ross Ihaka
- *  Copyright (C) 2004-2014 The R Foundation
  *
  *  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
@@ -75,8 +75,12 @@ double lchoose(double n, double k)
 #ifndef MATHLIB_STANDALONE
     R_CheckStack();
 #endif
-    if (fabs(k - k0) > 1e-7)
+#define non_INT_WARN_ROUNDING					\
+    /* warn "compatibly" with nmath.h's  R_nonint() : */	\
+    if (fabs(k - k0) > 1e-9 * fmax2(1., fabs(k0)))		\
 	MATHLIB_WARNING2(_("'k' (%.2f) must be integer, rounded to %.0f"), k0, k);
+
+    non_INT_WARN_ROUNDING;
     if (k < 2) {
 	if (k <	 0) return ML_NEGINF;
 	if (k == 0) return 0.;
@@ -118,8 +122,7 @@ double choose(double n, double k)
 #ifndef MATHLIB_STANDALONE
     R_CheckStack();
 #endif
-    if (fabs(k - k0) > 1e-7)
-	MATHLIB_WARNING2(_("'k' (%.2f) must be integer, rounded to %.0f"), k0, k);
+    non_INT_WARN_ROUNDING;
     if (k < k_small_max) {
 	int j;
 	if(n-k < k && n >= 0 && R_IS_INT(n))
diff --git a/src/dbinom.c b/src/dbinom.c
index d2f1cd9..f4020a2 100644
--- a/src/dbinom.c
+++ b/src/dbinom.c
@@ -4,7 +4,9 @@
  *   October 23, 2000.
  *
  *  Merge in to R and further tweaks :
- *	Copyright (C) 2000-2020 The R Core Team
+ *  notably using log1p() and pow1p(), thanks to Morten Welinder, PR#18642
+ *
+ *	Copyright (C) 2000-2024 The R Core Team
  *	Copyright (C) 2008 The R Foundation
  *
  *  This program is free software; you can redistribute it and/or modify
@@ -41,26 +43,63 @@
 #include "nmath.h"
 #include "dpq.h"
 
+/* Compute  (1+x)^y  accurately also for |x| << 1  */
+double pow1p(double x, double y)
+{
+    if(isnan(y))
+	return (x == 0) ? 1. : y; // (0+1)^NaN := 1  by standards
+    if(0 <= y && y == trunc(y) && y <= 4.) {
+	switch((int)y) {
+	case 0: return 1;
+	case 1: return x + 1.;
+	case 2: return x*(x + 2.) + 1.;
+	case 3: return x*(x*(x + 3.) + 3.) + 1.;
+	case 4: return x*(x*(x*(x + 4.) + 6.) + 4.) + 1.;
+	}
+    }
+    /* naive algorithm in two cases: (1) when 1+x is exact (compiler should not over-optimize !),
+     * and (2) when |x| > 1/2 and we have no better algorithm.
+     */
+    if ((x + 1) - 1 == x || fabs(x) > 0.5 || isnan(x))
+	return pow(1 + x, y);
+    else /* not perfect, e.g., for small |x|, non-huge y, use
+	    binom expansion 1 + y*x + y(y-1)/2 x^2 + .. */
+	return exp(y * log1p(x));
+}
+
 double dbinom_raw(double x, double n, double p, double q, int give_log)
 {
     if (p == 0) return((x == 0) ? R_D__1 : R_D__0);
     if (q == 0) return((x == n) ? R_D__1 : R_D__0);
 
-    double lc;
+    // NB: The smaller of p and q is the most accurate
     if (x == 0) {
 	if(n == 0) return R_D__1;
-	lc = (p < 0.1) ? -bd0(n,n*q) - n*p : n*log(q);
-	return( R_D_exp(lc) );
+	if (p > q)
+	    return give_log ? n * log(q)    : pow(q, n);
+	else // 0 < p <= 1/2
+	    return give_log ? n * log1p(-p) : pow1p(-p, n);
     }
-    if (x == n) {
-	lc = (q < 0.1) ? -bd0(n,n*p) - n*q : n*log(p);
-	return( R_D_exp(lc) );
+    if (x == n) { // r = p^x = p^n  -- accurately
+	if (p > q)
+	    return give_log ? n * log1p(-q) : pow1p(-q, n);
+	else
+	    return give_log ? n * log (p)   : pow (p, n);
     }
     if (x < 0 || x > n) return( R_D__0 );
 
+// TODO?  Improve accuracy in these cases:
+#ifdef _NO_LOG_DBINOM_
+    if(!give_log) { // more accurate *not* going via log when result is much much smaller than 1
+	if (x <= M || n-x <= M) { /* use "recursive" direct formula with
+				     k := min(x, n-x) multiplications */
+	}
+    }
+#endif
+
     /* n*p or n*q can underflow to zero if n and p or q are small.  This
        used to occur in dbeta, and gives NaN as from R 2.3.0.  */
-    lc = stirlerr(n) - stirlerr(x) - stirlerr(n-x) - bd0(x,n*p) - bd0(n-x,n*q);
+    double lc = stirlerr(n) - stirlerr(x) - stirlerr(n-x) - bd0(x,n*p) - bd0(n-x,n*q);
 
     /* f = (M_2PI*x*(n-x))/n; could overflow or underflow */
     /* Upto R 2.7.1:
diff --git a/src/dnorm.c b/src/dnorm.c
index c32b9c7..cc113e0 100644
--- a/src/dnorm.c
+++ b/src/dnorm.c
@@ -77,7 +77,7 @@ double dnorm4(double x, double mu, double sigma, int give_log)
      */
     if (x > sqrt(-2*M_LN2*(DBL_MIN_EXP + 1-DBL_MANT_DIG))) return 0.;
 
-    /* Now, to get full accurary, split x into two parts,
+    /* Now, to get full accuracy, split x into two parts,
      *  x = x1+x2, such that |x2| <= 2^-16.
      * Assuming that we are using IEEE doubles, that means that
      * x1*x1 is error free for x<1024 (but we have x < 38.6 anyway).
diff --git a/src/gamma.c b/src/gamma.c
index 3d0f5ee..0a2b854 100644
--- a/src/gamma.c
+++ b/src/gamma.c
@@ -1,7 +1,7 @@
 /*
  *  Mathlib : A C Library of Special Functions
- *  Copyright (C) 2000-2018 The R Core Team
- *  Copyright (C) 2002-2018 The R Foundation
+ *  Copyright (C) 2000-2024 The R Core Team
+ *  Copyright (C) 2002-2024 The R Foundation
  *  Copyright (C) 1998 Ross Ihaka
  *
  *  This program is free software; you can redistribute it and/or modify
@@ -89,10 +89,6 @@ double gammafn(double x)
 	-.5793070335782135784625493333333e-31
     };
 
-    int i, n;
-    double y;
-    double sinpiy, value;
-
 #ifdef NOMORE_FOR_THREADS
     static int ngam = 0;
     static double xmin = 0, xmax = 0., xsml = 0., dxrel = 0.;
@@ -128,7 +124,8 @@ double gammafn(double x)
 	return ML_NAN;
     }
 
-    y = fabs(x);
+    int i;
+    double y = fabs(x), value;
 
     if (y <= 10) {
 
@@ -136,7 +133,7 @@ double gammafn(double x)
 	 * Reduce the interval and find gamma(1 + y) for 0 <= y < 1
 	 * first of all. */
 
-	n = (int) x;
+	int n = (int) x;
 	if(x < 0) --n;
 	y = x - n;/* n = floor(x)  ==>	y in [ 0, 1 ) */
 	--n;
@@ -197,20 +194,21 @@ double gammafn(double x)
 	}
 	else { /* normal case */
 	    value = exp((y - 0.5) * log(y) - y + M_LN_SQRT_2PI +
-			((2*y == (int)2*y)? stirlerr(y) : lgammacor(y)));
+			((2*y == (int)2*y) ? stirlerr(y) : lgammacor(y)));
 	}
+
 	if (x > 0)
 	    return value;
+	// else:  x < 0, not an integer :
 
-	if (fabs((x - (int)(x - 0.5))/x) < dxrel){
-
+	if (fabs((x - (int)(x - 0.5))/x) < dxrel) {
 	    /* The answer is less than half precision because */
 	    /* the argument is too near a negative integer. */
 
 	    ML_WARNING(ME_PRECISION, "gammafn");
 	}
 
-	sinpiy = sinpi(y);
+	double sinpiy = sinpi(y);
 	if (sinpiy == 0) {		/* Negative integer arg - overflow */
 	    ML_WARNING(ME_RANGE, "gammafn");
 	    return ML_POSINF;
diff --git a/src/lgammacor.c b/src/lgammacor.c
index f0addc5..974e115 100644
--- a/src/lgammacor.c
+++ b/src/lgammacor.c
@@ -63,8 +63,6 @@ double attribute_hidden lgammacor(double x)
 	+.1276642195630062933333333333333e-30
     };
 
-    double tmp;
-
 /* For IEEE double precision DBL_EPSILON = 2^-52 = 2.220446049250313e-16 :
  *   xbig = 2 ^ 26.5
  *   xmax = DBL_MAX / 48 =  2^1020 / 3 */
@@ -72,15 +70,16 @@ double attribute_hidden lgammacor(double x)
 #define xbig  94906265.62425156
 #define xmax  3.745194030963158e306
 
-    if (x < 10)
+    if (x < 10) // possibly consider stirlerr()
 	ML_WARN_return_NAN
     else if (x >= xmax) {
 	ML_WARNING(ME_UNDERFLOW, "lgammacor");
 	/* allow to underflow below */
     }
     else if (x < xbig) {
-	tmp = 10 / x;
+	double tmp = 10 / x;
 	return chebyshev_eval(tmp * tmp * 2 - 1, algmcs, nalgm) / x;
     }
+    // else, xbig <= x < xmax :
     return 1 / (x * 12);
 }
diff --git a/src/mlutils.c b/src/mlutils.c
index 0151dd7..43dabfc 100644
--- a/src/mlutils.c
+++ b/src/mlutils.c
@@ -1,6 +1,6 @@
 /*
  *  Mathlib : A C Library of Special Functions
- *  Copyright (C) 1998-2006 The R Core Team
+ *  Copyright (C) 1998-2024 The R Core Team
  *
  *  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
@@ -57,7 +57,8 @@ double R_pow(double x, double y) /* = x ^ y */
 	return(1.);
     if(x == 0.) {
 	if(y > 0.) return(0.);
-	/* y < 0 */return(ML_POSINF);
+	else if(y < 0) return(ML_POSINF);
+	else return(y); /* y is NA or NaN, we assert */
     }
     if (R_FINITE(x) && R_FINITE(y))
 	return(pow(x,y));
@@ -73,7 +74,7 @@ double R_pow(double x, double y) /* = x ^ y */
 	    return((y < 0.)? 0. : ML_POSINF);
 	else {			/* (-Inf) ^ y */
 	    if(R_FINITE(y) && y == floor(y)) /* (-Inf) ^ n */
-		return((y < 0.) ? 0. : (myfmod(y,2.) ? x  : -x));
+		return (y < 0.) ? 0. : (myfmod(y,2.) != 0 ? x  : -x);
 	}
     }
     if(!R_FINITE(y)) {
diff --git a/src/nmath.h b/src/nmath.h
index 5f723ac..46f8e6b 100644
--- a/src/nmath.h
+++ b/src/nmath.h
@@ -1,6 +1,6 @@
 /*
  *  Mathlib : A C Library of Special Functions
- *  Copyright (C) 1998-2022  The R Core Team
+ *  Copyright (C) 1998-2024  The R Core Team
  *
  *  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
@@ -58,8 +58,10 @@ double	Rf_gamma_cody(double);
 #else
 # define R_forceint(x)   round(x)
 #endif
-//R >= 3.1.0:  previously, was defined as  (fabs((x) - R_forceint(x)) > 1e-7)
-# define R_nonint(x) 	  (fabs((x) - R_forceint(x)) > 1e-7*fmax2(1., fabs(x)))
+//R >= 3.1.0; previously: (fabs((x) - R_forceint(x)) > 1e-7)
+//R >= 4.4.0; previously: (fabs((x) - R_forceint(x)) > 1e-7 * fmax2(1., fabs(x)))
+# define R_nonint(x) 	  (fabs((x) - R_forceint(x)) > 1e-9 * fmax2(1., fabs(x)))
+/*						       .... maybe change even to ~ 1e-11 or 12 */
 
 #ifndef MATHLIB_STANDALONE
 
@@ -142,7 +144,7 @@ int R_finite(double);
 #define ME_PRECISION	8
 /*	does not have "full" precision */
 #define ME_UNDERFLOW	16
-/*	and underflow occured (important for IEEE)*/
+/*	and underflow occurred (important for IEEE)*/
 
 
 #define ML_WARN_return_NAN { ML_WARNING(ME_DOMAIN, ""); return ML_NAN; }
diff --git a/src/pbeta.c b/src/pbeta.c
index 98a98f5..28aba9a 100644
--- a/src/pbeta.c
+++ b/src/pbeta.c
@@ -1,6 +1,6 @@
 /*
  *  Mathlib : A C Library of Special Functions
- *  Copyright (C) 2022 The R Core Team
+ *  Copyright (C) 2024 The R Core Team
  *
  *  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
@@ -40,22 +40,24 @@
 attribute_hidden
 double pbeta_raw(double x, double a, double b, int lower_tail, int log_p)
 {
+    if (x >= 1) // may happen when called from qbeta()
+	return R_DT_1;
     // treat limit cases correctly here:
     if(a == 0 || b == 0 || !R_FINITE(a) || !R_FINITE(b)) {
-	// NB:  0 < x < 1 :
+	// NB:  0 <= x < 1 :
 	if(a == 0 && b == 0) // point mass 1/2 at each of {0,1} :
 	    return (log_p ? -M_LN2 : 0.5);
-	if (a == 0 || a/b == 0) // point mass 1 at 0 ==> P(X <= x) = 1, all x > 0
+	if (a == 0 || a/b == 0) // point mass 1 at 0 ==> P(X <= x) = 1, all x >= 0
 	    return R_DT_1;
 	if (b == 0 || b/a == 0) // point mass 1 at 1 ==> P(X <= x) = 0, all x < 1
 	    return R_DT_0;
 	// else, remaining case:  a = b = Inf : point mass 1 at 1/2
 	if (x < 0.5) return R_DT_0; else return R_DT_1;
     }
-    if (x >= 1) // may happen when called from qbeta()
-	return R_DT_1;
+    if (x <= 0)
+	return R_DT_0;
 
-    // Now:  0 < a < Inf;  0 < b < Inf
+    // Now:  0 < a < Inf;  0 < b < Inf  and  0 < x < 1
 
     double x1 = 0.5 - x + 0.5, w, wc;
     int ierr;
@@ -78,8 +80,5 @@ double pbeta(double x, double a, double b, int lower_tail, int log_p)
     if (a < 0 || b < 0) ML_WARN_return_NAN;
     // allowing a==0 and b==0  <==> treat as one- or two-point mass
 
-    if (x <= 0)
-	return R_DT_0;
-
     return pbeta_raw(x, a, b, lower_tail, log_p);
 }
diff --git a/src/pnorm.c b/src/pnorm.c
index 9815a01..5407d70 100644
--- a/src/pnorm.c
+++ b/src/pnorm.c
@@ -1,6 +1,6 @@
 /*
  *  Mathlib : A C Library of Special Functions
- *  Copyright (C) 2000-2020 The R Core Team
+ *  Copyright (C) 2000-2023 The R Core Team
  *  Copyright (C) 2003	    The R Foundation
  *  Copyright (C) 1998	    Ross Ihaka
  *
@@ -144,9 +144,6 @@ void pnorm_both(double x, double *cum, double *ccum, int i_tail, int log_p)
     };
 
     double xden, xnum, temp, del, eps, xsq, y;
-#ifdef NO_DENORMS
-    double min = DBL_MIN;
-#endif
     int i, lower, upper;
 
 #ifdef IEEE_754
@@ -226,9 +223,11 @@ void pnorm_both(double x, double *cum, double *ccum, int i_tail, int log_p)
  * Note that we do want symmetry(0), lower/upper -> hence use y
  */
     else if((log_p && y < 1e170) /* avoid underflow below */
-	/*  ^^^^^ MM FIXME: could speedup for log_p and  |x| >> 5.657 !
+	/*  ^^^^^ MM FIXME: could speed up for log_p and  |x| >> 5.657 !
 	 * Then, make use of  Abramowitz & Stegun, 26.2.13, p.932,  something like
 
+	 * Even smarter: work with   example(pnormAsymp, package="DPQ")
+
 	 xsq = x*x;
 
 	 if(xsq * DBL_EPSILON < 1.)
@@ -240,8 +239,10 @@ void pnorm_both(double x, double *cum, double *ccum, int i_tail, int log_p)
 
  	 swap_tail;
 
-	 Yes, but xsq might be infinite.;
+	 Yes, but xsq might be infinite;
  	 well, actually  x = -1.34..e154 = -sqrt(DBL_MAX) already overflows x^2
+	 The largest x for which  x/2*x is finite is
+	 x = +/- 1.89615038e154 ~= sqrt(2) * sqrt(.Machine$double.xmax)
 	*/
 	    || (lower && -37.5193 < x  &&  x < 8.2924)
 	    || (upper && -8.2924  < x  &&  x < 37.5193)
@@ -268,6 +269,7 @@ void pnorm_both(double x, double *cum, double *ccum, int i_tail, int log_p)
 
 #ifdef NO_DENORMS
     /* do not return "denormalized" -- we do in R */
+    double min = DBL_MIN;
     if(log_p) {
 	if(*cum > -min)	 *cum = -0.;
 	if(*ccum > -min)*ccum = -0.;
diff --git a/src/qnorm.c b/src/qnorm.c
index 1e08280..ec80d36 100644
--- a/src/qnorm.c
+++ b/src/qnorm.c
@@ -1,9 +1,8 @@
 /*
  *  Mathlib : A C Library of Special Functions
- *  Copyright (C) 2000--2020 The R Core Team
+ *  Copyright (C) 2000--2023 The R Core Team
  *  Copyright (C) 1998       Ross Ihaka
- *  based on AS 111 (C) 1977 Royal Statistical Society
- *  and   on AS 241 (C) 1988 Royal Statistical Society
+ *  based on AS 241 (C) 1988 Royal Statistical Society
  *
  *  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
@@ -29,21 +28,17 @@
  *
  *	Compute the quantile function for the normal distribution.
  *
- *	For small to moderate probabilities, algorithm referenced
- *	below is used to obtain an initial approximation which is
- *	polished with a final Newton step.
- *
- *	For very large arguments, an algorithm of Wichura is used.
+ *	The algorithm AS 241 of Wichura is used,
+ *      and has been improved for the very extreme tail (and log_p=TRUE)
  *
  *  REFERENCE
  *
- *	Beasley, J. D. and S. G. Springer (1977).
- *	Algorithm AS 111: The percentage points of the normal distribution,
- *	Applied Statistics, 26, 118-121.
- *
  *      Wichura, M.J. (1988).
  *      Algorithm AS 241: The Percentage Points of the Normal Distribution.
  *      Applied Statistics, 37, 477-484.
+ *
+ *      Maechler, M. (2022). Asymptotic tail formulas for gaussian quantiles;
+ *      https://CRAN.R-project.org/package=DPQ/vignettes/qnorm-asymp.pdf
  */
 
 #include "nmath.h"
@@ -96,13 +91,14 @@ double qnorm5(double p, double mu, double sigma, int lower_tail, int log_p)
     }
     else { /* closer than 0.075 from {0,1} boundary :
 	    *  r := log(p~);  p~ = min(p, 1-p) < 0.075 :  */
+	double lp;
 	if(log_p && ((lower_tail && q <= 0) || (!lower_tail && q > 0))) {
-	    r = p;
+	    lp = p;
 	} else {
-	    r = log( (q > 0) ? R_DT_CIv(p) /* 1-p */ : p_ /* = R_DT_Iv(p) ^=  p */);
+	    lp = log( (q > 0) ? R_DT_CIv(p) /* 1-p */ : p_ /* = R_DT_Iv(p) ^=  p */);
 	}
 	// r = sqrt( - log(min(p,1-p)) )  <==>  min(p, 1-p) = exp( - r^2 ) :
-        r = sqrt(-r);
+        r = sqrt(-lp);
 #ifdef DEBUG_qnorm
 	REprintf("\t close to 0 or 1: r = %7g\n", r);
 #endif
@@ -121,11 +117,11 @@ double qnorm5(double p, double mu, double sigma, int lower_tail, int log_p)
                      r + 1.6763848301838038494) * r +
                     2.05319162663775882187) * r + 1.);
         }
-        else if(r >= 816) { // p is *extremly* close to 0 or 1 - only possibly when log_p =TRUE
-	    // Using the asymptotical formula -- is *not* optimal but uniformly better than branch below
-	    val = r * M_SQRT2;
-        }
-	else { // p is very close to  0 or 1:  r > 5 <==> min(p,1-p) < exp(-25) = 1.3888..e-11
+	else if(r <= 27) { /* p is very close to  0 or 1: r in (5, 27] :
+		*  r >   5 <==> min(p,1-p)  < exp(-25) = 1.3888..e-11
+		*  r <= 27 <==> min(p,1-p) >= exp(-27^2) = exp(-729) ~= 2.507972e-317
+		* i.e., we are just barely in the range where min(p, 1-p) has not yet underflowed to zero.
+		*/
 	    // Wichura, p.478: minimax rational approx R_3(t) is for 5 <= t <= 27  (t :== r)
             r += -5.;
             val = (((((((r * 2.01033439929228813265e-7 +
@@ -141,10 +137,33 @@ double qnorm5(double p, double mu, double sigma, int lower_tail, int log_p)
                      * r + .13692988092273580531) * r +
                     .59983220655588793769) * r + 1.);
         }
-
+        else { // r > 27: p is *really* close to 0 or 1 .. practically only when log_p =TRUE
+	    if(r >= 6.4e8) { // p is *very extremely* close to 0 or 1
+		// Using the asymptotical formula ("0-th order"): qn = sqrt(2*s)
+		val = r * M_SQRT2;
+	    } else {
+		double s2 = -ldexp(lp, 1), // = -2*lp = 2s
+		    x2 = s2 - log(M_2PI * s2); // = xs_1
+		// if(r >= 36000.)  # <==> s >= 36000^2   use x2 = xs_1  above
+		if(r < 36000.) {
+		    x2 = s2 - log(M_2PI * x2) - 2./(2. + x2); // == xs_2
+		    if(r < 840.) { // 27 < r < 840
+			x2 = s2 - log(M_2PI * x2) + 2*log1p(- (1 - 1/(4 + x2))/(2. + x2)); // == xs_3
+			if(r < 109.) { // 27 < r < 109
+			  x2 = s2 - log(M_2PI * x2) +
+			      2*log1p(- (1 - (1 - 5/(6 + x2))/(4. + x2))/(2. + x2)); // == xs_4
+			  if(r < 55.) { // 27 < r < 55
+			    x2 = s2 - log(M_2PI * x2) +
+			      2*log1p(- (1 - (1 - (5 - 9/(8. + x2))/(6. + x2))/(4. + x2))/(2. + x2)); // == xs_5
+			  }
+			}
+		    }
+		}
+                val = sqrt(x2);
+	    }
+	}
 	if(q < 0.0)
 	    val = -val;
-        /* return (q >= 0.)? r : -r ;*/
     }
     return mu + sigma * val;
 }
diff --git a/src/rnchisq.c b/src/rnchisq.c
index e873dee..a7657e6 100644
--- a/src/rnchisq.c
+++ b/src/rnchisq.c
@@ -35,7 +35,7 @@
   mixture of central chisquares with integer degrees of freedom),
   see Formula (29.5b-c) in Johnson, Kotz, Balakrishnan (1995).
 
-  The noncentral chisquare with arbitary degrees of freedom is of interest
+  The noncentral chisquare with arbitrary degrees of freedom is of interest
   for simulating the Cox-Ingersoll-Ross model for interest rates in
   finance.
 
diff --git a/src/signrank.c b/src/signrank.c
index 84289c3..2cb5f5b 100644
--- a/src/signrank.c
+++ b/src/signrank.c
@@ -1,6 +1,6 @@
 /*
  *  Mathlib : A C Library of Special Functions
- *  Copyright (C) 1999-2014  The R Core Team
+ *  Copyright (C) 1999-2024  The R Core Team
  *
  *  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
@@ -124,7 +124,7 @@ double dsignrank(double x, double n, int give_log)
     if (n <= 0)
 	ML_WARN_return_NAN;
 
-    if (fabs(x - R_forceint(x)) > 1e-7)
+    if (R_nonint(x))
 	return(R_D__0);
     x = R_forceint(x);
     if ((x < 0) || (x > (n * (n + 1) / 2)))
diff --git a/src/stirlerr.c b/src/stirlerr.c
index c9b3cc6..7788e19 100644
--- a/src/stirlerr.c
+++ b/src/stirlerr.c
@@ -4,7 +4,7 @@
  *    October 23, 2000.
  *
  *  Merge in to R:
- *	Copyright (C) 2000-2021, The R Core Team
+ *	Copyright (C) 2000-2024, The R Core Team
  *
  *  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
@@ -42,8 +42,14 @@
  *
  * see also lgammacor() in ./lgammacor.c  which computes almost the same!
  *
- * NB: stirlerr(n/2) is called from dt() *and* gamma(n/2) when n is integer and n/2 <= 50
+ * NB: stirlerr(n/2) & stirlerr((n+1)/2) are called from dt(x,n) for all real n > 0 ;
+ *     stirlerr(x) from gammafn(x) when |x| > 10,  2|x| is integer, but |x| is *not* in {11:50}
+ *     stirlerr(x) from dpois_raw(x, lam) for any x > 0  which itself is called by many,
+ *                                        including pgamma(), hence ppois(), ..
+
+ *     stirlerr(n), stirlerr(x), stirlerr(n-x) from binom_raw(x, n, ..) for all possible 0 < x < n
  */
+
 double attribute_hidden stirlerr(double n)
 {
 
@@ -52,6 +58,19 @@ double attribute_hidden stirlerr(double n)
 #define S2 0.00079365079365079365079365  /* 1/1260 */
 #define S3 0.000595238095238095238095238 /* 1/1680 */
 #define S4 0.0008417508417508417508417508/* 1/1188 */
+#define S5 0.0019175269175269175269175262 // 691/360360
+#define S6 0.0064102564102564102564102561 // 1/156
+#define S7 0.029550653594771241830065352  // 3617/122400
+#define S8 0.17964437236883057316493850   // 43867/244188
+#define S9 1.3924322169059011164274315    // 174611/125400
+#define S10 13.402864044168391994478957 // 77683/5796
+#define S11 156.84828462600201730636509 // 236364091/1506960
+#define S12 2193.1033333333333333333333 // 657931/300
+#define S13 36108.771253724989357173269 // 3392780147/93960
+#define S14 691472.26885131306710839498 // 1723168255201/2492028
+#define S15 15238221.539407416192283370 // 7709321041217/505920
+#define S16 382900751.39141414141414141 // 151628697551/396
+/* #define S17 10882266035.784391089015145 // 26315271553053477373/2418179400 */
 
 /*
   exact values for 0, 0.5, 1.0, 1.5, ..., 14.5, 15.0.
@@ -91,16 +110,43 @@ double attribute_hidden stirlerr(double n)
     };
     double nn;
 
-    if (n <= 15.0) {
+    if (n <= 23.5) {
 	nn = n + n;
-	if (nn == (int)nn) return(sferr_halves[(int)nn]);
-	return(lgammafn(n + 1.) - (n + 0.5)*log(n) + n - M_LN_SQRT_2PI);
-    }
+	if (n <= 15. && (nn == (int)nn)) return sferr_halves[(int)nn];
+	// else:
+	if (n <= 5.25) {
+	    if(n >= 1.) { // "MM2"; slightly more accurate than direct form
+		double l_n = log(n);	      // ldexp(u, -1) == u/2
+		return lgamma(n) + n*(1 - l_n) + ldexp(l_n - M_LN_2PI, -1);
+	    }
+	    else // n < 1
+		return lgamma1p(n) - (n + 0.5)*log(n) + n - M_LN_SQRT_2PI;
+	}
+	// else  5.25 < n <= 23.5
+	nn = n*n;
+	if (n > 12.8)	return (S0-(S1-(S2-(S3-(S4-(S5 -S6/nn)/nn)/nn)/nn)/nn)/nn)/n;		// k = 7
+	if (n > 12.3)	return (S0-(S1-(S2-(S3-(S4-(S5-(S6 -S7/nn)/nn)/nn)/nn)/nn)/nn)/nn)/n;	// k = 8
+	if (n >  8.9)	return (S0-(S1-(S2-(S3-(S4-(S5-(S6-(S7 -S8/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/n;	// k = 9
+	/* if (n >  7.9)	return (S0-(S1-(S2-(S3-(S4-(S5-(S6-(S7-(S8 -S9/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/n; skip k = 10 */
+	if (n >  7.3)	return (S0-(S1-(S2-(S3-(S4-(S5-(S6-(S7-(S8-(S9-S10/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/n; // 11
+	/* if (n >  6.5)	return (S0-(S1-(S2-(S3-(S4-(S5-(S6-(S7-(S8-(S9-(S10-S11/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/n; skip k=12*/
+	if (n >  6.6)	return (S0-(S1-(S2-(S3-(S4-(S5-(S6-(S7-(S8-(S9-(S10-(S11-S12/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/n;
+	/* if (n >  5.7)	return (S0-(S1-(S2-(S3-(S4-(S5-(S6-(S7-(S8-(S9-(S10-(S11-(S12-S13/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/n; skip k= 14 */
+	if (n >  6.1)	return (S0-(S1-(S2-(S3-(S4-(S5-(S6-(S7-(S8-(S9-(S10-(S11-(S12-(S13-S14/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/n; // k = 15
+	/* ....		return (S0-(S1-(S2-(S3-(S4-(S5-(S6-(S7-(S8-(S9-(S10-(S11-(S12-(S13-(S14-S15/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/n;
+	 * skip order k=16 : never "good" for double prec */
+	/* 6.1 >= n > 5.25 */
+	return (S0-(S1-(S2-(S3-(S4-(S5-(S6-(S7-(S8-(S9-(S10-(S11-(S12-(S13-(S14-(S15-S16/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/n;
+	/* return (S0-(S1-(S2-(S3-(S4-(S5-(S6-(S7-(S8-(S9-(S10-(S11-(S12-(S13-(S14-(S15-(S16-S17/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/n; */
 
-    nn = n*n;
-    if (n>500) return((S0-S1/nn)/n);
-    if (n> 80) return((S0-(S1-S2/nn)/nn)/n);
-    if (n> 35) return((S0-(S1-(S2-S3/nn)/nn)/nn)/n);
-    /* 15 < n <= 35 : */
-    return((S0-(S1-(S2-(S3-S4/nn)/nn)/nn)/nn)/n);
+     } else { // n > 23.5
+	nn = n*n;
+	if (n > 15.7e6)	return S0/n;
+	if (n > 6180)	return (S0 -S1/nn)/n;
+	if (n > 205)    return (S0-(S1 -S2/nn)/nn)/n;
+	if (n > 86)	return (S0-(S1-(S2 -S3/nn)/nn)/nn)/n;
+	if (n > 27)	return (S0-(S1-(S2-(S3 -S4/nn)/nn)/nn)/nn)/n;
+	/* 23.5 < n <= 27 */
+	return (S0-(S1-(S2-(S3-(S4 -S5/nn)/nn)/nn)/nn)/nn)/n;
+    }
 }
diff --git a/src/toms708.c b/src/toms708.c
index 3136fa5..ab9bdbb 100644
--- a/src/toms708.c
+++ b/src/toms708.c
@@ -5,6 +5,21 @@
    ancillary routines, without investigating the error analysis as we
    do need very high relative accuracy.  This version has about
    14 digits accuracy.
+
+   More specifically,  Brown & Levy (1994) "Certification of Algorithm 708" write
+   "
+    The number of significant digits of accuracy [..] was calculated [..] as
+
+                       - log10 (2 RelativeError),
+    [....]
+    Accuracy ranged from 9.64 significant digits to 15.65 with a
+    median of 14.65 and a lower quartile of 13.81.
+    [...]
+    ... overall accuracy increases slightly as a/b moves away from 1.
+        Linear regression indicates that
+    (1) an average of 13.71 significant digits are obtained in cases in which a = b and
+    (2) the number increases 0.14 significant digits for each unit change in log10(a/b).
+   "
 */
 
 #undef min
@@ -21,7 +36,7 @@
  *
  *	make CFLAGS='-DDEBUG_bratio  ...'
  *MM (w/ Debug, w/o Optimization):
- (cd `R-devel-pbeta-dbg RHOME`/src/nmath ; gcc -I. -I../../src/include -I../../../R/src/include  -DHAVE_CONFIG_H -fopenmp -g -pedantic -Wall --std=gnu99 -DDEBUG_q -DDEBUG_bratio -Wcast-align -Wclobbered  -c ../../../R/src/nmath/toms708.c -o toms708.o; cd ../..; make R)
+ (cd `R-devel-pbeta-dbg RHOME`/src/nmath ; gcc -I. -I../../src/include -I../../../R/src/include  -DHAVE_CONFIG_H -fopenmp  -g -pedantic -Wall -DDEBUG_bratio -DDEBUG_q -Wcast-align -Wconversion -fno-common -Wno-sign-conversion -Wstrict-prototypes -Wclobbered -Werror=implicit-function-declaration  -c ../../../R/src/nmath/toms708.c -o toms708.o; cd ../..; make R)
 */
 #ifdef DEBUG_bratio
 # define R_ifDEBUG_printf(...) REprintf(__VA_ARGS__)
@@ -121,7 +136,7 @@ bratio(double a, double b, double x, double y, double *w, double *w1,
 /*  eps is a machine dependent constant: the smallest
  *      floating point number for which   1. + eps > 1.
  * NOTE: for almost all purposes it is replaced by 1e-15 (~= 4.5 times larger) below */
-    double eps = 2. * Rf_d1mach(3); /* == DBL_EPSILON (in R, Rmath) */
+    double eps = 2. * Rf_d1mach(3); // == DBL_EPSILON (in R, Rmath), but then set to  1e-15  below
 
 /* ----------------------------------------------------------------------- */
     *w  = R_D__0;
@@ -151,9 +166,9 @@ bratio(double a, double b, double x, double y, double *w, double *w1,
     if (a == 0.) goto L211;
     if (b == 0.) goto L201;
 
-    eps = max(eps, 1e-15);
+    eps = max(eps, 1e-15); // = 1e-15 (for IEEE 754)
     Rboolean a_lt_b = (a < b);
-    if (/* max(a,b) */ (a_lt_b ? b : a) < eps * .001) { /* procedure for a and b < 0.001 * eps */
+    if (/* max(a,b) */ (a_lt_b ? b : a) < eps * .001) { /* procedure for a and b < 0.001 * eps = 1e-18 */
 	// L230:  -- result *independent* of x (!)
 	// *w  = a/(a+b)  and  w1 = b/(a+b) :
 	if(log_p) {
@@ -264,7 +279,7 @@ bratio(double a, double b, double x, double y, double *w, double *w1,
     }
     else { /* L30: -------------------- both  a, b > 1  {a0 > 1  &  b0 > 1} ---*/
 
-	/* lambda := a y - b x  =  (a + b)y  =  a - (a+b)x    {using x + y == 1},
+	/* lambda := a y - b x  =  (a + b)y - b  =  a - (a+b)x    {using x + y == 1},
 	 * ------ using the numerically best version : */
 	lambda = R_FINITE(a+b)
 	    ? ((a > b) ? (a + b) * y - b : a - (a + b) * x)
@@ -342,11 +357,11 @@ bratio(double a, double b, double x, double y, double *w, double *w1,
     *w = bup(b0, a0, y0, x0, n, eps, FALSE);
 
     if(*w < DBL_MIN && log_p) { /* do not believe it; try bpser() : */
-	R_ifDEBUG_printf(" L140: bup(b0=%g,..)=%.15g < DBL_MIN - not used; ", b0, *w);
+	R_ifDEBUG_printf(" L140: bup(b0=%g,..)=%.15g < DBL_MIN - not used;\n --> ", b0, *w);
 	/*revert: */ b0 += n;
 	/* which is only valid if b0 <= 1 || b0*x0 <= 0.7 */
 	goto L_w_bpser;
-    }
+    } // else :
     R_ifDEBUG_printf(" L140: *w := bup(b0=%g,..) = %.15g; ", b0, *w);
     if (x0 <= 0.7) {
 	/* log_p :  TODO:  w = bup(.) + bpser(.)  -- not so easy to use log-scale */
@@ -510,36 +525,36 @@ static double bpser(double a, double b, double x, double eps, int log_p)
 /* -----------------------------------------------------------------------
  * Power SERies expansion for evaluating I_x(a,b) when
  *	       b <= 1 or b*x <= 0.7.   eps is the tolerance used.
- * NB: if log_p is TRUE, also use it if   (b < 40  & lambda > 650)
+ * NB: if log_p is TRUE, also use it if   (b < 40  & lambda > 650)   where
+ *                       lambda := a y - b x = (a + b)y - b  = a - (a+b)x   {x + y == 1}
  * ----------------------------------------------------------------------- */
 
-    int i, m;
-    double ans, c, t, u, z, a0, b0, apb;
-
     if (x == 0.) {
 	return R_D__0;
     }
 /* ----------------------------------------------------------------------- */
 /*	      compute the factor  x^a/(a*Beta(a,b)) */
 /* ----------------------------------------------------------------------- */
-    a0 = min(a,b);
+    double ans, c, z, a0 = min(a,b);
     if (a0 >= 1.) { /*		 ------	 1 <= a0 <= b0  ------ */
 	z = a * log(x) - betaln(a, b);
 	ans = log_p ? z - log(a) : exp(z) / a;
     }
     else {
-	b0 = max(a,b);
-
+	double t, u, apb, b0 = max(a,b);
 	if (b0 < 8.) {
 
-	    if (b0 <= 1.) { /*	 ------	 a0 < 1	 and  b0 <= 1  ------ */
+	    if (b0 <= 1.) { /*	------  a0 < 1  and  a0 <= b0 <= 1  ------ */
 
 		if(log_p) {
 		    ans = a * log(x);
 		} else {
 		    ans = pow(x, a);
-		    if (ans == 0.) /* once underflow, always underflow .. */
+		    if (ans == 0.) { /* once underflow, always underflow .. */
+			R_ifDEBUG_printf(" bpser(a=%g, b=%g, x=%g): x^a underflows to 0",
+					 a,b,x);
 			return ans;
+		    }
 		}
 		apb = a + b;
 		if (apb > 1.) {
@@ -558,10 +573,10 @@ static double bpser(double a, double b, double x, double eps, int log_p)
 	    } else { /* 	------	a0 < 1 < b0 < 8	 ------ */
 
 		u = gamln1(a0);
-		m = (int)(b0 - 1.);
+		int m = (int)(b0 - 1.);
 		if (m >= 1) {
 		    c = 1.;
-		    for (i = 1; i <= m; ++i) {
+		    for (int i = 1; i <= m; ++i) {
 			b0 += -1.;
 			c *= b0 / (a0 + b0);
 		    }
@@ -595,11 +610,16 @@ static double bpser(double a, double b, double x, double eps, int log_p)
 		ans = a0 / a * exp(z);
 	}
     }
-    R_ifDEBUG_printf(" bpser(a=%g, b=%g, x=%g, log=%d): prelim.ans = %.14g;\n",
-		     a,b,x, log_p, ans);
+    R_ifDEBUG_printf(" bpser(a=%g, b=%g, x=%g, log=%d, eps=%g): %s = %.14g;",
+		     a,b,x, log_p, eps,
+		     log_p ? "log(x^a/(a*B(a,b)))" : "x^a/(a*B(a,b))",  ans);
     if (ans == R_D__0 || (!log_p && a <= eps * 0.1)) {
+	R_ifDEBUG_printf(" = final answer\n");
 	return ans;
     }
+#ifdef DEBUG_bratio
+    else REprintf("\n");
+#endif
 
 /* ----------------------------------------------------------------------- */
 /*		       COMPUTE THE SERIES */
@@ -622,9 +642,9 @@ static double bpser(double a, double b, double x, double eps, int log_p)
 		" bpser(a=%g, b=%g, x=%g,...) did not converge (n=1e7, |w|/tol=%g > 1; A=%g)",
 		a,b,x, fabs(w)/tol, ans);
     }
-    R_ifDEBUG_printf("  -> n=%.0f iterations, |w|=%g %s %g=tol:=eps/a ==> a*sum=%g\n",
-		     n, fabs(w), (fabs(w) > tol) ? ">!!>" : "<=",
-		     tol, a*sum);
+    R_ifDEBUG_printf("  -> n=%.0f iterations, |w|=%g %s %g=tol:=eps/a ==> a*sum=%g %s -1\n",
+		     n, fabs(w), (fabs(w) > tol) ? ">!!>" : "<=", tol,
+		     a*sum, (a*sum > -1.) ? ">" : "<=");
     if(log_p) {
 	if (a*sum > -1.) ans += log1p(a * sum);
 	else {
@@ -662,7 +682,7 @@ static double bup(double a, double b, double x, double y, int n, double eps,
 	k = (int) exparg(0);
 	if (mu > k)
 	    mu = k;
- 	d = exp(-(double) mu);
+ 	d = exp(-(double) mu); // = exp(-709) = 1.216780751..e-308  nowadays
     }
     else {
 	mu = 0;
@@ -1120,7 +1140,9 @@ static void bgrat(double a, double b, double x, double y, double *w,
 	    "bgrat(a=%g, b=%g, x=%g, y=%g): z=%g, b*z == 0 underflow, hence inaccurate pbeta()",
 	    a,b,x,y, z);
 	/* L_Error:    THE EXPANSION CANNOT BE COMPUTED */
-	 *ierr = 1; return;
+	 *ierr = 1;
+	 /* TODO (?):  do not yet return, trying further using log-scale below */
+	 return;
     }
 
 /*                 COMPUTATION OF THE EXPANSION */
@@ -1519,6 +1541,7 @@ static double rlog1(double x)
 {
 /* -----------------------------------------------------------------------
  *             Evaluation of the function  x - ln(1 + x)
+ *                                        ~=~ - log1pmx(x)  in ./pgamma.c
  * ----------------------------------------------------------------------- */
 
     static double a = .0566749439387324;
@@ -1626,8 +1649,9 @@ static double erfc1(int ind, double x)
 /* ----------------------------------------------------------------------- */
 /*         EVALUATION OF THE COMPLEMENTARY ERROR FUNCTION */
 
-/*          ERFC1(IND,X) = ERFC(X)            IF IND = 0 */
-/*          ERFC1(IND,X) = EXP(X*X)*ERFC(X)   OTHERWISE */
+/*          ERFC1(ind,X) = ERFC(X)            if ind = 0
+ *          ERFC1(ind,X) = EXP(X*X)*ERFC(X)   otherwise, the *scaled* erfc(),
+ *	    				      (the only one used here)     */
 /* ----------------------------------------------------------------------- */
 
     /* Initialized data */
@@ -1682,12 +1706,13 @@ static double erfc1(int ind, double x)
 	    return ret_val;
 	}
 	if (ind == 0 && (x > 100. || x * x > -exparg(1))) {
-	    // LIMIT VALUE FOR LARGE POSITIVE X   WHEN IND = 0
+	    // nowadays: -exparg(1) = 709.0825.. : above <===> |x| > 26.6286
+	    // Underflow to limit for large positive x   when ind = 0
 	    // L60:
 	    return 0.;
 	}
 
-	// L30:
+	// L30:  -5.6 < x < -4  or  4 < x <= 26.6286..
 	t = 1. / (x * x);
 	top = (((r[0] * t + r[1]) * t + r[2]) * t + r[3]) * t + r[4];
 	bot = (((s[0] * t + s[1]) * t + s[2]) * t + s[3]) * t + 1.;
@@ -2010,11 +2035,9 @@ static double psi(double x)
 static double betaln(double a0, double b0)
 {
 /* -----------------------------------------------------------------------
- *     Evaluation of the logarithm of the beta function  ln(beta(a0,b0))
+ *     Evaluation of the logarithm of the beta function  ln(beta(a0,b0))   -- in R:  lbeta(a0, b0)
  * ----------------------------------------------------------------------- */
 
-    static double e = .918938533204673;/* e == 0.5*LN(2*PI) */
-
     double
 	a = min(a0 ,b0),
 	b = max(a0, b0);
@@ -2086,7 +2109,7 @@ static double betaln(double a0, double b0)
 /* ----------------------------------------------------------------------- */
 	// L60:			A >= 8
 /* ----------------------------------------------------------------------- */
-
+	static double e = .918938533204673;/* e == 0.5*LN(2*PI) */
 	double
 	    w = bcorr(a, b),
 	    h = a / b,
@@ -2138,9 +2161,6 @@ static double bcorr(double a0, double b0)
     static double c4 = 8.37308034031215e-4;
     static double c5 = -.00165322962780713;
 
-    /* System generated locals */
-    double ret_val, r1;
-
     /* Local variables */
     double a, b, c, h, t, w, x, s3, s5, x2, s7, s9, s11;
 /* ------------------------ */
@@ -2152,8 +2172,7 @@ static double bcorr(double a0, double b0)
     x = 1. / (h + 1.);
     x2 = x * x;
 
-/*                SET SN = (1 - X^N)/(1 - X) */
-
+/*                SET s<n> := (1 - x^n)/(1 - x) */
     s3 = x + x2 + 1.;
     s5 = x + x2 * s3 + 1.;
     s7 = x + x2 * s5 + 1.;
@@ -2162,21 +2181,14 @@ static double bcorr(double a0, double b0)
 
 /*                SET W = DEL(B) - DEL(A + B) */
 
-/* Computing 2nd power */
-    r1 = 1. / b;
-    t = r1 * r1;
+    t = 1. / b; t *= t; // t := 1 / b^2
     w = ((((c5 * s11 * t + c4 * s9) * t + c3 * s7) * t + c2 * s5) * t + c1 *
 	    s3) * t + c0;
     w *= c / b;
 
 /*                   COMPUTE  DEL(A) + W */
-
-/* Computing 2nd power */
-    r1 = 1. / a;
-    t = r1 * r1;
-    ret_val = (((((c5 * t + c4) * t + c3) * t + c2) * t + c1) * t + c0) / a +
-	    w;
-    return ret_val;
+    t = 1. / a; t *= t; // t:= 1 / a^2
+    return (((((c5 * t + c4) * t + c3) * t + c2) * t + c1) * t + c0) / a + w;
 } /* bcorr */
 
 static double algdiv(double a, double b)
diff --git a/src/wilcox.c b/src/wilcox.c
index 27f0416..785ef6b 100644
--- a/src/wilcox.c
+++ b/src/wilcox.c
@@ -1,6 +1,6 @@
 /*
   Mathlib : A C Library of Special Functions
-  Copyright (C) 1999-2014  The R Core Team
+  Copyright (C) 1999-2024  The R Core Team
 
   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
@@ -174,7 +174,7 @@ double dwilcox(double x, double m, double n, int give_log)
     if (m <= 0 || n <= 0)
 	ML_WARN_return_NAN;
 
-    if (fabs(x - R_forceint(x)) > 1e-7)
+    if (R_nonint(x))
 	return(R_D__0);
     x = R_forceint(x);
     if ((x < 0) || (x > m * n))

From 048bd39a675cda03d9cf348601ca262bc0da130e Mon Sep 17 00:00:00 2001
From: "Viral B. Shah" <ViralBShah@users.noreply.github.com>
Date: Sat, 18 May 2024 23:35:29 -0400
Subject: [PATCH 2/3] Update CI

---
 .github/workflows/CI.yml | 11 +++++++----
 1 file changed, 7 insertions(+), 4 deletions(-)

diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index 7329650..fc558ef 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -17,16 +17,19 @@ jobs:
         version:
           - '1.0'
           - '1'
-          - 'nightly'
         os:
           - ubuntu-latest
+          - windows-latest
         arch:
           - x64
+        include:
+          - os: macOS-14
+            arch: aarch64
     steps:
-      - uses: actions/checkout@v2
-      - uses: julia-actions/setup-julia@v1
+      - uses: actions/checkout@v4
+      - uses: julia-actions/setup-julia@latest
         with:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}
-      - run: make -j 3
+      - run: make -j
       - run: julia test.jl

From ea01888d9276abae98b593a464ee3a6aaf8a6981 Mon Sep 17 00:00:00 2001
From: "Viral B. Shah" <ViralBShah@users.noreply.github.com>
Date: Sat, 18 May 2024 23:36:51 -0400
Subject: [PATCH 3/3] Fix mac aarch64 build

---
 .github/workflows/CI.yml | 1 +
 1 file changed, 1 insertion(+)

diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index fc558ef..7d95cc4 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -25,6 +25,7 @@ jobs:
         include:
           - os: macOS-14
             arch: aarch64
+            version: '1'
     steps:
       - uses: actions/checkout@v4
       - uses: julia-actions/setup-julia@latest