lib/58237: modfl returns wrong answers on ld128 architectures

[campbell+netbsd%mumble.net@localhost]

>Number:         58237
>Category:       lib
>Synopsis:       modfl returns wrong answers on ld128 architectures
>Confidential:   no
>Severity:       serious
>Priority:       medium
>Responsible:    lib-bug-people
>State:          open
>Class:          sw-bug
>Submitter-Id:   net
>Arrival-Date:   Wed May 08 22:50:00 +0000 2024
>Originator:     Taylor R Campbell
>Release:        current, 10, 9, ...
>Organization:
The NetBSD Modflation
>Environment:
>Description:
f = modfl(x, &i) is supposed to set f to be a fractional part lying in [0,1) and i to be an integer.

But on ld128 architectures such as aarch64, it fails to do this in some cases.  For example, for 1 + FLT_EPSILON (that is, 1 + 2^{-23}), it returns f = 0 and i = 1 + 2^{-23}.
>How-To-Repeat:
long double i, f;

f = modfl(1 + FLT_EPSILON, &i);

For example:

https://releng.netbsd.org/b5reports/evbarm-aarch64/2024/2024.05.08.02.08.11/test.html#lib_libm_t_modf_modfl
>Fix:
Yes, please!