Note that this works fine in the WASM version:
?-A is nexttoward(1, 10).✎✕
A = 1.0000000000000002.
That doesn’t help as @ridgeworks says. I think the only way out is to use a C library that implements software floats with correct rounding. That has been on the table with @ridgeworks . I forgot the name of the library. If I recall well though, the portability was non-trivial (which will hurt here) and the library is not maintained AFAIK. The other possibility is to use LibBF, which is already used in the WASM version for big ints and rationals because it is so much smaller than libGMP The library also implements arbitrary precision floats and, if I recall correctly, rounding modes. It implements the usual +, -, * and / and a couple of functions:
int bf_const_log2(bf_t *T, limb_t prec, bf_flags_t flags);
int bf_const_pi(bf_t *T, limb_t prec, bf_flags_t flags);
int bf_exp(bf_t *r, const bf_t *a, limb_t prec, bf_flags_t flags);
int bf_log(bf_t *r, const bf_t *a, limb_t prec, bf_flags_t flags);
#define BF_POW_JS_QUIRKS (1 << 16) /* (+/-1)^(+/-Inf) = NaN, 1^NaN = NaN */
int bf_pow(bf_t *r, const bf_t *x, const bf_t *y, limb_t prec, bf_flags_t flags)
;
int bf_cos(bf_t *r, const bf_t *a, limb_t prec, bf_flags_t flags);
int bf_sin(bf_t *r, const bf_t *a, limb_t prec, bf_flags_t flags);
int bf_tan(bf_t *r, const bf_t *a, limb_t prec, bf_flags_t flags);
int bf_atan(bf_t *r, const bf_t *a, limb_t prec, bf_flags_t flags);
int bf_atan2(bf_t *r, const bf_t *y, const bf_t *x,
limb_t prec, bf_flags_t flags);
int bf_asin(bf_t *r, const bf_t *a, limb_t prec, bf_flags_t flags);
int bf_acos(bf_t *r, const bf_t *a, limb_t prec, bf_flags_t flags);
That could be used to implement double operations for a subset of the SWI-Prolog float functions with rounding. A little macro layer can than be used to select these rather than the C library versions. I.e, just introduce #define PL_SIN(x) sin(x), etc.
I think that route keep the code fairly elegant and is quite feasible.