Re: [Fortran] half-cycle trig functions and atan[d] fixes
On Wed, Jan 24, 2024 at 11:13:05AM +0200, Janne Blomqvist wrote: > On Wed, Jan 24, 2024 at 10:28 AM FX Coudert wrote: > > > Now, if > > > the OS adds cospi() to libm and it's in libm's symbol map, then the > > > cospi() used by gfortran depends on the search order of the loaded > > > libraries. > > > > We only include the fallback math functions in libgfortran when they are > > not available on the system. configure detects what is present in the libc > > being targeted, and conditionally compiles the necessary fallback functions > > (and only them). > > Exactly. However, there is the (corner?) case when libgfortran has > been compiled, and cospi() not found and thus the fallback > implementation is included, and then later libc is updated to a > version that does provide cospi(). I believe in that case which > version gets used is down to the library search order (i.e. the order > that "ldd /path/to/binary" prints the libs), it will use the first > symbol it finds. Also, it's not necessary to do some ifdef tricks > with gfortran.map, if a symbol listed there isn't found in the library > it's just ignored. So the *pi() trig functions can be unconditionally > added there, and then depending on whether the target libm includes > those or not they are then included in the exported symbol list. > > It's possible to override this to look for specific symbol versions > etc., but that probably goes deep into the weeds of target-specific > stuff (e.g. are we looking for cospi@FBSD_1.7, cospi@GLIBC_X.Y.Z, or > something else?). I'm sure you don't wanna go there. > Ah, so that's the part I was missing. I was under the impression that if a symbol appears in a libraries symbol map, then the library had to contain a function by that name. If the loader ignores symbols for a missing function, then yes, I think I can get rid of the indirection via _gfortran_cospi_r4(). It will take a few days for me to redesign this, which shouldn't be too problematic in that GCC is in stage 4 and this is neither a regression or doc fix. Janne, FX, Harald, thanks for taking a peek. -- Steve
Re: [Fortran] half-cycle trig functions and atan[d] fixes
Am 24.01.24 um 10:13 schrieb Janne Blomqvist: On Wed, Jan 24, 2024 at 10:28 AM FX Coudert wrote: Now, if the OS adds cospi() to libm and it's in libm's symbol map, then the cospi() used by gfortran depends on the search order of the loaded libraries. We only include the fallback math functions in libgfortran when they are not available on the system. configure detects what is present in the libc being targeted, and conditionally compiles the necessary fallback functions (and only them). Exactly. However, there is the (corner?) case when libgfortran has been compiled, and cospi() not found and thus the fallback implementation is included, and then later libc is updated to a version that does provide cospi(). I believe in that case which version gets used is down to the library search order (i.e. the order that "ldd /path/to/binary" prints the libs), it will use the first symbol it finds. Also, it's not necessary to do some ifdef tricks with gfortran.map, if a symbol listed there isn't found in the library it's just ignored. So the *pi() trig functions can be unconditionally added there, and then depending on whether the target libm includes those or not they are then included in the exported symbol list. It's possible to override this to look for specific symbol versions etc., but that probably goes deep into the weeds of target-specific stuff (e.g. are we looking for cospi@FBSD_1.7, cospi@GLIBC_X.Y.Z, or something else?). I'm sure you don't wanna go there. Isn't this something that could be addressed by __attribute__(__weakref__)?
Re: [Fortran] half-cycle trig functions and atan[d] fixes
On Wed, Jan 24, 2024 at 10:28 AM FX Coudert wrote: > > Now, if > > the OS adds cospi() to libm and it's in libm's symbol map, then the > > cospi() used by gfortran depends on the search order of the loaded > > libraries. > > We only include the fallback math functions in libgfortran when they are not > available on the system. configure detects what is present in the libc being > targeted, and conditionally compiles the necessary fallback functions (and > only them). Exactly. However, there is the (corner?) case when libgfortran has been compiled, and cospi() not found and thus the fallback implementation is included, and then later libc is updated to a version that does provide cospi(). I believe in that case which version gets used is down to the library search order (i.e. the order that "ldd /path/to/binary" prints the libs), it will use the first symbol it finds. Also, it's not necessary to do some ifdef tricks with gfortran.map, if a symbol listed there isn't found in the library it's just ignored. So the *pi() trig functions can be unconditionally added there, and then depending on whether the target libm includes those or not they are then included in the exported symbol list. It's possible to override this to look for specific symbol versions etc., but that probably goes deep into the weeds of target-specific stuff (e.g. are we looking for cospi@FBSD_1.7, cospi@GLIBC_X.Y.Z, or something else?). I'm sure you don't wanna go there. -- Janne Blomqvist
Re: [Fortran] half-cycle trig functions and atan[d] fixes
Hi, > Hopefully, FX sees this as my emails to gmail bounce. I am seeing this email. > Now, if > the OS adds cospi() to libm and it's in libm's symbol map, then the > cospi() used by gfortran depends on the search order of the loaded > libraries. We only include the fallback math functions in libgfortran when they are not available on the system. configure detects what is present in the libc being targeted, and conditionally compiles the necessary fallback functions (and only them). FX
Re: [Fortran] half-cycle trig functions and atan[d] fixes
On Tue, Jan 23, 2024 at 01:37:54PM +0200, Janne Blomqvist wrote: > > > - If I get this right, to take one example, the Fortran front-end will emit > > a call to gfortran_acospi_r4(), libgfortran provides this as a wrapper > > calling acospif(), which is called either from libm or from libgfortran. > > This is different from other math library functions, like ACOS() where the > > acosf() call is generated directly from the front-end (and then the > > implementation comes either from libm or from libgfortran). Why not follow > > our usual way of doing things? > > Good point, that's what we've done in c99_functions.c in libgfortran. > We should probably do this so we can skip jumping via libgfortran when > libm implements these directly. > Hopefully, FX sees this as my emails to gmail bounce. I don't see how this can work or I don't understand symbol versioning in libraries. If an OS does not supply cospi() in its libm, then gfortran will fallback to a cospi() in libgfortran. With the indirection I currently implement, cospi() is not in the symbol map of libgfortran and _gfortran_cospi_r4() will use the libgfortran version. Now, if the OS adds cospi() to libm and it's in libm's symbol map, then the cospi() used by gfortran depends on the search order of the loaded libraries. Consider on FreeBSD, I have % nm --dynamic /lib/libm.so.5 | grep cospi 00025b60 T cospi@@FBSD_1.7 00025fe0 T cospif@@FBSD_1.7 0002e230 T cospil@@FBSD_1.7 % nm --dynamic work/lib/libgfortran.so.5 | grep cospi 002b1e60 T _gfortran_acospi_r10@@GFORTRAN_14 002b4590 T _gfortran_acospi_r16@@GFORTRAN_14 002b1d80 T _gfortran_acospi_r4@@GFORTRAN_14 002b1df0 T _gfortran_acospi_r8@@GFORTRAN_14 002b1ea0 T _gfortran_cospi_r10@@GFORTRAN_14 002b46d0 T _gfortran_cospi_r16@@GFORTRAN_14 002b1dc0 T _gfortran_cospi_r4@@GFORTRAN_14 002b1e30 T _gfortran_cospi_r8@@GFORTRAN_14 U cospi@FBSD_1.7 U cospif@FBSD_1.7 U cospil@FBSD_1.7 The FE generates code for _gfortran_cospi_rXX. If FreeBSD adds say acospif() and libgfortran also exported, then how does the linker choose between acospif@@FBSD_1.7 and acospif@@GFORTRAN_14?. -- Steve
Re: [Fortran] half-cycle trig functions and atan[d] fixes
On Tue, Jan 23, 2024 at 01:37:54PM +0200, Janne Blomqvist wrote: > On Tue, Jan 23, 2024 at 11:09 AM FX Coudert wrote: > > > > Hi Steve, > > Hello, long time no see. Time is short and we're all busy with life, but it is nice to see familiar names! > > > Thanks for the patch. I’ll take time to do a proper review, but after a > > first read I had the following questions: > > > > - "an OS's libm may/will contain cospi(), etc.”: do those functions conform > > to any standard? Are there plans to implement them outside FreeBSD at this > > point? > > acospi, asinpi, atan2pi, atanpi, cospi, sinpi, tanpi are all in IEEE > 754-2019, and are slated to be part of C23. I would assume actively > developed libm's will eventually catch up and implement them. > > > - If I get this right, to take one example, the Fortran front-end will emit > > a call to gfortran_acospi_r4(), libgfortran provides this as a wrapper > > calling acospif(), which is called either from libm or from libgfortran. > > This is different from other math library functions, like ACOS() where the > > acosf() call is generated directly from the front-end (and then the > > implementation comes either from libm or from libgfortran). Why not follow > > our usual way of doing things? > > Good point, that's what we've done in c99_functions.c in libgfortran. > We should probably do this so we can skip jumping via libgfortran when > libm implements these directly. > > > - On most targets, cospi() and friends are not available. Therefore, we end > > up doing the fallback (with limited precision as you noted) but with a lot > > of indirection. We could generate that code directly in the front-end, > > couldn’t we? > > The frontend generally doesn't know what the target libm implements, > hence it's better to just generate the call, and if necessary have a > barebones implementation in libgfortran if libm doesn't implement it > properly. > I suppose I could add a gfc_conv_trigpi() to trans-intrinsic.c, which could be used emit the function calls directly. Give me until Saturday to look at the issue. I've finally understood what Harald has been trying to tell me about the MPFR version issue in bugzilla, so I need to revamp simplification. -- Steve
Re: [Fortran] half-cycle trig functions and atan[d] fixes
On Tue, Jan 23, 2024 at 06:40:27AM -0800, Steve Kargl wrote:
> On Tue, Jan 23, 2024 at 10:08:43AM +0100, FX Coudert wrote:
> > Hi Steve,
> >
> > Thanks for the patch. I’ll take time to do a proper review, but
> > after a first read I had the following questions:
> >
> > - "an OS's libm may/will contain cospi(), etc.”: do those functions
> > conform to any standard? Are there plans to implement them outside
> >FreeBSD at this point?
>
> Yes. cospi, sinpi, and tanpi are in IEEE754-2008. These are
> also in ISO/IEC TS 18661-4 along with acospi, asinpi, and atanpi,
> which I believe where merged into C23 (see n3096.pdf). I've
I haven't checked.
> checked if atan2pi is in C23, but it is in F2023.
>
> AFAIK, FreeBSD's libm is the only OS that contains cospi, sinpi,
> and tanpi. The arc functions haven't been written/committed,
> because something like cospi(x) = cos(x) / pi with pi split
> into hi-lo parts is good enough.
Whoops. acospi(x) = acos(x) / pi. The fallback implements this
as acospi(x) = acos(x) * invpihi + acos(x) * invpilo with invpihi
the leading digits(x)/2 bits of 1/pi and invpilo the trailing
digits(x) of 1/pi.
For most libm's, acos(x) with |x| <= 1 have good numerical accuracy
and handle subnormal, |x| > 1, +-inf, and NaN by setting appropriate
exceptions. With FreeBSD and exhaustive testing in the interval, I see
% tlibm acos -f -x 0x1p-9 -X 1 -PED
Interval tested for acosf: [0.00195312,1]
ulp <= 0.5: 96.432% 72803871 | 96.432% 72803871
0.5 < ulp < 0.6: 3.110% 2348017 | 99.542% 75151888
0.6 < ulp < 0.7: 0.419%316134 | 99.961% 75468022
0.7 < ulp < 0.8: 0.039% 29450 | 100.000% 75497472
Max ulp: 0.796035 at 4.99814630e-01
Exhaustive testing of the fallback routine for acospi(x) gives
program foo
implicit none
integer n
real(4) f4, x, xm
real(8) f8, u, v
u = -1
x = scale(1., -9)
do
f4 = acospi(x)
f8 = acospi(real(x,8))
n = exponent(f8)
v = abs(f8 - f4)
v = scale(v, digits(f4) - n)
if (v > u) then
u = v
xm = x
end if
x = nearest(x, 2.)
if (x > 1) exit
end do
print *, u, acospi(xm), acospi(real(xm,8))
end program foo
% gfcx -o z -O a.f90 && ./z
1.9551682807505131 0.337791681 0.33779173955822828
so a max ulp of 1.955. Hopefully, OS libm's will catch up and
provide an implementation that gets the extra 1+ bit of precision.
Admittedly, the fallback routines for cospi(), sinpi(), and
tanpi() could be better, but much slower. For now, I do for
example,
float
cospif (float x)
{
return cosf (x * pihi_f + x * pilo_f);
}
A better routine would be,
/* cospi(x) = cos(pi*x) = cos(pi*n+pi*r) = +-cos(pi*r) with 0 <= r < 1i
and n integer. */
float
cospif (float x)
{
int s;
float ax, n, r;
ax = fabsf(ax);
if (ax > 0x1p23) {
return 1;
} else {
/* FreeBSD uses bit twiddling. See
https://cgit.freebsd.org/src/tree/lib/msun/src/s_cospif.c */
n = floor(ax); /* integer part */
r = ax - n; /* remainder */
s = floor(n/2) * 2 - n == 0 : 1 : -1; /* even n? */
return s * cosf (r * pi);
}
}
> > - On systems where libquadmath is used for _Float128, does
> > libquadmath contain the implementation of the q() variants for
> > those functions?
>
> AFAIK, libquadmath does not have the half-cycle functions. I
> wrote the function trig_fallback2.F90 to deal with REAL(16) (and
> maybe REAL17).
>
> > - If I get this right, to take one example, the Fortran front-end
> > will emit a call to gfortran_acospi_r4(), libgfortran provides this
> > as a wrapper calling acospif(), which is called either from libm
> > or from libgfortran.
>
> Yes, that is correct. I tried to install __builtin_acospif in
> trans-intrinsic to generate a direct call to acospif, but that
> led to many hours/days of frustration. I gave up.
>
> > This is different from other math library functions, like ACOS()
> > where the acosf() call is generated directly from the front-end
> > (and then the implementation comes either from libm or from
> > libgfortran). Why not follow our usual way of doing things?
>
> I gave up on that approach when I got some real obscure error
> about some math function which I did not touch.
>
I tried adding
DEFINE_MATH_BUILTIN_C (COSPI, "cospi", 0)
or
DEFINE_MATH_BUILTIN (COSPI, "cospi", 0)
or
OTHER_BUILTIN (COSPI, "cospi", 1, false) (and variations).
to mathbuiltins.def but this led to some obscure error message
about some other unrelated intrinsic subprogram.
--
Steve
Re: [Fortran] half-cycle trig functions and atan[d] fixes
On Tue, Jan 23, 2024 at 10:08:43AM +0100, FX Coudert wrote: > Hi Steve, > > Thanks for the patch. I’ll take time to do a proper review, but > after a first read I had the following questions: > > - "an OS's libm may/will contain cospi(), etc.”: do those functions > conform to any standard? Are there plans to implement them outside >FreeBSD at this point? Yes. cospi, sinpi, and tanpi are in IEEE754-2008. These are also in ISO/IEC TS 18661-4 along with acospi, asinpi, and atanpi, which I believe where merged into C23 (see n3096.pdf). I've checked if atan2pi is in C23, but it is in F2023. AFAIK, FreeBSD's libm is the only OS that contains cospi, sinpi, and tanpi. The arc functions haven't been written/committed, because something like cospi(x) = cos(x) / pi with pi split into hi-lo parts is good enough. > - On systems where libquadmath is used for _Float128, does > libquadmath contain the implementation of the q() variants for > those functions? AFAIK, libquadmath does not have the half-cycle functions. I wrote the function trig_fallback2.F90 to deal with REAL(16) (and maybe REAL17). > - If I get this right, to take one example, the Fortran front-end > will emit a call to gfortran_acospi_r4(), libgfortran provides this > as a wrapper calling acospif(), which is called either from libm > or from libgfortran. Yes, that is correct. I tried to install __builtin_acospif in trans-intrinsic to generate a direct call to acospif, but that led to many hours/days of frustration. I gave up. > This is different from other math library functions, like ACOS() > where the acosf() call is generated directly from the front-end > (and then the implementation comes either from libm or from > libgfortran). Why not follow our usual way of doing things? I gave up on that approach when I got some real obscure error about some math function which I did not touch. > - On most targets, cospi() and friends are not available. Therefore, > we end up doing the fallback (with limited precision as you noted) > but with a lot of indirection. We could generate that code directly > in the front-end, couldn’t we? The precision of the fallback routines isn't too bad. More later today. Late for a doctor's approintment. -- Steve
Re: [Fortran] half-cycle trig functions and atan[d] fixes
On Tue, Jan 23, 2024 at 11:09 AM FX Coudert wrote: > > Hi Steve, Hello, long time no see. > Thanks for the patch. I’ll take time to do a proper review, but after a first > read I had the following questions: > > - "an OS's libm may/will contain cospi(), etc.”: do those functions conform > to any standard? Are there plans to implement them outside FreeBSD at this > point? acospi, asinpi, atan2pi, atanpi, cospi, sinpi, tanpi are all in IEEE 754-2019, and are slated to be part of C23. I would assume actively developed libm's will eventually catch up and implement them. > - If I get this right, to take one example, the Fortran front-end will emit a > call to gfortran_acospi_r4(), libgfortran provides this as a wrapper calling > acospif(), which is called either from libm or from libgfortran. This is > different from other math library functions, like ACOS() where the acosf() > call is generated directly from the front-end (and then the implementation > comes either from libm or from libgfortran). Why not follow our usual way of > doing things? Good point, that's what we've done in c99_functions.c in libgfortran. We should probably do this so we can skip jumping via libgfortran when libm implements these directly. > - On most targets, cospi() and friends are not available. Therefore, we end > up doing the fallback (with limited precision as you noted) but with a lot of > indirection. We could generate that code directly in the front-end, couldn’t > we? The frontend generally doesn't know what the target libm implements, hence it's better to just generate the call, and if necessary have a barebones implementation in libgfortran if libm doesn't implement it properly. -- Janne Blomqvist
Re: [Fortran] half-cycle trig functions and atan[d] fixes
Hi Steve, Thanks for the patch. I’ll take time to do a proper review, but after a first read I had the following questions: - "an OS's libm may/will contain cospi(), etc.”: do those functions conform to any standard? Are there plans to implement them outside FreeBSD at this point? - On systems where libquadmath is used for _Float128, does libquadmath contain the implementation of the q() variants for those functions? - If I get this right, to take one example, the Fortran front-end will emit a call to gfortran_acospi_r4(), libgfortran provides this as a wrapper calling acospif(), which is called either from libm or from libgfortran. This is different from other math library functions, like ACOS() where the acosf() call is generated directly from the front-end (and then the implementation comes either from libm or from libgfortran). Why not follow our usual way of doing things? - On most targets, cospi() and friends are not available. Therefore, we end up doing the fallback (with limited precision as you noted) but with a lot of indirection. We could generate that code directly in the front-end, couldn’t we? Best, FX
[Fortran] half-cycle trig functions and atan[d] fixes
All, I have attached a new patch to https://gcc.gnu.org/bugzilla/show_bug.cgi?id=113152 which addresses the following issues. PR113152 -- implement half-cycle trigonometric functions PR113412 -- better error message for atan(y,x) PR113413 -- implement atand(y,x) The patch clocks in at 3488 lines, so will not be attached here. There are a few items that may need some explanation. First, some (if not all) of the half-cycle functions are contained in ISO/IEC TS 18661-4, IEEE754-2008, and now Fortran 2023. The patch assumes that an OS's libm may/will contain cospi(), etc. (see for example FreeBSD's libm). If these are not in libm, then fallback routines are provided by libgfortran. The fallback routines have seen limited testing, so your mileage may vary. Second, I have no idea what REAL17 entails; however, I provide untested routines in trigpi_fallback2.F90 where I use REAL(16) in the declarations (yikes). Third, MPFR provides half-cycle functions in version 4.2.0 (and later). For simplicifation, the patch checks for an adequate MPFR version, if an older version is found fallback algorithms are used. These have had limited testing. The patch appears to pass regression testing. Hopefully, someone will find time to test and commit the patch. The remainder is an attempt to write ChangeLog entries. gcc/fortran: * gfortran.h (gfc_isym_id): Add GFC_ISYM_ACOSPI, GFC_ISYM_ASINPI, GFC_ISYM_ATANPI, GFC_ISYM_ATAN2PI, GFC_ISYM_COSPI, GFC_ISYM_SINPI, and GFC_ISYM_TANPI * intrinsic.cc (do_check): Whitespace and typo in comments. (add_functions): Add two-argument form of ATAND. Add half-cycle trigonometric functions ACOSPI, ASINPI, ATANPI, ATAN2PI, COSPI, SINPI, and TANPI. (sort_actual): Generate sensible error messages for two argument versions of ATAN, ATAND, and ATANPI. * intrinsic.h: New prototypes for gfc_simplify_acospi, gfc_simplify_asinpi, gfc_simplify_atanpi, gfc_simplify_atan2pi, gfc_simplify_cospi, gfc_simplify_sinpi, gfc_simplify_tanpi, gfc_resolve_acospi, gfc_resolve_asinpi, gfc_resolve_atanpi, gfc_resolve_atan2pi, gfc_resolve_cospi, gfc_resolve_sinpi, and gfc_resolve_tanpi * intrinsic.texi: Document new functions ACOSPI, ASINPI, ATANPI, ATAN2PI, COSPI, SINPI, and TANPI. Put the ATAN* family of functions in lexigraphical order. * iresolve.cc (gfc_resolve_acospi, gfc_resolve_asinpi, gfc_resolve_atanpi, gfc_resolve_atan2pi, gfc_resolve_cospi, gfc_resolve_sinpi, gfc_resolve_tanpi): New functions. * simplify.cc (gfc_simplify_acospi, gfc_simplify_asinpi, gfc_simplify_atanpi, gfc_simplify_atan2pi, gfc_simplify_cospi, gfc_simplify_sinpi, gfc_simplify_tanpi): New functions. Introduce MPFR_HALF_CYCLE macros to use MPFR half-cycle functions if available. * trans-intrinsic.cc: Declare asinpi, acospi, atanpi, atan2pi, sinpi, cospi, and tanpi as LIB_FUNCTION's. libgfortran: *Makefile.am: New files trigpi.c, trigpi_fallback1.c, and trigpi_fallback2.F90 * configure.ac: Check libm for float, double, and long double functions for asinpi, acospi, atanpi, atan2pi, sinpi, cospi, and tanpi * Makefile.in: Regenerated. * config.h.in: Ditto. * configure: Ditto. * gfortran.map: Add GFORTRAN_14 section to expose the library symbols. * trigpi.c: Fallback implementations of the half-cycle trigonometric functions. * trigpi_fallback1.c: Fallback functions for float, double, and long double routines if libm does provide them. * trigpi_fallback2.F90: REAL(16)/REAL(17?) fallback implementations of the half-cycle trigonometric functions. -- Steve
