On Thu, Jul 03, 2025 at 02:43:43PM +0200, Michael Matz wrote:
> Hello,
>
> On Thu, 3 Jul 2025, Yuao Ma wrote:
>
> > This patch adds the required function for Fortran trigonometric functions to
> > work with glibc versions prior to 2.26. It's based on glibc source commit
> > 632d895f3e5d98162f77b9c3c1da4ec19968b671.
> >
> > I've built it successfully on my end. Documentation is also included.
> >
> > Please take a look when you have a moment.
>
> +__float128
> +cospiq (__float128 x)
> +{
> ...
> + if (__builtin_islessequal (x, 0.25Q))
> + return cosq (M_PIq * x);
>
> Isn't the whole raison d'etre for the trig-pi functions that the internal
> argument reduction against multiples of pi becomes trivial and hence (a)
> performant, and (b) doesn't introduce rounding artifacts? Expressing the
> trig-pi functions in terms of their counterparts completely defeats this
> purpose. The other way around would be more sensible for the cases where
> it works, but the above doesn't seem very attractive.
glibc has
FLOAT
M_DECL_FUNC (__cospi) (FLOAT x)
{
if (isless (M_FABS (x), M_EPSILON))
return M_LIT (1.0);
if (__glibc_unlikely (isinf (x)))
__set_errno (EDOM);
x = M_FABS (x - M_LIT (2.0) * M_SUF (round) (M_LIT (0.5) * x));
if (islessequal (x, M_LIT (0.25)))
return M_SUF (__cos) (M_MLIT (M_PI) * x);
else if (x == M_LIT (0.5))
return M_LIT (0.0);
else if (islessequal (x, M_LIT (0.75)))
return M_SUF (__sin) (M_MLIT (M_PI) * (M_LIT (0.5) - x));
else
return -M_SUF (__cos) (M_MLIT (M_PI) * (M_LIT (1.0) - x));
}
for this case, so if it is incorrect, has too bad precision
or there are better ways to do it, then perhaps it should be changed
on the glibc side first and then ported to libquadmath.
Jakub
