mdebug

This is a tool for debugging over-sensitivity issues of Windows programs (primarily R packages) to the accuracy of C Math library functions on Windows. Mingw-w64 v12 switched a number of C Math functions from an internal implementation in mingw-w64 to the C runtime, so to UCRT. This made a number of R packages to fail their tests. When debugging these issues, it is important to find out which of the Math functions is (are) causing the issue and how far off the correctly rounded (correct) result it is. If it were very inaccurate, one would report a bug against mingw-w64. If it were just a very small difference (say 1 ULP), clearly the package test/code is wrong (over-sensitive) and needs to be adjusted. For more background, see blog posts Sensitivity to C math library and mingw-w64 v12 and Sensitivity to C math library and mingw-w64 v12 - part 2.

In principle, a C Math function can be used from R itself (so the package may be invoking it from R code), or from a dependent package, or from the failing package itself, or from an external library used by some of the involved packages.

This tools provides a DLL with Math function wrappers. One then builds R from source, and all involved R packages, linking them against this DLL. Wrappers are created for all Math functions that in mingw-w64 v12 switched from an internal implementation to UCRT. The wrappers are controlled by environment variables. One can simply log calls to these functions, to know which are used by a specific test. And then one can selectively choose which functions should be used from UCRT (so, finding the function causing the failure by bisecting). The wrapper DLL should be built with mingw-w64 v11 (so that the internal implementations are available, the UCRT functions (if used), are used dynamically).

Building the wrapper DLL

Using Rtools45 (so still mingw-w64 v11), one would run

dlltool --input-def mdebug.def --output-lib libm.a
gcc -o mdebug.dll -shared -I. mdebug.c wrappers.c mdebug.def

In the above, libm.a is a static library, so called "interface" or "import" library on Windows. The library is created from the definition file mdebug.f, which starts with

LIBRARY mdebug.dll
EXPORTS 
  acos = acos_debug
  acosh = acosh_debug
  asin = asin_debug
...

so for example the import library will provide function acos, which will call into acos_debug in mdebug.dll. Then, mdebug.dll is a shared library with implementations of functions like acos_debug. In practice, one would build R and R packages, linking to libm.a. One would then arrange for mdebug.dll to be discovered by the OS. One can modify the implementation of the wrappers, rebuild mdebug.dll and replace it, without rebuilding R and packages.

Building R

To build R (and then packages) against the wrapper import library libm.a, one needs to modify the make file in R sources accordingly. An example patch is provided in mdebug.patch (it assumes libm.a is in c:/mdebug):

cd trunk
patch -p0 \<mdebug.patch

The building of R includes also running R itself (e.g. when building packages). That already needs mdebug.dll. So, one way to proceed is to start the build without it. The build will fail, but create (among others) target directory bin/x64. One can then place a copy of mdebug.dll there and re-start the build.

Building R packages

R packages can be built using R built above - the patch includes also the necessary modification of the linking of packages. One needs to re-build all packages that use native code from source, but ideally even those not using native code (not very likely, but their build could still in theory capture some Math library results).

Testing the build

Try this example in R:

x1 <- -log(-expm1(-0.5))
x2 <- -log1p(-exp(-x1))/0.5
x2 == 1
sprintf("%a", x2)

x2 == 1 was true with mingw-w64 v11, but false in v12 (on at least a number of Windows 10 machines and Windows 2022). It should be noted that with v12, as UCRT is used, the results are OS-version specific.

So, run the example as above in the newly built R, and expect a result like

> x2 == 1
[1] TRUE
> sprintf("%a", x2)
[1] "0x1p+0"

Now, run with Math functions from UCRT (as v12 would do; the command assumes running from Msys2 shell, if running otherwise, set the environment variable accordingly):

env MDEBUG_UCRT=all R

and expect a result like (this is what I get on Windows 10 build 19045):

> x2 == 1
[1] FALSE
> sprintf("%a", x2)
[1] "0x1.0000000000001p+0"

Logging calls to Math functions

