So hereās a somewhat extreme solution but I think it covers a lot of the bases. ācrlibm
ā is:
a portable, proven, correctly rounded, and efficient mathematical library (libm) for double precision. ā¦ The ultimate goal of the crlibm project is to push towards the standardization of correctlyrounded elementary functions.
(reference: https://halenslyon.archivesouvertes.fr/ensl01529804/file/crlibm.pdf contains way more than you want to know but first chapter or two is useful background).
Key attributes:
 correctly rounded  i.e., āas perfect as can beā
 portable 
crlibm
is written in C and will be compiled by any compiler fulfilling basic requirements of the ISO/IEC 9899:1999 (C99) standard.
 proven 
crlibm
intends to provide a comprehensive proof of the theoretical possibility of correct rounding, the algorithms used, and the implementation, assuming C99 and IEEE754 compliance.
 efficient  performance and resource usage of
crlibm
should be comparable to existing libm
implementations, both in average and in the worst case.
A corollary is that given C99 and IEEE754 compliance, the results are the same on all platforms; a useful property for composing test cases. And the need for floating point ādisclaimersā is minimized, if required at all.
crlibm
is an INRIA sponsored project and has been around for over 15 years. I found additional language bindings for Python (among the top 5% of packages on PyPi), Julia, and Rust with a quick search. The main argument against using it as a general replacement for libm is that it only supports doubles, but thatās all thatās required for our purposes.
crlibm
is not a replacement library for libm
; rather it provides a set of four functions (1 for each rounding direction) for each double elementary function (e.g., exp, sin, etc.)  see https://github.com/taschini/crlibm/blob/eb3063791aa75bc9705b49283bf14250465220a7/crlibm.h. So one attractive possibility is to make crlibm
an optional library in SWIP, analogous to the GMP library for rationals. To take advantage of its capabilities would require changes to the various āar_
ā elementary functions in pl_arith.c
, but this could largely be done with a macro to replace the existing calls to ālibmā. Iām not a C programmer but I assume it would look something like:
#define CALL_ELEM_FUNCTION1(func,n,r) \
#ifdef CRLIBM \
int rounding = fegetround(); \
if (rounding == FE_TONEAREST) \
r>value.f = func_rn(n>value.f); \
else {
fesetround(FE_TONEAREST); \
switch(rounding) { \
FE_UPWARD : r>value.f = func_ru(n>value.f); \
FE_DOWNWARD : r>value.f = func_rd(n>value.f); \
FE_TOWARDZERO : r>value.f = func_rz(n>value.f); \
} \
fesetround(rounding); \
} \
#else \
r>value.f = func(n>value.f); \
#endif
If the flag CRLIBM
is defined, this just calls the version of the function corresponding to the current rounding mode after setting the rounding mode to FE_TONEAREST
(expected by crlibm
) and afterwards restoring the mode. If the flag is not defined, the corresponding funtion in libm
is used. I think this would handle most of the elementary function calls, but there may be a few outliers, e.g., ar_pow
. In addition, crlibm
does not provide direct implementations of atan2
, erf
, erfc
, and lgamma
, so these need further investigation. (Worst case, these just revert to libm
with a disclaimer.)
The big advantage of this is to isolate SWIP from the differences (including bugs such as the one that initiated this thread) in the libm
elementary functions between the various platforms, which should improve reliability and simplify testing. In addition, the results should be directly comparable with other language systems using crlibm
. Most users are unlikely to notice any changes, but the affected function results are guaranteed precise and correctly follow the IEEE rounding mode semantics, contrary to my experience with existing math libraries. Iāve done some superficial testing using the Python implementation and crlibm
seems to fulfill its promise.
One example using SWISH (8.5.9):
? Xl is roundtoward(exp(log(2)),to_negative), Xh is roundtoward(exp(log(2)),to_negative), W is XhXl.
W = 0.0,
Xh = Xl, Xl = 2.0
On macOS (8.4.1):
? Xl is roundtoward(exp(log(2)),to_negative), Xh is roundtoward(exp(log(2)),to_negative), W is XhXl.
Xl = Xh, Xh = 1.9999999999999998,
W = 0.0.
Note that rounding up and rounding down the same expression yields an identical answer, which should never be the case for correct rounding. With Python crlibm
:
>>> format(crlibm.exp_rd(crlibm.log_rd(2)),'.20g')
'1.999999999999999778'
>>> format(crlibm.exp_ru(crlibm.log_ru(2)),'.20g')
'2.0000000000000004441'
>>> format(crlibm.exp_ru(crlibm.log_ru(2))crlibm.exp_rd(crlibm.log_rd(2)),'.20g')
'6.6613381477509392425e16'
Outstanding issues:

I assume that given the source here https://github.com/taschini/crlibm, supporting crlibm
in SWIP in much the same way as is currently done for GMP is not a problem, i.e. ājustā incremental work; is this a valid assumption?

Licensing: Like GMP, crlibm
has an LGPL license so I assume it would positioned the same way, i.e., as a build option (as reflected in the draft CALL_ELEM_FUNCTION1 macro above).

Environment flag: any need? (Itās unlikely the average user will see any difference.)

Who?: I donāt feel comfortable tinkering with the SWIP build process but I can help with PRās for the detailed changes to plarith.c
and additional test cases as required.