Rmath-julia

This is a slightly modified version of the standalone Rmath library from R, built to be used with the Rmath.jl Julia package.

The main difference is that it is built to allow defining custom random number generating functions via C function pointers (see include/callback.h). When using the library, these should be defined before calling any of the random functions.

Build instructions

Rmath-julia requires GNU Make (https://www.gnu.org/software/make). Just run make to compile the library.

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

make update re-extracts the upstream nmath sources over src/ and then re-applies patches/thread-local.patch, so please check that the generators' per-thread state survives the update — see the section below and #50.

The patch must contain only files the extraction actually overwrites — currently the seven in src/ that R ships under src/nmath/. Files we keep locally are not re-extracted, so a hunk for one of them would be re-applied to an already-patched file: patch reports Reversed (or previously applied) patch detected, and non-interactively it assumes -R and undoes the change. src/sunif.c is the one to watch, since R ships it under src/nmath/standalone/ which the extraction excludes; its _Thread_local seed is committed directly and must stay out of the patch. To regenerate, diff a pristine extraction against src/ and keep only the overwritten files.

Per-thread state

Rmath keeps the random generators' setup constants in function-local static variables: a memoisation cache so that repeated calls with unchanged parameters skip the setup arithmetic. That is not thread-safe, so this fork makes the state per-thread — but it must stay off the static TLS segment. On glibc, a dlopen'd module's TLS block is drawn from a fixed ~1664-byte per-process surplus shared with every other shared library, and it is only tunable (glibc.rtld.optional_static_tls) from glibc 2.32. Declaring the state _Thread_local cost 856 bytes, over half the surplus, which made unrelated libraries fail to load with cannot allocate memory in static TLS block — see #56.

So the state lives on the heap, in one Rmath_tls container reached through a single thread-local pointer (src/rmath_tls.h). The budget is now 16 bytes: that pointer, plus sunif.c's seed. CI asserts the segment stays small, so a regression here fails the build rather than surfacing later as somebody else's load failure.

Each generator reaches its own slice of the container at the point where it first needs it, and then aliases the fields back to upstream's names:

    struct rbeta_state *st = &Rmath_tls_get()->rbeta;
#define beta	st->beta
#define olda	st->olda

That keeps the function bodies byte-identical to R's, so an R release that edits the algorithm does not conflict with this fork: the only upstream lines the patch removes are the static declarations themselves. If a future R release introduces a local variable with the same name as one of the aliased fields, the alias expands into it and the compile fails at that line — noisy, but not silent. Fix it by renaming the alias, never the upstream local.

Three rules when touching this code:

1. New mutable per-thread state goes in Rmath_tls, not in a new _Thread_local. Add a field to the relevant struct <name>_state, and if it needs a non-zero initial value, set it in that file's Rmath_<name>_state_init() hook — the hook is the mechanical image of the initialisers Rmath writes on the declaration, so static double olda = -1.0; becomes st->olda = -1.0;. Keeping the hook in the same file as the code that reads the value means an update that changes a sentinel shows both halves in one patch hunk. Rmath_tls_alloc() callocs, so fields that are zero upstream need no hook entry — but write them out anyway to keep the mapping one-for-one.

2. Fetch the state after the early returns, not on entry. Every generator has escape paths — argument validation, degenerate parameters, an INT_MAX fallback that delegates elsewhere — that touch none of the cached state. Fetching below them means a caller that only ever hits those paths never allocates the container at all.

3. Only per-thread state belongs there at all. Read-only coefficient tables are already thread-safe; making them thread-local advertises per-thread state that does not exist, and at -O3 the compiler folds them anyway. Leave them exactly as upstream has them — which is const static for fact[] in src/rpois.c but plain static for the coefficients in rexpm1() in src/toms708.c; both are fine, and neither is worth a diff. A quick way to check whether a variable is genuinely written: add const to its declaration and see whether the compiler rejects an assignment (use -ferror-limit=0, or clang stops at 20 diagnostics and you will miss some).

Cleanup is automatic — src/rmath_tls.c registers a thread-exit destructor via pthread_key_create() on POSIX and FlsAlloc() on Windows — so there is nothing for callers to free and no public API for it. Windows deliberately avoids pthreads: on mingw-w64 that resolves through winpthreads and adds a runtime dependency on libwinpthread-1.dll, which Rmath_jll does not ship.