math consistency: fallback code used for sub_/add_epsilon on systems without 'have_nextafter' weak for subnormal values?
nextafter( x, go_ninf )
and ~go_pinf
are used in goffice/goffice/math/go-math.c/go_sub_epsilon
and ~add_epsilon
to produce a 'grace' rounding for INT, CEIL, FLOOR, CEILING, TRUNC, ROUNDDOWN, ROUNDUP, ROUND and evtl. others. For systems lacking HAVE_NEXTAFTER
there is fallback code with frexp and ldexp.
I tested this code by de-activating ( comenting out ) the part
#ifdef HAVE_NEXTAFTER
return x == 0 ? x : nextafter (x, go_ninf);
#else
and saw it weak - do nothing - for subnormal values. That produces different results to HAVE_NEXTAFTER systems where nextafter( x, go_pinf ) changes the origin.
Despite I can't tell whether grace / fuzzy rounding is meaningful in this range because the granularity is increasing dramatically for very small values, I can say that we should try to have the code consistent for systems with or without 'HAVE_NEXTAFTER'. Either de-activate sub/add_epsilon for 'normal' code, or implement it for subnormals too.
Both could be easy to code, either exclude the range, or add / subtract 5E-324 to / from subnormals instead of other action ( 4E-4951 for long double version ).
When reworking this it could evtl. be meaningful / faster to calculate the nextafters by multiplying / dividing the origin with 0.9999999999999999 ( 0.99999999999999999995 for long doubles ). Tested that for some 64k random values, nextafter
, sub/add_epsilon
and *0.99999999999999999
/ /0.9999999999999999
resp. the corresponding long double value do match for normal values.