From e54f12ddc60e96e0386c079b1616a6ec31b409f1 Mon Sep 17 00:00:00 2001
From: Andreas Noack <andreas@noack.dk>
Date: Tue, 10 Sep 2024 11:21:06 +0200
Subject: [PATCH 1/2] Update to version 4.4.1

---
 RVERSION               |  2 +-
 src/qDiscrete_search.h | 28 +++++++++++++++++-----------
 src/qpois.c            |  2 +-
 3 files changed, 19 insertions(+), 13 deletions(-)

diff --git a/RVERSION b/RVERSION
index fdc6698..cca25a9 100644
--- a/RVERSION
+++ b/RVERSION
@@ -1 +1 @@
-4.4.0
+4.4.1
diff --git a/src/qDiscrete_search.h b/src/qDiscrete_search.h
index 0bdbcb4..df7bdf1 100644
--- a/src/qDiscrete_search.h
+++ b/src/qDiscrete_search.h
@@ -1,6 +1,6 @@
 /*
  *  Mathlib : A C Library of Special Functions
- *  Copyright (C) 2000-2021 The R Core Team
+ *  Copyright (C) 2000-2024 The R Core Team
  *  Copyright (C) 2005-2021 The R Foundation
  *
  *  This program is free software; you can redistribute it and/or modify
@@ -18,7 +18,7 @@
  *  https://www.R-project.org/Licenses/
  */
 
-/* This is #included  from ./qnbinom.c , .........
+/* This is #included  from ./qpois.c ./qbinom.c, ./qnbinom{,_mu}.c
 */
 
 #define PST_0(a, b) a ## b
@@ -31,7 +31,7 @@
 #define _qDIST_  PASTE(q, _thisDIST_)
 /**
    For a discrete distribution on  the integers,
-   For P(x) := <pDist>(x, <distPars>),  find p-quantile  y(p)   :<==>   P(y) < p <= P(y)
+   For P(x) := <pDist>(x, <distPars>),  find p-quantile  y = y(p)   :<==>   P(y-1) < p <= P(y)
 
    @param y    current guess
    @param *z   := <pDist>(y, ..)
@@ -90,38 +90,44 @@ static double DO_SEARCH_FUN(_dist_PARS_DECL_)
     }
     else { // (lower_tail, *z < p)  or  (upper tail, *z >= p): search to the __right__
 	for(int iter = 0; ; iter++) {
+	    double prevy = y;
+	    double newz = -1.; // -Wall
 #ifndef MATHLIB_STANDALONE
 	    if(iter % 10000 == 0) R_CheckUserInterrupt();
 #endif
 	    y += incr;
 #ifdef _dist_MAX_y
 	    if(y < _dist_MAX_y)
-		*z = P_DIST(y, _dist_PARS_);
+		newz = P_DIST(y, _dist_PARS_);
 	    else if(y > _dist_MAX_y)
 		y = _dist_MAX_y;
 #else
-	    *z = P_DIST(y, _dist_PARS_);
+	    newz = P_DIST(y, _dist_PARS_);
 #endif
 
 	    if(
 #ifdef _dist_MAX_y
 		y == _dist_MAX_y ||
 #endif
-		ISNAN(*z) || (lower_tail ? (*z >= p) : (*z < p)))
+		ISNAN(newz) || (lower_tail ? (newz >= p) : (newz < p)))
 	    {
 		R_DBG_printf("  new y=%.15g, z=%g = " AS_CHAR(_pDIST_) "(y,*) %s;"
-			     " ==> search() returns after %d iter.\n", y, *z,
-			     ISNAN(*z) ? "is NaN" : (lower_tail ? ">= p" : "< p"), iter);
-		return y;
+			     " ==> search() returns after %d iter.\n", y, newz,
+			     ISNAN(newz) ? "is NaN" : (lower_tail ? ">= p" : "< p"), iter);
+		if (incr <= 1) {
+		    *z = newz;
+		    return y;
+		}
+		return prevy;
 	    }
+	    *z = newz;
 	}
     }
 } // do_search()
 
 
 /*
-* Note : "same" code in qbinom.c, qnbinom.c __FIXME__ NOT YET for qpois() ??
-* FIXME: This is far from optimal [cancellation for p ~= 1, etc]:
+* Note : called in qbinom.c, qnbinom.c but not (yet) qpois.c -- NB: only DBG_print()ing; *no other* effect
 */
 #define q_DISCRETE_01_CHECKS() do {					\
     double p_n; /* p  in the "normal" (lower_tail=TRUE, log.p=FALSE) scale */ \
diff --git a/src/qpois.c b/src/qpois.c
index e6f87d3..b5c81e9 100644
--- a/src/qpois.c
+++ b/src/qpois.c
@@ -64,7 +64,7 @@ double qpois(double p, double lambda, int lower_tail, int log_p)
     double
 	mu = lambda,
 	sigma = sqrt(lambda),
-	// had gamma = sigma; PR#8058 should be kurtosis which is mu^-0.5 = 1/sigma
+	// had gamma = sigma; PR#8058 should be skewness which is mu^-0.5 = 1/sigma
 	gamma = 1.0/sigma;
 
      R_DBG_printf("qpois(p=%.12g, lambda=%.15g, l.t.=%d, log=%d):"

From f9c3f1c91941f5254a39c85392f9a8147935fe0f Mon Sep 17 00:00:00 2001
From: Andreas Noack <andreas@noack.dk>
Date: Tue, 10 Sep 2024 20:16:23 +0200
Subject: [PATCH 2/2] Add README comment about _Thread_local

---
 README.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/README.md b/README.md
index e45c5c6..bdcfa3e 100644
--- a/README.md
+++ b/README.md
@@ -22,4 +22,4 @@ Updating
 To update to the latest version of R, bump the `RVERSION` file, and run `make
 update`. Some additional manual changes to the headers may be necessary: these should go
 in `include/Rconfig.h` (this would typically be generated by autotools, but we try to
-simplify the build process).
+simplify the build process). Please also check that declarations with `_Thread_local` are not dropped as part of the update, see https://github.com/JuliaStats/Rmath-julia/pull/50.