Calls to Math functions can be logged (traced) to a file. The path is 'c:/temp/mdebug.log' by default and can be changed by environment variable MDEBUG_LOGFILE. Lets run the example above (in file t.r):

env MDEBUG_LOGFILE=t.log MDEBUG_LOG=all Rscript.exe t.r

and expect a result like

$ head t.log
log10 ( 0.000000(0x1.0000000000000p-1022) ) -> -307.652656(-0x1.33a7146f72a42p+8)
floor ( -307.652656(-0x1.33a7146f72a42p+8) ) -> -308.000000(-0x1.3400000000000p+8)
floor ( 0.000000(0x0.0000000000000p+0) ) -> 0.000000(0x0.0000000000000p+0)
floor ( 3.125000(0x1.9000000000001p+1) ) -> 3.000000(0x1.8000000000000p+1)
floor ( 0.000000(0x0.0000000000000p+0) ) -> 0.000000(0x0.0000000000000p+0)
floor ( 3.125000(0x1.9000000000001p+1) ) -> 3.000000(0x1.8000000000000p+1)
floor ( 0.000000(0x0.0000000000000p+0) ) -> 0.000000(0x0.0000000000000p+0)
floor ( 3.125000(0x1.9000000000001p+1) ) -> 3.000000(0x1.8000000000000p+1)
floor ( 0.000000(0x0.0000000000000p+0) ) -> 0.000000(0x0.0000000000000p+0)
floor ( 3.125000(0x1.9000000000001p+1) ) -> 3.000000(0x1.8000000000000p+1)

so, the log file includes names of the functions called, including the arguments and results. To get the count of functions called

$ cat t.log | cut -d' ' -f1 | sort | uniq -c | sort -n
      2 exp
      2 expm1
      2 log1p
      4 log
     16 log10
    270 floor

Bisecting Math function causing a failure

Lets find which of the Math functions causes the difference in the test output. Lets use only half of functions from UCRT:

env MDEBUG_UCRT="exp expm1 log1p" Rscript.exe t.r

the output is FALSE, so one of these is causing the problem. So lets try each of those, and we eventually find that expm1 alone causes the problem (and none of the other ones):

env MDEBUG_UCRT="log1p" Rscript.exe t.r
[1] TRUE
[1] "0x1p+0"

Measuring the inaccuracy

Lets check how much different are results of expm1 in UCRT compared to v11.

env MDEBUG_LOGFILE=d.log MDEBUG_LOG=all MDEBUG_UCRT="expm1" Rscript.exe t.r

The log file with this simple example only has one line:

-1 ulp expm1 ( -0.500000(-0x1.0000000000000p-1) ) -> -0.393469(-0x1.92e9a0720d3edp-2)

where -1 is the difference between the results using v11 and UCRT in ULPs (units of least precision).

As a check, lets run the same example with all functions from v11 and compare the log outputs:

env MDEBUG_LOGFILE=dd.log  MDEBUG_LOG=expm1 ../../bin/Rscript.exe t.r
expm1 ( -0.500000(-0x1.0000000000000p-1) ) -> -0.393469(-0x1.92e9a0720d3ecp-2)

So we can se that the result of expm1(-0.5) changed from -0x1.92e9a0720d3ecp-2 to -0x1.92e9a0720d3edp-2), so there is a change in the last hex digit from c to d, hence in the very last bit, hence 1 ULP

Limitations

The tool doesn't have any synchronization for writing to the log file. Hence, tracing code with concurrency may cause some lines in the file to be garbled - but, in practice, it has never been a big enough problem so far, it was still possible to use the tool to narrow down the Math function causing a change in test outputs (in principle, indeed, these can be multiple functions, it may be that one is not enough).

In practice, I've sometimes modified the code of selected wrappers, to narrow down issues in more complicated packages. For example, when one runs into a situation when there are many calls to the function causing the difference, and some of the calls provide very different results (many ULPs), but - are those calls really the ones causing the package to fail its checks? To know, I would have modified the sources to say use UCRT result only when the difference is 1 ULP. So, one might have to do such things when using this tool.