On Mon, Mar 24, 2025 at 08:46:50PM +0200, Janne Blomqvist wrote: > On Mon, Mar 24, 2025 at 8:25 PM Steve Kargl > <s...@troutmask.apl.washington.edu> wrote: > > > > I've already written a prototype of the half-cycle trig > > functions. > > > > https://gcc.gnu.org/bugzilla/show_bug.cgi?id=113152 > > > > There are two issues that need to be address. First, some > > operating systems provide half-cycle trig functions in their > > libm. The initial patch tries to use libm functions if > > found, but in hindsight I think gfortran should simply > > provide its own implementations. ” > > Hello Steve, long time no see. AFAIU the motivation for these > half-cycle functions is to avoid the non-trivial gymnastics associated > with argument reduction of a multiple of pi. Both in terms of cpu > time, and in the possibility of getting it wrong and not producing a > correctly rounded result. It seems to me that if libm provides said > functions it would be preferably to call those directly rather than > always using some fallback implementation that ends up calling the > "normal" trig functions. Not sure how that should be implemented, > maybe some kind of weak symbol trickery to use the libgfortran > fallback implementations in case libm doesn't provide it? >
Hi Janne, Indeed, it's been awhile. I wrote FreeBSD's sinpi(), cospi(), and tanpi(). I haven't contributed asinpi() and friends to FreeBSD, yet. The argument reduction isn't too tricky. One simply breaks 'x' into 'x = n + r' with 'n' some integer and '0 <= r < 1'. The algorithm description for sinpi() is in the leading comment here https://cgit.freebsd.org/src/tree/lib/msun/src/s_sinpi.c Similar comments are in cospi() and tanpi(). The nice features of these functions are the special cases, e.g., 'sinpi(int(x)) = +-0' exactly. The sign is chosen based on whether int(x) is even or odd. I think the biggest issues will come with REAL(10) and REAL(16) on X86 hardware. The former is typically C's 'long double' and the functions are likely present in libm. REAL(16) is __Float128, which I have not used, so would/will need to learn how to use. For non-X86 hardware, REAL(16) is either IEEE binary128 type or IBM double-double representation. If the OS's libm does not have an appropriate function, we'll need a fallback (and I would likely write it in Fortran). One other issue to consider is that weak symbol trickery does not work with static linking. -- Steve
