Sorry for the long post but here’s a Plan B which, while not quite as good as Plan A, may be feasible in the short term. To recap:
As most know, unlike arithmetic using rational numbers (which include integers), arithmetic using floating point values does not obey the standard mathematical laws (associativity, distributivity, etc.) due to rounding errors. But arithmetic over reals can be made sound using sets of real values defined by floating point bounds; various schemes include interval arithmetic, ball arithmetic, and enum’s, but they’re all based on the same principle. However this requires control over rounding as defined by the IEEE floating point standard. This control is only mandated to be correct over the normal arithmetic operations and sqrt
and this means there are no such guarantees over other real functions like exp
, log
, trig. functions, etc. and this can be easily demonstrated using a simple predicate like:
round_value(Exp,[Rn,Rd,Ru,Rz]) :-
Rn is roundtoward(Exp,to_nearest),
Rd is roundtoward(Exp,to_negative),
Ru is roundtoward(Exp,to_positive),
Rz is roundtoward(Exp,to_zero).
A few examples using SWISH:
?- X is pi/4, round_value(sin(X),[Rn,Rd,Ru,Rz]).
Rd = Rn, Rn = Ru, Ru = Rz, Rz = 0.7071067811865475,
X = 0.7853981633974483
?- X is -pi/4, round_value(sin(X),[Rn,Rd,Ru,Rz]).
Rd = Rn, Rn = Ru, Ru = Rz, Rz = -0.7071067811865475,
X = -0.7853981633974483
?- X is pi/2, round_value(sin(X),[Rn,Rd,Ru,Rz]).
Rd = Rn, Rn = Ru, Ru = Rz, Rz = 1.0,
X = 1.5707963267948966
?- X is -pi/2, round_value(sin(X),[Rn,Rd,Ru,Rz]).
Rd = Rn, Rn = Ru, Ru = Rz, Rz = -1.0,
X = -1.5707963267948966
?- X is pi,round_value(X**X,[Rn,Rd,Ru,Rz]).
Rd = Rn, Rn = Ru, Ru = Rz, Rz = 36.4621596072079,
X = 3.141592653589793
So the math library used by SWIP on SWISH just seems to ignore the rounding mode. At least it’s consistent, although it’s not correct because, except for exact values (which these are not), rounding up and rounding down should never produce the same value.
The other platform I have access to is MacOS which generates quite different results:
?- X is pi/4, round_value(sin(X),[Rn,Rd,Ru,Rz]).
X = 0.7853981633974483,
Rn = Rd, Rd = Rz, Rz = 0.7071067811865475,
Ru = 0.7071067811865476.
?- X is -pi/4, round_value(sin(X),[Rn,Rd,Ru,Rz]).
X = -0.7853981633974483,
Rn = Rd, Rd = Ru, Ru = -0.7071067811865475,
Rz = -0.7071067811865476.
?- X is pi/2, round_value(sin(X),[Rn,Rd,Ru,Rz]).
X = 1.5707963267948966,
Rn = 1.0,
Rd = Rz, Rz = 0.9999999999999999,
Ru = 1.000000000251173.
?- X is -pi/2, round_value(sin(X),[Rn,Rd,Ru,Rz]).
X = -1.5707963267948966,
Rn = Ru, Ru = -1.0,
Rd = -1.000000000251173,
Rz = -0.9999999999999999.
?- X is pi, round_value(X**X,[Rn,Rd,Ru,Rz]).
X = 3.141592653589793,
Rn = Rd, Rd = Rz, Rz = 36.4621596072079,
Ru = 36.46215960720791.
The result for pi/4
is actually correct, but -pi/4
is definitely wrong since rounding towards zero actually produces a more negative value than any of the other modes. And, as in the incorrect SWISH results, rounding up and rounding down produces the same value. The tests around ±pi/2
demonstrates a significant bug on the order of the 6 least significant digits and produces an answer outside the ±1.0 range of the sin
function so definitely broken there.
Using crlibm
on MacOS produces correct results:
?- X is pi/4, round_value(sin(X),[Rn,Rd,Ru,Rz]).
X = 0.7853981633974483,
Rn = Rd, Rd = Rz, Rz = 0.7071067811865475,
Ru = 0.7071067811865476.
?- X is -pi/4, round_value(sin(X),[Rn,Rd,Ru,Rz]).
X = -0.7853981633974483,
Rn = Ru, Ru = Rz, Rz = -0.7071067811865475,
Rd = -0.7071067811865476.
?- X is pi/2, round_value(sin(X),[Rn,Rd,Ru,Rz]).
X = 1.5707963267948966,
Rn = Ru, Ru = 1.0,
Rd = Rz, Rz = 0.9999999999999999.
?- X is -pi/2, round_value(sin(X),[Rn,Rd,Ru,Rz]).
X = -1.5707963267948966,
Rn = Rd, Rd = -1.0,
Ru = Rz, Rz = -0.9999999999999999.
?- X is pi, round_value(X**X,[Rn,Rd,Ru,Rz]).
X = 3.141592653589793,
Rn = 36.4621596072079,
Rd = Rz, Rz = 36.462159607207894,
Ru = 36.46215960720791.
So here’s three different implementations of the floating point elementary function, all with different answers. And other than crlibm
the answers are wrong depending on the rounding mode. But the thing to note is they all agree for the default “to_nearest
” rounding mode which is almost universally used except by those of us interested in correctness.
So if we can make the (perhaps dubious) assumption that the various math libraries get the correct “to_nearest
” answer for these functions within 1ulp (unit of least precision), rounding can be implemented using the nexttoward
C function. Such rounding is not as precise as possible (normally 2ulp’s between round up and round down values instead of 1) and some care must be taken to ensure rounding doesn’t violate the mathematical definition of the function (e.g., sin
values > 1.0) or to round infinite values to finite ones, but that’s pretty manageable. It only requires an addition to pl-arith.c
using the redefinition of library functions as for the latest crlibm
version (no new libraries to import or build config issues) and produces the following results:
?- X is pi/4, round_value(sin(X),[Rn,Rd,Ru,Rz]).
X = 0.7853981633974483,
Rn = 0.7071067811865475,
Rd = Rz, Rz = 0.7071067811865474,
Ru = 0.7071067811865476.
?- X is -pi/4, round_value(sin(X),[Rn,Rd,Ru,Rz]).
X = -0.7853981633974483,
Rn = -0.7071067811865475,
Rd = -0.7071067811865476,
Ru = Rz, Rz = -0.7071067811865474.
?- X is pi/2, round_value(sin(X),[Rn,Rd,Ru,Rz]).
X = 1.5707963267948966,
Rn = Ru, Ru = 1.0,
Rd = Rz, Rz = 0.9999999999999999.
?- X is -pi/2, round_value(sin(X),[Rn,Rd,Ru,Rz]).
X = -1.5707963267948966,
Rn = Rd, Rd = -1.0,
Ru = Rz, Rz = -0.9999999999999999.
?- X is pi,round_value(X**X,[Rn,Rd,Ru,Rz]).
X = 3.141592653589793,
Rn = 36.4621596072079,
Rd = Rz, Rz = 36.462159607207894,
Ru = 36.46215960720791.
The good:
- produces correctly rounded results with a fairly straight forward addition to
pl-arith.c
in one place; no other code changes, and no conditional compilation, required
- no changes to build process (i.e., the
cmake
“quagmire” - my term)
- no software licensing issues
The bad:
- depends on the possibly dubious assumption that the math libraries generate answers within 1ulp in (default)
to_nearest
mode. (Some simple testing similar to the above examples would help with this determination.)
- produce slightly larger range of rounding values (2ulp’s instead of 1)
The ugly:
- nothing really; the vast majority of users don’t go near rounding modes and the executed code in those cases remains the practically the same other than a quick check to see what rounding is in effect when these elementary functions are called.
Summary - not as good a solution as Plan A using crlibm
but probably good enough, and much better than the status quo.