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 correctly-rounded elementary functions.
(reference: https://hal-ens-lyon.archives-ouvertes.fr/ensl-01529804/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 IEEE-754 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 IEEE-754 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 Xh-Xl.
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 Xh-Xl.
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.6613381477509392425e-16'
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 pl-arith.c and additional test cases as required.
Alternatively we leave the test failing and possibly add a bug note to the round function. Finally we can provide our own correct lround for systems where it is broken.