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

[Steve Kargl <sgk@troutmask.apl.washington.edu>]

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