Merge pull request #45 from JuliaStats/vs/r-4.4.0

Update to R 4.4.0

.github/workflows/CI.yml

@@ -17,16 +17,20 @@ jobs:
         version:
           - '1.0'
           - '1'
-          - 'nightly'
         os:
           - ubuntu-latest
+          - windows-latest
         arch:
           - x64
+        include:
+          - os: macOS-14
+            arch: aarch64
+            version: '1'
     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

RVERSION

@@ -1 +1 @@
-4.2.2
+4.4.0

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);
 
 

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

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

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

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)

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))

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:

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).