[Fortran] half-cycle trig functions and atan[d] fixes

[Fortran] half-cycle trig functions and atan[d] fixes

[Steve Kargl]

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

Re: [Fortran] half-cycle trig functions and atan[d] fixes

[FX Coudert]

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

Re: [Fortran] half-cycle trig functions and atan[d] fixes

[Janne Blomqvist]

On Tue, Jan 23, 2024 at 11:09 AM FX Coudert <fxcoudert@gmail.com> 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

[Steve Kargl]

On Tue, Jan 23, 2024 at 01:37:54PM +0200, Janne Blomqvist wrote:
> On Tue, Jan 23, 2024 at 11:09 AM FX Coudert <fxcoudert@gmail.com> 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

[Steve Kargl]

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
0000000000025b60 T cospi@@FBSD_1.7
0000000000025fe0 T cospif@@FBSD_1.7
000000000002e230 T cospil@@FBSD_1.7

% nm --dynamic work/lib/libgfortran.so.5 | grep cospi
00000000002b1e60 T _gfortran_acospi_r10@@GFORTRAN_14
00000000002b4590 T _gfortran_acospi_r16@@GFORTRAN_14
00000000002b1d80 T _gfortran_acospi_r4@@GFORTRAN_14
00000000002b1df0 T _gfortran_acospi_r8@@GFORTRAN_14
00000000002b1ea0 T _gfortran_cospi_r10@@GFORTRAN_14
00000000002b46d0 T _gfortran_cospi_r16@@GFORTRAN_14
00000000002b1dc0 T _gfortran_cospi_r4@@GFORTRAN_14
00000000002b1e30 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

[FX Coudert]

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

[Janne Blomqvist]

On Wed, Jan 24, 2024 at 10:28 AM FX Coudert <fxcoudert@gmail.com> 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

[Harald Anlauf]

Am 24.01.24 um 10:13 schrieb Janne Blomqvist:
> On Wed, Jan 24, 2024 at 10:28 AM FX Coudert <fxcoudert@gmail.com> 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

[Steve Kargl]

On Wed, Jan 24, 2024 at 11:13:05AM +0200, Janne Blomqvist wrote:
> On Wed, Jan 24, 2024 at 10:28 AM FX Coudert <fxcoudert@gmail.com> 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

[Steve Kargl]

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

[Steve Kargl]

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