Re: PATCH -- Fix degree trignometric functions

classic Classic list List threaded Threaded
21 messages Options
12
Reply | Threaded
Open this post in threaded view
|

Re: PATCH -- Fix degree trignometric functions

Steve Kargl
On Fri, Mar 06, 2020 at 03:18:19PM -0800, Steve Kargl wrote:
>
>   3. Simplification routines do the following mappings:
>      sind(x) = sin((pi/180) * x)     asind(x) = (180/pi) * asin(x)
>      cosd(x) = cos((pi/180) * x)     acosd(x) = (180/pi) * acos(x)
>      tand(x) = tan((pi/180) * x)     atand(x) = (180/pi) * atan(x)
>      atan2d(y,x) = (180/pi) * atan2(y,x)
>      cotand(x) = cotan((pi/180) * x)
>      All computations are carried out by MPFR or MPC.

In looking at some basic tests, the above simplification will
need to be modified to do what ...


>   5. New functions have been added to libgfortran to handle sind, cosd,
>      and tand.

... these functions do.  Otherwise, things like cos(real(60+123*360))
give wrong values.  Modifications are so easy even a lurker can do
them.

--
steve
Reply | Threaded
Open this post in threaded view
|

Re: PATCH -- Fix degree trignometric functions

Fritz Reese
On Sat, Mar 7, 2020 at 1:32 PM Steve Kargl
<[hidden email]> wrote:
>
> Fix the simplification and handling of the degree trigonometric functions.
>  This includes fixing a number of ICEs.  See PR 93871.
>
>  ChangeLog and patch attached.

As the author of the original degree-trig functions I intend to review
this patch soon, and...

> On Fri, Mar 06, 2020 at 03:18:19PM -0800, Steve Kargl wrote:
> >
> >   3. Simplification routines do the following mappings:
> >      sind(x) = sin((pi/180) * x)     asind(x) = (180/pi) * asin(x)
> >      cosd(x) = cos((pi/180) * x)     acosd(x) = (180/pi) * acos(x)
> >      tand(x) = tan((pi/180) * x)     atand(x) = (180/pi) * atan(x)
> >      atan2d(y,x) = (180/pi) * atan2(y,x)
> >      cotand(x) = cotan((pi/180) * x)
> >      All computations are carried out by MPFR or MPC.
>
> In looking at some basic tests, the above simplification will
> need to be modified to do what ...
>
>
> >   5. New functions have been added to libgfortran to handle sind, cosd,
> >      and tand.
>
> ... these functions do.  Otherwise, things like cos(real(60+123*360))
> give wrong values.  Modifications are so easy even a lurker can do
> them.

... extend the patch to include these changes (unless someone
enthusiastic gets around to these mods before I do).

I should be able to start on this next week (around 16 March).

---
Fritz Reese
Reply | Threaded
Open this post in threaded view
|

Re: PATCH -- Fix degree trignometric functions

gcc - fortran mailing list
Steve,

> > On Fri, Mar 06, 2020 at 03:18:19PM -0800, Steve Kargl wrote:
[...]
> > >   5. New functions have been added to libgfortran to handle sind, cosd,
> > >      and tand.

Much of the patch looks good. Thanks for your work on this.

However, after applying your patch, the new libgfortran functions
appear to have a calling-convention mismatch between the frontend and
the library, at least for my platform.

$ nl sind_intrinsic.f90
     1    implicit none

     2    real(8) :: xout
     3    real(8) :: xin = 33.5_8

     4    xout = dsind(xin)
     5    print *, xout

     6    end
$ gfortran -g -O0 sind_intrinsic.f90
$ gfortran -dumpmachine
x86_64-pc-linux-gnu
$ ./a.out

Program received signal SIGSEGV: Segmentation fault - invalid memory reference.

Backtrace for this error:
#0  0x7fd442ece3af in ???
#1  0x40076d in MAIN__
    at /data/foreese/src/gcc-git/gcc/testsuite/gfortran.dg/sind_intrinsic.f90:6
#2  0x40080f in main
    at /data/foreese/src/gcc-git/gcc/testsuite/gfortran.dg/sind_intrinsic.f90:9
Segmentation fault (core dumped)

$ gdb -ex "set breakpoint pending on" -ex "b _gfortran_sind_r8" ./a.out
[...]
(gdb) run
[...]
Breakpoint 1, _gfortran_sind_r8 (x=0x19) at
/data/foreese/src/gcc-git/libgfortran/intrinsics/trigd_inc.c:44
44      if (isfinite (*x))
(gdb) p x
$2 = (GFC_REAL_8 *) 0x19
(gdb) disas $rip,+16
Dump of assembler code from 0x7ffff7b91a50 to 0x7ffff7b91a60:
=> 0x00007ffff7b91a50 <_gfortran_sind_r8+0>:    movsd  (%rdi),%xmm1
   0x00007ffff7b91a54 <_gfortran_sind_r8+4>:    movsd
0x14bcc(%rip),%xmm2        # 0x7ffff7ba6628
   0x00007ffff7b91a5c <_gfortran_sind_r8+12>:    pxor   %xmm3,%xmm3
End of assembler dump.
(gdb) p/x $rdi
0x19
(gdb) up
#1  0x000000000040076e in MAIN__ () at sind_intrinsic.f90:6
6    xout = dsind(xin)
(gdb) disas $rip-24,+24
Dump of assembler code from 0x400756 to 0x40076e:
   0x0000000000400756 <MAIN__+4>:    sub    $0x220,%rsp
   0x000000000040075d <MAIN__+11>:    mov    0x200904(%rip),%rax
 # 0x601068 <xin.3898>
   0x0000000000400764 <MAIN__+18>:    movq   %rax,%xmm0
   0x0000000000400769 <MAIN__+23>:    callq  0x400640 <_gfortran_sind_r8@plt>
End of assembler dump.


Seems like for some reason the calling code passes the floating-point
parameter directly by-value in the SSE register %xmm0, while the
library routine expects a pointer to the value in %rdi. This is the IR
with  -fdump-tree-original:

$ head -n6 sind_intrinsic.f90.004t.original
MAIN__ ()
{
  static real(kind=8) xin = 3.35e+1;
  real(kind=8) xout;

  xout = _gfortran_sind_r8 (xin);
$

An interesting note is that packaging the call into a function appears
to work at first...

$ nl sind_function.f90
     1    module mod_sind
     2      contains
     3      function lib_sind(x)
     4        implicit none
     5        real(8) :: lib_sind
     6        real(8), intent(in) :: x
     7        lib_sind = SIND(x)
     8      end function
     9    end module

    10    use mod_sind
    11    implicit none

    12    real(8) :: xout
    13    real(8) :: xin = 33.5_8

    14    xout = lib_sind(xin)
    15    print *, xout

    16    end
$ gfortran -g -O0 sind_function.f90 -o sind_function -fdump-tree-original
$ ./sind_function
  0.55193698531205815

The IR for comparison:

$ head -n15 sind_function.f90.004t.original
lib_sind (real(kind=8) & restrict x)
{
  real(kind=8) __result_lib_sind;

  __result_lib_sind = _gfortran_sind_r8 (*x);
  return __result_lib_sind;
}


MAIN__ ()
{
  static real(kind=8) xin = 3.35e+1;
  real(kind=8) xout;

  xout = lib_sind (&xin);
$

... BUT by looking at the disassembly you can see the same error is
actually present, and it is an accident that the first parameter to
the wrapper function lib_sind is accidentally leftover in %rdi when
entering _gfortran_sind_r8.

$ objdump -d sind_function | grep -B5 'call.*gfortran_sind'
  400756:    48 83 ec 20              sub    $0x20,%rsp
  40075a:    48 89 7d e8              mov    %rdi,-0x18(%rbp)
  40075e:    48 8b 45 e8              mov    -0x18(%rbp),%rax
  400762:    48 8b 00                 mov    (%rax),%rax
  400765:    66 48 0f 6e c0           movq   %rax,%xmm0
  40076a:    e8 d1 fe ff ff           callq  400640 <_gfortran_sind_r8@plt>

You can confirm this by changing lib_sind to accept two parameters and
pass the second to sind(). The first value ends up in
_gfortran_sind_r8. (Of course the same issue occurs for COSD and TAND,
both with REAL(4) and REAL(8).)


Unfortunately, I'm not that familiar with middle-end translation yet,
so it's not immediately clear to me why the calling code is passing
by-value instead of by-reference to the library function. I don't see
significant differences between the way these functions are wrapped
compared to other libgfortran routines which do pass by-reference. The
biggest difference I see is the LIB_FUNCTION() lines in the
gfc_intrinsic_map in trans-intrinsic.c: why are these present for the
trigd functions and not for other libgfortran functions?

$ grep -C4 'LIB_FUNCTION (SIND' trans-intrinsic.c | nl -ba -v 119
   119    #include "mathbuiltins.def"
   120
   121      /* Functions in libgfortran.  */
   122      LIB_FUNCTION (ERFC_SCALED, "erfc_scaled", false),
   123      LIB_FUNCTION (SIND, "sind", false),
   124      LIB_FUNCTION (COSD, "cosd", false),
   125      LIB_FUNCTION (TAND, "tand", false),
   126
   127      /* End the list.  */


Steve, do you have any ideas? I will continue to poke around.

---
Fritz Reese
Reply | Threaded
Open this post in threaded view
|

Re: PATCH -- Fix degree trignometric functions

gcc - fortran mailing list
On Mon, Mar 16, 2020 at 08:27:25PM -0400, Fritz Reese wrote:

> Steve,
>
> > > On Fri, Mar 06, 2020 at 03:18:19PM -0800, Steve Kargl wrote:
> [...]
> > > >   5. New functions have been added to libgfortran to handle sind, cosd,
> > > >      and tand.
>
> Much of the patch looks good. Thanks for your work on this.
>
> However, after applying your patch, the new libgfortran functions
> appear to have a calling-convention mismatch between the frontend and
> the library, at least for my platform.
>
> $ nl sind_intrinsic.f90
>      1    implicit none
>
>      2    real(8) :: xout
>      3    real(8) :: xin = 33.5_8
>
>      4    xout = dsind(xin)
>      5    print *, xout
>
>      6    end
> $ gfortran -g -O0 sind_intrinsic.f90
> $ gfortran -dumpmachine
> x86_64-pc-linux-gnu
> $ ./a.out
>
> Program received signal SIGSEGV: Segmentation fault - invalid memory reference.
>
> Backtrace for this error:
> #0  0x7fd442ece3af in ???
> #1  0x40076d in MAIN__
>     at /data/foreese/src/gcc-git/gcc/testsuite/gfortran.dg/sind_intrinsic.f90:6
> #2  0x40080f in main
>     at /data/foreese/src/gcc-git/gcc/testsuite/gfortran.dg/sind_intrinsic.f90:9
> Segmentation fault (core dumped)
>
> $ gdb -ex "set breakpoint pending on" -ex "b _gfortran_sind_r8" ./a.out
> [...]
> (gdb) run
> [...]
> Breakpoint 1, _gfortran_sind_r8 (x=0x19) at
> /data/foreese/src/gcc-git/libgfortran/intrinsics/trigd_inc.c:44
> 44      if (isfinite (*x))
> (gdb) p x
> $2 = (GFC_REAL_8 *) 0x19
> (gdb) disas $rip,+16
> Dump of assembler code from 0x7ffff7b91a50 to 0x7ffff7b91a60:
> => 0x00007ffff7b91a50 <_gfortran_sind_r8+0>:    movsd  (%rdi),%xmm1
>    0x00007ffff7b91a54 <_gfortran_sind_r8+4>:    movsd
> 0x14bcc(%rip),%xmm2        # 0x7ffff7ba6628
>    0x00007ffff7b91a5c <_gfortran_sind_r8+12>:    pxor   %xmm3,%xmm3
> End of assembler dump.
> (gdb) p/x $rdi
> 0x19
> (gdb) up
> #1  0x000000000040076e in MAIN__ () at sind_intrinsic.f90:6
> 6    xout = dsind(xin)
> (gdb) disas $rip-24,+24
> Dump of assembler code from 0x400756 to 0x40076e:
>    0x0000000000400756 <MAIN__+4>:    sub    $0x220,%rsp
>    0x000000000040075d <MAIN__+11>:    mov    0x200904(%rip),%rax
>  # 0x601068 <xin.3898>
>    0x0000000000400764 <MAIN__+18>:    movq   %rax,%xmm0
>    0x0000000000400769 <MAIN__+23>:    callq  0x400640 <_gfortran_sind_r8@plt>
> End of assembler dump.
>
>
> Seems like for some reason the calling code passes the floating-point
> parameter directly by-value in the SSE register %xmm0, while the
> library routine expects a pointer to the value in %rdi. This is the IR
> with  -fdump-tree-original:
>
> $ head -n6 sind_intrinsic.f90.004t.original
> MAIN__ ()
> {
>   static real(kind=8) xin = 3.35e+1;
>   real(kind=8) xout;
>
>   xout = _gfortran_sind_r8 (xin);
> $
>
> An interesting note is that packaging the call into a function appears
> to work at first...
>
> $ nl sind_function.f90
>      1    module mod_sind
>      2      contains
>      3      function lib_sind(x)
>      4        implicit none
>      5        real(8) :: lib_sind
>      6        real(8), intent(in) :: x
>      7        lib_sind = SIND(x)
>      8      end function
>      9    end module
>
>     10    use mod_sind
>     11    implicit none
>
>     12    real(8) :: xout
>     13    real(8) :: xin = 33.5_8
>
>     14    xout = lib_sind(xin)
>     15    print *, xout
>
>     16    end
> $ gfortran -g -O0 sind_function.f90 -o sind_function -fdump-tree-original
> $ ./sind_function
>   0.55193698531205815
>
> The IR for comparison:
>
> $ head -n15 sind_function.f90.004t.original
> lib_sind (real(kind=8) & restrict x)
> {
>   real(kind=8) __result_lib_sind;
>
>   __result_lib_sind = _gfortran_sind_r8 (*x);
>   return __result_lib_sind;
> }
>
>
> MAIN__ ()
> {
>   static real(kind=8) xin = 3.35e+1;
>   real(kind=8) xout;
>
>   xout = lib_sind (&xin);
> $
>
> ... BUT by looking at the disassembly you can see the same error is
> actually present, and it is an accident that the first parameter to
> the wrapper function lib_sind is accidentally leftover in %rdi when
> entering _gfortran_sind_r8.
>
> $ objdump -d sind_function | grep -B5 'call.*gfortran_sind'
>   400756:    48 83 ec 20              sub    $0x20,%rsp
>   40075a:    48 89 7d e8              mov    %rdi,-0x18(%rbp)
>   40075e:    48 8b 45 e8              mov    -0x18(%rbp),%rax
>   400762:    48 8b 00                 mov    (%rax),%rax
>   400765:    66 48 0f 6e c0           movq   %rax,%xmm0
>   40076a:    e8 d1 fe ff ff           callq  400640 <_gfortran_sind_r8@plt>
>
> You can confirm this by changing lib_sind to accept two parameters and
> pass the second to sind(). The first value ends up in
> _gfortran_sind_r8. (Of course the same issue occurs for COSD and TAND,
> both with REAL(4) and REAL(8).)
>
>
> Unfortunately, I'm not that familiar with middle-end translation yet,
> so it's not immediately clear to me why the calling code is passing
> by-value instead of by-reference to the library function. I don't see
> significant differences between the way these functions are wrapped
> compared to other libgfortran routines which do pass by-reference. The
> biggest difference I see is the LIB_FUNCTION() lines in the
> gfc_intrinsic_map in trans-intrinsic.c: why are these present for the
> trigd functions and not for other libgfortran functions?
>
> $ grep -C4 'LIB_FUNCTION (SIND' trans-intrinsic.c | nl -ba -v 119
>    119    #include "mathbuiltins.def"
>    120
>    121      /* Functions in libgfortran.  */
>    122      LIB_FUNCTION (ERFC_SCALED, "erfc_scaled", false),
>    123      LIB_FUNCTION (SIND, "sind", false),
>    124      LIB_FUNCTION (COSD, "cosd", false),
>    125      LIB_FUNCTION (TAND, "tand", false),
>    126
>    127      /* End the list.  */
>
>
> Steve, do you have any ideas? I will continue to poke around.
>

I may have misread the -fdump-tree-original as I was trying
to get things to work.  After you apply the patch, you should
have a libgfortran/intrinsics/trigd_inc.c file.  It can probaby
be adjusted to take pass value.  Does this patch fix things?
 
--- trigd_inc.c.orig 2020-03-16 17:44:02.756817000 -0700
+++ trigd_inc.c 2020-03-16 17:45:40.204464000 -0700
@@ -26,30 +26,30 @@
 
 /* Compute sind(x) = sin(x * pi / 180). */
 
-extern WTYPE SIND (WTYPE *);
+extern WTYPE SIND (WTYPE);
 export_proto(SIND);
 
-extern WTYPE COSD (WTYPE *);
+extern WTYPE COSD (WTYPE);
 export_proto(COSD);
 
-extern WTYPE TAND (WTYPE *);
+extern WTYPE TAND (WTYPE);
 export_proto(TAND);
 
 WTYPE
-SIND (WTYPE *x)
+SIND (WTYPE x)
 {
   int n;
   WTYPE ax, s;
 
-  if (isfinite(*x))
+  if (isfinite(x))
     {
       /* sin(-x) = - sin(x).  */
-      s = CPYSGN(*x);
-      ax = FABS(*x);
+      s = CPYSGN(x);
+      ax = FABS(x);
 
       /* In radians, sin(x) = x as x -> 0.  */
       if (ax < SMALL)
- return (*x * PIO180);
+ return (x * PIO180);
 
       /* Reduce angle to ax in [0,360].  */
       ax = FMOD(ax);
@@ -93,7 +93,7 @@
 /* Compute cosd(x) = cos(x * pi / 180).  */
 
 WTYPE
-COSD (WTYPE *x)
+COSD (WTYPE x)
 {
 #ifdef TINYF
   static const volatile WTYPE tiny = 1.e-30f;
@@ -106,9 +106,9 @@
   int n;
   WTYPE ax;
 
-  if (isfinite(*x))
+  if (isfinite(x))
     {
-      ax = FABS(*x);
+      ax = FABS(x);
 
       /* No spurious underflows!.  In radians, cos(x) = 1-x*x/2 as x -> 0.  */
       if (ax < SMALL)
@@ -156,20 +156,20 @@
 /* Compute tand(x) = tan(x * pi / 180).  */
 
 WTYPE
-TAND (WTYPE *x)
+TAND (WTYPE x)
 {
   int n;
   WTYPE ax, s;
 
-  if (isfinite(*x))
+  if (isfinite(x))
     {
       /* tan(-x) = - tan(x). */
-      s = CPYSGN(*x);
-      ax = FABS(*x);
+      s = CPYSGN(x);
+      ax = FABS(x);
 
       /* In radians, tan(x) = x as x -> 0. */
       if (ax < SMALL)
- return (*x * PIO180);
+ return (x * PIO180);
 
       ax = FMOD(ax);
 

--
Steve
Reply | Threaded
Open this post in threaded view
|

Re: PATCH -- Fix degree trignometric functions

gcc - fortran mailing list
On Mon, Mar 16, 2020 at 8:55 PM Steve Kargl
<[hidden email]> wrote:
>
> On Mon, Mar 16, 2020 at 08:27:25PM -0400, Fritz Reese wrote:
[...]
> > > > On Fri, Mar 06, 2020 at 03:18:19PM -0800, Steve Kargl wrote:
> > [...]
> > > > >   5. New functions have been added to libgfortran to handle sind, cosd,
> > > > >      and tand.
[...]
> > However, after applying your patch, the new libgfortran functions
> > appear to have a calling-convention mismatch between the frontend and
> > the library
[...]
> I may have misread the -fdump-tree-original as I was trying
> to get things to work.  After you apply the patch, you should
> have a libgfortran/intrinsics/trigd_inc.c file.  It can probaby
> be adjusted to take pass value.  Does this patch fix things?
[...]

Yes, this works. I'd still like to know what determines PBV versus PBR
-- does this have to do with CLASS_ELEMENTAL and PURE-ity?

I have a few other questions about the library functions. Please
forgive any silly questions which may be simply due to ignorance.

For COSD from trigd_inc.c...

WTYPE
COSD (WTYPE x)
{
#ifdef TINYF
  static const volatile WTYPE tiny = 1.e-30f;
#elif TINY
  static const volatile WTYPE tiny = 1.e-300f;
#else
  static const volatile WTYPE tiny = 0x1p-10000L;
#endif
[...]
      /* No spurious underflows!.  In radians, cos(x) = 1-x*x/2 as x -> 0.  */
      if (ax < SMALL)
        return (ax == 0 ? 1 : 1 - tiny);

AFAICT neither TINYF nor TINY are defined by trigd.c. Even if they
were, 1.e-300f is a single precision literal and will always be zero
-- this is probably meant to be 1.e-300 without the single precision
suffix. Why use this method to allegedly avoid underflow in COSD but
not SIND or TAND? SIND uses x * PIO180 -- and why use multiplication
here instead of FMA like in D2R(x)?

Additionally, from gfortran/intrinsics/trigd.c:

#ifdef HAVE_GFC_REAL_10
[...]
#define SMALL       0x1.p-27L
#ifdef HAVE_GFC_REAL_16
[...]
#ifdef GFC_REAL_16_IS_FLOAT128   /* libquadmath.  */
#define SMALL       3.4e-17Q
[...]
#else
#define SMALL       3.4e-17L
[...]

Is there a reason to use base 10 for the r16 versions, as opposed to
the base 2 notation used by r4, r8, and r10?

Regarding the angle wrapping - I've done a comparison against the libm
implementation, and it appears there are a few inconsistencies and
mishandled edge conditions. I used the attached program
trigd_compare.f90 to check these -- output attached both for your
original patch (trigd1.out) and the next version (trigd2.out). I also
attached a diff (trigd2.patch) which is applied on top of your
pass-by-value-patched version. This includes some formatting fixes,
fixes for TINY, and fixes for below. Please glance over it and let me
know if you think the changes look okay.

* cosd(180) gives +1, should be -1.
    - Fixed in the attached diff so cosd(180) returns -1.

* cos(~ 3pi/2) gives ~ -0, whereas cosd(270) gives +0.
    - Fixed in the attached diff so cosd(270) returns -0.

* sin(~ -2pi) gives ~ +0, sin(~ 2pi) gives ~ -0; whereas sind(-360)
gives -0 and sind(360) gives +0.
    - I've not changed this in the attached diff, do you think this matters?

* After fixing the above sign issues, tand has opposite sign for every
multiple of 180 other than 0, compared to tan() for multiples of pi.
    - tand(+-180) -> +- 0, where tan(+-pi) -> -+0
    - tand(+-360) -> +- 0, where tan(+-pi) -> -+0
    - I've not changed this in the attached diff, do you think this matters?

As for speed, I was curious so I ran some primitive benchmarks. SIND
appears to be within +-10% of libm's SIN speed for both r4 and r8, so
that's good news.

---
Fritz

trigd2.out (55K) Download Attachment
trigd_compare.f90 (6K) Download Attachment
trigd1.out (55K) Download Attachment
trigd2.patch (19K) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: PATCH -- Fix degree trignometric functions

gcc - fortran mailing list
On Wed, Mar 18, 2020 at 05:31:49PM -0400, Fritz Reese wrote:
> On Mon, Mar 16, 2020 at 8:55 PM Steve Kargl
> <[hidden email]> wrote:
> >

(deleting stuff to shorten to only pertinent info).

> > I may have misread the -fdump-tree-original as I was trying
> > to get things to work.  After you apply the patch, you should
> > have a libgfortran/intrinsics/trigd_inc.c file.  It can probaby
> > be adjusted to take pass value.  Does this patch fix things?
>
> Yes, this works. I'd still like to know what determines PBV versus PBR
> -- does this have to do with CLASS_ELEMENTAL and PURE-ity?

Like you, I struggle with the trans-* files.  If you look in
trans-intrinsic.c, you'll find gfc_conv_intrinsic_function(),
which contains a large switch statement that controls the
conversion of function calls.  GFC_ISYM_SIND, GFC_ISYM_COSD,
and GFC_ISYM_TANH are not explicitly handles, so we hit the
default.  The default case simply calls gfc_conv_intrinsic_lib_function()
to generate a library function call, and that call is assumed
to follow C semantics.  The signatures for these functions are
built by lines 123-125, which invokes the LIB_FUNCTION macro.
From here, it's magic to me. :-)


> I have a few other questions about the library functions. Please
> forgive any silly questions which may be simply due to ignorance.
>
> For COSD from trigd_inc.c...
>
> WTYPE
> COSD (WTYPE x)
> {
> #ifdef TINYF
>   static const volatile WTYPE tiny = 1.e-30f;
> #elif TINY
>   static const volatile WTYPE tiny = 1.e-300f;
> #else
>   static const volatile WTYPE tiny = 0x1p-10000L;
> #endif

The functions in trigd_inc.c evolved from 3 seperate C files, where
I developed everything for float and then used macros to test double
and long double.  The above was used to get 1 - tiny = 1 for x << 1.
tiny needs to be volatile to prevent gcc from constant-folding the
subtraction, which is needed to raise the FE_INEXACT exception.
Unfortunately, 1.e-30f isn't small enough for IEEE128 floating point
so we need at least two values.  This can be changed to  

#ifdef TINYF
   static const volatile WTYPE tiny = 1.e-30f;
#else
   static const volatile WTYPE tiny = 1.e-300f;
#endif

and TINYF should be defined in trigd.c for GFC_REAL_4
and undefined for the other types.

> [...]
>       /* No spurious underflows!.  In radians, cos(x) = 1-x*x/2 as x -> 0.  */
>       if (ax < SMALL)
>         return (ax == 0 ? 1 : 1 - tiny);
>
> AFAICT neither TINYF nor TINY are defined by trigd.c. Even if they
> were, 1.e-300f is a single precision literal and will always be zero
> -- this is probably meant to be 1.e-300 without the single precision
> suffix. Why use this method to allegedly avoid underflow in COSD but
> not SIND or TAND? SIND uses x * PIO180 -- and why use multiplication
> here instead of FMA like in D2R(x)?

cosd(x) = 1 as x -> 0 without underflow.  For a subnormal number
and normal numbers near the transition from normal to subnormal,
cosd(x) = 1 - x*x/2 will raise FE_UNDERFLOW.  We want to avoid
that signal, but we should FE_INEXACT.

For sind(x) and tand(x), these have sind(x) = tand(x) = x * pi / 180.
So, if FE_UNDERFLOW is raised (and I think it is but haven't
checked by x * PIO180), then that's ok.

>
> Additionally, from gfortran/intrinsics/trigd.c:
>
> #ifdef HAVE_GFC_REAL_10
> [...]
> #define SMALL       0x1.p-27L
> #ifdef HAVE_GFC_REAL_16
> [...]
> #ifdef GFC_REAL_16_IS_FLOAT128   /* libquadmath.  */
> #define SMALL       3.4e-17Q
> [...]
> #else
> #define SMALL       3.4e-17L
> [...]
>
> Is there a reason to use base 10 for the r16 versions, as opposed to
> the base 2 notation used by r4, r8, and r10?

Laziness.  SMALL is a sloppy threshold where for x*x/2 < EPSILON,
so needs to be something like 2**(p/2) with p the precision.  For
speed, these should be as large as possible to give cosd(x) = 1 =
1 - tiny with FE_INEXACT to avoid a multiplication and division.

> Regarding the angle wrapping - I've done a comparison against the libm
> implementation, and it appears there are a few inconsistencies and
> mishandled edge conditions.

Which libm?  gfortran should agree with DEC Fortran, but I don't
have access to that hardware/software.  Your fixes below are
probably correct.

> I used the attached program
> trigd_compare.f90 to check these -- output attached both for your
> original patch (trigd1.out) and the next version (trigd2.out). I also
> attached a diff (trigd2.patch) which is applied on top of your
> pass-by-value-patched version. This includes some formatting fixes,
> fixes for TINY, and fixes for below. Please glance over it and let me
> know if you think the changes look okay.
>
> * cosd(180) gives +1, should be -1.
>     - Fixed in the attached diff so cosd(180) returns -1.
>
> * cos(~ 3pi/2) gives ~ -0, whereas cosd(270) gives +0.
>     - Fixed in the attached diff so cosd(270) returns -0.
>
> * sin(~ -2pi) gives ~ +0, sin(~ 2pi) gives ~ -0; whereas sind(-360)
> gives -0 and sind(360) gives +0.
>     - I've not changed this in the attached diff, do you think this matters?

I have been in an email exchange with someone much smarter
than I when it comes to floating point.  He suggested that for
sin(x), and this should apply for sind(x), that it may be
appropriate to have sin(n*pi) = 0 for even n and sin(n*pi)=-0
for odd.  I think he based this on the derivative of sin(x)
is cos(x) and n*pi gives 1 for even and -1 for odd.  For Fortran,
it likely doesn't matter as +-0 is normally treated as just 0.
If you want to pursue this, in SIND(WTYPE x) you can probably do

      /* Special cases.  */
      n = (int)ax;
      if (n%30 == 0 && ax == n)
        {
           s = s * (n & 1 ? -1 : 1);  /* Flip sign for n*pi even/odd. */

or some variation.


> * After fixing the above sign issues, tand has opposite sign for every
> multiple of 180 other than 0, compared to tan() for multiples of pi.
>     - tand(+-180) -> +- 0, where tan(+-pi) -> -+0
>     - tand(+-360) -> +- 0, where tan(+-pi) -> -+0
>     - I've not changed this in the attached diff, do you think this matters?

I could have gotten the symmetry wrong for special cases.
We could probably choose a sign based on n&1 (see above discussion).

> As for speed, I was curious so I ran some primitive benchmarks. SIND
> appears to be within +-10% of libm's SIN speed for both r4 and r8, so
> that's good news.

I did put effort into making this fast and correct.  Unfortunately,
fmod() is going to do argument checking (e.g., for +-inf and NaN),
which the these routine also do, and the conversion from degree
to radian uses fma().  We could use a poor man's rint by replacing
ax = fmod(ax) with ax = ax + 0x1.pN - 0x1.p-N where N > P is chosen
appropriately for each type.  This only works for ax < 0x1.pP where
P is precision.

--
steve
Reply | Threaded
Open this post in threaded view
|

Re: PATCH -- Fix degree trignometric functions

gcc - fortran mailing list
On Wed, Mar 18, 2020 at 6:58 PM Steve Kargl
<[hidden email]> wrote:
>
> On Wed, Mar 18, 2020 at 05:31:49PM -0400, Fritz Reese wrote:
> > On Mon, Mar 16, 2020 at 8:55 PM Steve Kargl
[...]
> [...] If you look in
> trans-intrinsic.c, you'll find gfc_conv_intrinsic_function(),
> [...]

Ah, thanks!

> [...] This can be changed to
>
> #ifdef TINYF
>    static const volatile WTYPE tiny = 1.e-30f;
> #else
>    static const volatile WTYPE tiny = 1.e-300f;
> #endif
>
> and TINYF should be defined in trigd.c for GFC_REAL_4
> and undefined for the other types.
[...]
> cosd(x) = 1 as x -> 0 without underflow.  For a subnormal number
> and normal numbers near the transition from normal to subnormal,
> cosd(x) = 1 - x*x/2 will raise FE_UNDERFLOW.  We want to avoid
> that signal, but we should FE_INEXACT.

Ah, now it clicks. I mixed the definitions of UNDERFLOW versus INEXACT
in my mind. That makes sense.


> For sind(x) and tand(x), these have sind(x) = tand(x) = x * pi / 180.
> So, if FE_UNDERFLOW is raised (and I think it is but haven't
> checked by x * PIO180), then that's ok.

What about using sind(x) = tand(x) = D2R(x) instead of x * PIO180, as
is done elsewhere?

> [...] SMALL is a sloppy threshold where for x*x/2 < EPSILON,
> so needs to be something like 2**(p/2) with p the precision.  For
> speed, these should be as large as possible to give cosd(x) = 1 =
> 1 - tiny with FE_INEXACT to avoid a multiplication and division.


> > Regarding the angle wrapping - I've done a comparison against the libm
> > implementation, and it appears there are a few inconsistencies and
> > mishandled edge conditions.
>
> Which libm?  gfortran should agree with DEC Fortran, but I don't
> have access to that hardware/software.  Your fixes below are
> probably correct.

This is with libm from glibc-2.17.el7.x86_64.

> I have been in an email exchange with someone much smarter
> than I when it comes to floating point.  He suggested that for
> sin(x), and this should apply for sind(x), that it may be
> appropriate to have sin(n*pi) = 0 for even n and sin(n*pi)=-0
> for odd.  I think he based this on the derivative of sin(x)
> is cos(x) and n*pi gives 1 for even and -1 for odd.  For Fortran,
> it likely doesn't matter as +-0 is normally treated as just 0.
> If you want to pursue this, in SIND(WTYPE x) you can probably do
>
>       /* Special cases.  */
>       n = (int)ax;
>       if (n%30 == 0 && ax == n)
>         {
>            s = s * (n & 1 ? -1 : 1);  /* Flip sign for n*pi even/odd. */
>
> or some variation.

Sounds good, I'll implement that.

> I did put effort into making this fast and correct.  Unfortunately,
> fmod() is going to do argument checking (e.g., for +-inf and NaN),
> which the these routine also do, and the conversion from degree
> to radian uses fma().  We could use a poor man's rint by replacing
> ax = fmod(ax) with ax = ax + 0x1.pN - 0x1.p-N where N > P is chosen
> appropriately for each type.  This only works for ax < 0x1.pP where
> P is precision.

I'll leave this as an exercise for the future if somebody is unhappy
with the performance.

I'll send another draft shortly.

---
Fritz
Reply | Threaded
Open this post in threaded view
|

Re: PATCH -- Fix degree trignometric functions

gcc - fortran mailing list
On Thu, Mar 19, 2020 at 02:59:05PM -0400, Fritz Reese wrote:
> On Wed, Mar 18, 2020 at 6:58 PM Steve Kargl
> <[hidden email]> wrote:
> >
> > #ifdef TINYF
> >    static const volatile WTYPE tiny = 1.e-30f;
> > #else
> >    static const volatile WTYPE tiny = 1.e-300f;

Note, I kept the typo!   1.e-300f should not have
the f suffice.  It should be 1.e-300.

> > For sind(x) and tand(x), these have sind(x) = tand(x) = x * pi / 180.
> > So, if FE_UNDERFLOW is raised (and I think it is but haven't
> > checked by x * PIO180), then that's ok.
>
> What about using sind(x) = tand(x) = D2R(x) instead of x * PIO180, as
> is done elsewhere?

This could be done, but leads to a call to fma().  So,
the simple multiplication is a micro-optimization.

> > [...] SMALL is a sloppy threshold where for x*x/2 < EPSILON,
> > so needs to be something like 2**(p/2) with p the precision.  For

Again a typo, should have have been 2**(-p/2), e.g., 2**(-12)
for REAL_4.  The threshold is sloppy, because at least with
REAL_4 one can explicit search for the largest value that
SMALL can be to still get cosd(x) = 1 for x << 1.

--
Steve
Reply | Threaded
Open this post in threaded view
|

Re: PATCH -- Fix degree trignometric functions

gcc - fortran mailing list
On Thu, Mar 19, 2020 at 3:36 PM Steve Kargl
<[hidden email]> wrote:

>
> On Thu, Mar 19, 2020 at 02:59:05PM -0400, Fritz Reese wrote:
> > On Wed, Mar 18, 2020 at 6:58 PM Steve Kargl
> > <[hidden email]> wrote:
> > >
> > > #ifdef TINYF
> > >    static const volatile WTYPE tiny = 1.e-30f;
> > > #else
> > >    static const volatile WTYPE tiny = 1.e-300f;
>
> Note, I kept the typo!   1.e-300f should not have
> the f suffice.  It should be 1.e-300.

Right.

> > > For sind(x) and tand(x), these have sind(x) = tand(x) = x * pi / 180.
> > > So, if FE_UNDERFLOW is raised (and I think it is but haven't
> > > checked by x * PIO180), then that's ok.
> >
> > What about using sind(x) = tand(x) = D2R(x) instead of x * PIO180, as
> > is done elsewhere?
>
> This could be done, but leads to a call to fma().  So,
> the simple multiplication is a micro-optimization.

OK.



> > > [...] SMALL is a sloppy threshold where for x*x/2 < EPSILON,
> > > so needs to be something like 2**(p/2) with p the precision.  For
>
> Again a typo, should have have been 2**(-p/2), e.g., 2**(-12)
> for REAL_4.  The threshold is sloppy, because at least with
> REAL_4 one can explicit search for the largest value that
> SMALL can be to still get cosd(x) = 1 for x << 1.

Some tests show 1 is the closest representable value for cosd(x)
without loss of precision when x is less than 2**-{7, 21, 26, 51} for
real({4,8,10,16}).

This threshold (SMALL) is also used for sind(x) ~= x for x << 1 (->
0). This is true without loss of precision when x < 2**-5 for real(4)
and x < 2**-20 for real(8), but for with real(10) and real(16) there
are no values for which this is true to within one (or even several
ULP). There is a loss of precision after about the 17th place (out of
32) when approximating with multiplication for x < 2**-20. For
comparison, [an old version of] Intel's DEC Fortran compiler appears
to favor precision (for abitrarily small x). I think we should mimic
this behavior.

I'll post a revised patch soon.

---
Fritz
Reply | Threaded
Open this post in threaded view
|

Re: PATCH -- Fix degree trignometric functions

gcc - fortran mailing list
On Tue, Mar 24, 2020 at 03:45:29PM -0400, Fritz Reese wrote:

> I wrote:
> >
> > Again a typo, should have have been 2**(-p/2), e.g., 2**(-12)
> > for REAL_4.  The threshold is sloppy, because at least with
> > REAL_4 one can explicit search for the largest value that
> > SMALL can be to still get cosd(x) = 1 for x << 1.
>
> Some tests show 1 is the closest representable value for cosd(x)
> without loss of precision when x is less than 2**-{7, 21, 26, 51} for
> real({4,8,10,16}).
>
> This threshold (SMALL) is also used for sind(x) ~= x for x << 1 (->
> 0). This is true without loss of precision when x < 2**-5 for real(4)
> and x < 2**-20 for real(8), but for with real(10) and real(16) there
> are no values for which this is true to within one (or even several
> ULP). There is a loss of precision after about the 17th place (out of
> 32) when approximating with multiplication for x < 2**-20. For
> comparison, [an old version of] Intel's DEC Fortran compiler appears
> to favor precision (for abitrarily small x). I think we should mimic
> this behavior.
>

I would need to redo some testing for sind(x) with x << 1.  If one is
trying to avoid spurious floating exceptions, one can use

      if (ax < SMALL)
        return (sin(D2R(x)));
or
      if (ax < SMALL)
        {
          WTYPE z = D2R(x);
          return (z * (1 - z * a) * (1 + z * a));   /* a = sqrt(1/6.)
        }

The former has a speed penalty from sin() checking for +-inf and NaN
and dealing with the sign of its argument.  The latter uses the 1st
2 terms of sin(x) = x - x**3/6 = x*(1-x**2/6).

Thanks for carry this to a finish line.

--
Steve
Reply | Threaded
Open this post in threaded view
|

Re: PATCH -- Fix degree trignometric functions

gcc - fortran mailing list
I've discovered my DEC compiler (Intel's ifort, albeit an old version)
has the following behaviors:

* No trig intrinsics ever return negative zero -- this includes
radian-valued as well as degree-valued intrinsics.

I don't think we need to sweat over this too much, since glibc's
radian-trig functions return negative zeroes, and Fortran doesn't
usually care.

* TAND(x) returns sgn(x) * +Infinity for x % 360 = 90, and returns
sgn(x) * -Infinity for x % 360 = 270 -- our current implementation
returns NaN for these cases.

Though NaN is mathematically correct, this result allows 1 /
tand(+-90) to return zero instead of NaN while providing +/- symmetry
across the number line. Perhaps we should follow suit.

* COTAND(x) returns sgn(x) * +Infinity for x % 180 == 0, and returns
sgn(x) * -Infinity for x % 180 == 90 -- our current implementation
returns a large number.

The DEC behavior here seems preferable. If we change COTAND to fold as
1/TAND instead of 1/TAN(x*pi/180), we should at least get Inf for free
at 180 since TAND(180) === 0. Then we would need to change TAND as
above to return +/- Inf rather than NaN, so that COTAND(90) === 0
instead of NaN. This should be easier than implementing
_gfortran_cotand_*.

Thoughts?

---
Fritz Reese
Reply | Threaded
Open this post in threaded view
|

Re: PATCH -- Fix degree trignometric functions

gcc - fortran mailing list
On Wed, Mar 25, 2020 at 07:58:50PM -0400, Fritz Reese wrote:
> I've discovered my DEC compiler (Intel's ifort, albeit an old version)
> has the following behaviors:
>
> * No trig intrinsics ever return negative zero -- this includes
> radian-valued as well as degree-valued intrinsics.
>
> I don't think we need to sweat over this too much, since glibc's
> radian-trig functions return negative zeroes, and Fortran doesn't
> usually care.

sin(-0) should return -0, so I think sind(-0) should also
return -0.  For x near n*pi, sin(x) will never return 0 as
n*pi is not exactly representable.  For sind(x) and x = n * 180,
we should follow what DEC does (if we can determin the behavior).


> * TAND(x) returns sgn(x) * +Infinity for x % 360 = 90, and returns
> sgn(x) * -Infinity for x % 360 = 270 -- our current implementation
> returns NaN for these cases.
>
> Though NaN is mathematically correct, this result allows 1 /
> tand(+-90) to return zero instead of NaN while providing +/- symmetry
> across the number line. Perhaps we should follow suit.

This one is tricky. n * pi / 2 is not exactly representable so
tan(x) will never hit the NaN condition (ie., x is either larger
or smaller than n*pi/2, so +-inf seems correct).  For tand(x),
x = n * 90 is exactly representable, so we are exactly on the
asymptote.  To me, NaN seems sensible.  But, again, gfortran
should probably do what DEC does.  Afterall, this is a compatibility
extension.

> * COTAND(x) returns sgn(x) * +Infinity for x % 180 == 0, and returns
> sgn(x) * -Infinity for x % 180 == 90 -- our current implementation
> returns a large number.
>
> The DEC behavior here seems preferable. If we change COTAND to fold as
> 1/TAND instead of 1/TAN(x*pi/180), we should at least get Inf for free
> at 180 since TAND(180) === 0. Then we would need to change TAND as
> above to return +/- Inf rather than NaN, so that COTAND(90) === 0
> instead of NaN. This should be easier than implementing
> _gfortran_cotand_*.

For cotand(), I in-lined it with gfc_conv_intrinsic_cotand().  We
can change the in-lining to generate 1/tand(x).  Are you up for
changing trans-intrinsic.c or would you like me to do the change?

--
Steve
Reply | Threaded
Open this post in threaded view
|

Re: PATCH -- Fix degree trignometric functions

gcc - fortran mailing list
On Wed, Mar 25, 2020 at 9:21 PM Steve Kargl
<[hidden email]> wrote:
>
> On Wed, Mar 25, 2020 at 07:58:50PM -0400, Fritz Reese wrote:
> > I've discovered my DEC compiler (Intel's ifort, albeit an old version)
> > has the following behaviors:
[...]



> sin(-0) should return -0, so I think sind(-0) should also
> return -0.  For x near n*pi, sin(x) will never return 0 as
> n*pi is not exactly representable.  For sind(x) and x = n * 180,
> we should follow what DEC does (if we can determin the behavior).

DEC always returns positive zero -- even for SIN(-0) and SIND(-0). But
I prefer to follow SIND(-0) = -0 and SIND(n*180) based on the parity
of n as we've previously established.

> For cotand(), I in-lined it with gfc_conv_intrinsic_cotand().  We
> can change the in-lining to generate 1/tand(x).  Are you up for
> changing trans-intrinsic.c or would you like me to do the change?

I'll take care of it. Thanks for the offer. Will post an updated patch soon.

---
Fritz
Reply | Threaded
Open this post in threaded view
|

Re: PATCH -- Fix degree trignometric functions

gcc - fortran mailing list
In reply to this post by Steve Kargl
On Fri, Mar 6, 2020 at 6:18 PM Steve Kargl
<[hidden email]> wrote:
[...]
> TL;DR version.
>
>   Fix the simplification and handling of the degree trigonometric functions.
>   This includes fixing a number of ICEs.  See PR 93871.

An updated version of the patch is attached. Regression tests (and new
test cases) are pending.

Changes since Steve's patch of note:

* libgfortran/intrinsics/trigd.c now indirectly includes trigd.inc
(formerly trigd_inc.c) through trigd_lib.inc. trigd.inc is now written
using GMP/MPFR functions, so that the same file can be included from
the front-end for simplification. The GMP/MPFR functions are
pre-processed into native code by trigd_lib.inc. This ensures that
both the FE and the library are using the same code for resolving
these functions, preventing future maintenance woes.

* TAND(90, 270) returns +Infinity and -Infinity, instead of NaN. This
is compatible with behavior of (at least some) DEC compilers.

* COTAN[D] are implemented as -TAN[D](x + 90 degrees) rather than 1 /
COTAN[D] for speed and to avoid singularities (though the new
implementation of TAND allows 1 / TAND to return zero when TAND
returns infinity).

* SMALL thresholds and signs of some specific values are corrected.
For REAL(10) and REAL(16), the SMALL threshold for SIND(x) = D2R(x) is
eliminated, since there would be loss of precision. The other
thresholds achieve COSD(x) = 1 and SIND(x) = D2R(x) without loss of
precision.

ChangeLogs are below. I will post an update soon after I perform
regression tests, which include some new tests which I will add.

gcc/fortran/ChangeLog:

        PR fortran/93871

        * gfortran.h (GFC_ISYM_ACOSD, GFC_ISYM_ASIND, GFC_ISYM_ATAN2D,
        GFC_ISYM_ATAND, GFC_ISYM_COSD, GFC_ISYM_COTAND, GFC_ISYM_SIND,
        GFC_ISYM_TAND): New.
        * intrinsic.c (add_functions): Remove check for flag_dec_math.
        Give degree trig functions simplification and name resolution
        functions (e.g, gfc_simplify_atrigd () and gfc_resolve_atrigd ()).
        (do_simplify): Remove special casing of degree trig functions.
        * intrinsic.h (gfc_simplify_acosd, gfc_simplify_asind,
        gfc_simplify_atand, gfc_simplify_cosd, gfc_simplify_cotand,
        gfc_simplify_sind, gfc_simplify_tand, gfc_resolve_trigd2): Add new
        prototypes.
        (gfc_simplify_atrigd, gfc_simplify_trigd, gfc_resolve_cotan,
        resolve_atrigd): Remove prototypes of deleted functions.
        * iresolve.c (is_trig_resolved, copy_replace_function_shallow,
        gfc_resolve_cotan, get_radians, get_degrees, resolve_trig_call,
        gfc_resolve_atrigd, gfc_resolve_atan2d): Delete functions.
        (gfc_resolve_trigd, gfc_resolve_trigd2): Resolve to library functions.
        * simplify.c (rad2deg, deg2rad, gfc_simplify_acosd, gfc_simplify_asind,
        gfc_simplify_atand, gfc_simplify_atan2d, gfc_simplify_cosd,
        gfc_simplify_sind, gfc_simplify_tand, gfc_simplify_cotand): New
        functions.
        (gfc_simplify_atan2): Fix error message.
        (simplify_trig_call, gfc_simplify_trigd, gfc_simplify_atrigd,
        radians_f): Delete functions.
        * trans-intrinsic.c: Add LIB_FUNCTION decls for sind, cosd, tand.
        (rad2deg, gfc_conv_intrinsic_atrigd, gfc_conv_intrinsic_cotan,
        gfc_conv_intrinsic_cotand, gfc_conv_intrinsic_atan2d): New functions.
        (gfc_conv_intrinsic_function): Handle ACOSD, ASIND, ATAND, COTAN,
        COTAND, ATAN2D.
        * trigd_fe.inc: New file. Included by simplify.c to implement
        simplify_sind, simplify_cosd, simplify_tand with code common to the
        libgfortran implementation.

libgfortran/ChangeLog:

        PR fortran/93871
        * Makefile.am, Makefile.in: New make rule for intrinsics/trigd.c.
        * gfortran.map: New routines for {sind, cosd, tand}X{r4, r8, r10, r16}.
        * intrinsics/trigd.c, intrinsics/trigd_lib.inc, intrinsics/trigd.inc:
        New files. Defines native degree-valued trig functions.

trigd_v3.patch (93K) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: PATCH -- Fix degree trignometric functions

Tobias Burnus
Two remarks:

* As the file trigd_lib.inc is shared between libgfortran
   and gcc/fortran, I wonder whether it should be placed
   under include/ (under the GCC toplevel directroy)

* All included files need dependency; I do not quickly
   see whether that' the case.

Cheers,

Tobias

On 3/28/20 12:36 AM, Fritz Reese via Fortran wrote:

> On Fri, Mar 6, 2020 at 6:18 PM Steve Kargl
> <[hidden email]> wrote:
> [...]
>> TL;DR version.
>>
>>    Fix the simplification and handling of the degree trigonometric functions.
>>    This includes fixing a number of ICEs.  See PR 93871.
> An updated version of the patch is attached. Regression tests (and new
> test cases) are pending.
>
> Changes since Steve's patch of note:
>
> * libgfortran/intrinsics/trigd.c now indirectly includes trigd.inc
> (formerly trigd_inc.c) through trigd_lib.inc. trigd.inc is now written
> using GMP/MPFR functions, so that the same file can be included from
> the front-end for simplification. The GMP/MPFR functions are
> pre-processed into native code by trigd_lib.inc. This ensures that
> both the FE and the library are using the same code for resolving
> these functions, preventing future maintenance woes.
>
> * TAND(90, 270) returns +Infinity and -Infinity, instead of NaN. This
> is compatible with behavior of (at least some) DEC compilers.
>
> * COTAN[D] are implemented as -TAN[D](x + 90 degrees) rather than 1 /
> COTAN[D] for speed and to avoid singularities (though the new
> implementation of TAND allows 1 / TAND to return zero when TAND
> returns infinity).
>
> * SMALL thresholds and signs of some specific values are corrected.
> For REAL(10) and REAL(16), the SMALL threshold for SIND(x) = D2R(x) is
> eliminated, since there would be loss of precision. The other
> thresholds achieve COSD(x) = 1 and SIND(x) = D2R(x) without loss of
> precision.
>
> ChangeLogs are below. I will post an update soon after I perform
> regression tests, which include some new tests which I will add.
>
> gcc/fortran/ChangeLog:
>
>          PR fortran/93871
>
>          * gfortran.h (GFC_ISYM_ACOSD, GFC_ISYM_ASIND, GFC_ISYM_ATAN2D,
>          GFC_ISYM_ATAND, GFC_ISYM_COSD, GFC_ISYM_COTAND, GFC_ISYM_SIND,
>          GFC_ISYM_TAND): New.
>          * intrinsic.c (add_functions): Remove check for flag_dec_math.
>          Give degree trig functions simplification and name resolution
>          functions (e.g, gfc_simplify_atrigd () and gfc_resolve_atrigd ()).
>          (do_simplify): Remove special casing of degree trig functions.
>          * intrinsic.h (gfc_simplify_acosd, gfc_simplify_asind,
>          gfc_simplify_atand, gfc_simplify_cosd, gfc_simplify_cotand,
>          gfc_simplify_sind, gfc_simplify_tand, gfc_resolve_trigd2): Add new
>          prototypes.
>          (gfc_simplify_atrigd, gfc_simplify_trigd, gfc_resolve_cotan,
>          resolve_atrigd): Remove prototypes of deleted functions.
>          * iresolve.c (is_trig_resolved, copy_replace_function_shallow,
>          gfc_resolve_cotan, get_radians, get_degrees, resolve_trig_call,
>          gfc_resolve_atrigd, gfc_resolve_atan2d): Delete functions.
>          (gfc_resolve_trigd, gfc_resolve_trigd2): Resolve to library functions.
>          * simplify.c (rad2deg, deg2rad, gfc_simplify_acosd, gfc_simplify_asind,
>          gfc_simplify_atand, gfc_simplify_atan2d, gfc_simplify_cosd,
>          gfc_simplify_sind, gfc_simplify_tand, gfc_simplify_cotand): New
>          functions.
>          (gfc_simplify_atan2): Fix error message.
>          (simplify_trig_call, gfc_simplify_trigd, gfc_simplify_atrigd,
>          radians_f): Delete functions.
>          * trans-intrinsic.c: Add LIB_FUNCTION decls for sind, cosd, tand.
>          (rad2deg, gfc_conv_intrinsic_atrigd, gfc_conv_intrinsic_cotan,
>          gfc_conv_intrinsic_cotand, gfc_conv_intrinsic_atan2d): New functions.
>          (gfc_conv_intrinsic_function): Handle ACOSD, ASIND, ATAND, COTAN,
>          COTAND, ATAN2D.
>          * trigd_fe.inc: New file. Included by simplify.c to implement
>          simplify_sind, simplify_cosd, simplify_tand with code common to the
>          libgfortran implementation.
>
> libgfortran/ChangeLog:
>
>          PR fortran/93871
>          * Makefile.am, Makefile.in: New make rule for intrinsics/trigd.c.
>          * gfortran.map: New routines for {sind, cosd, tand}X{r4, r8, r10, r16}.
>          * intrinsics/trigd.c, intrinsics/trigd_lib.inc, intrinsics/trigd.inc:
>          New files. Defines native degree-valued trig functions.
Reply | Threaded
Open this post in threaded view
|

Re: PATCH -- Fix degree trignometric functions

gcc - fortran mailing list
On Sat, Mar 28, 2020 at 08:10:38AM +0100, Tobias Burnus wrote:
> Two remarks:
>
> * As the file trigd_lib.inc is shared between libgfortran
>   and gcc/fortran, I wonder whether it should be placed
>   under include/ (under the GCC toplevel directroy)

The original name and location were chosen to match
intrinsics/erfc_scaled_inc.c, which is included in
intrinsics/erfc_scaled.c.  I think Fritz's re-use
in simplification makes since given the ugly nested
if-elseif-else constructs in the file.  As this is
Fortran specific, I do not believe it should be moved
up to a toplevel diretory.

>
> * All included files need dependency; I do not quickly
>   see whether that' the case.


They do?  See intrinsics/erfc_scaled*.c.

% grep erfc_ Makefile.am
intrinsics/erfc_scaled.c \

--
steve
Reply | Threaded
Open this post in threaded view
|

Re: PATCH -- Fix degree trignometric functions

gcc - fortran mailing list
On Sat, Mar 28, 2020 at 12:37 PM Steve Kargl
<[hidden email]> wrote:

>
> On Sat, Mar 28, 2020 at 08:10:38AM +0100, Tobias Burnus wrote:
> > Two remarks:
> >
> > * As the file trigd_lib.inc is shared between libgfortran
> >   and gcc/fortran, I wonder whether it should be placed
> >   under include/ (under the GCC toplevel directroy)
>
> The original name and location were chosen to match
> intrinsics/erfc_scaled_inc.c, which is included in
> intrinsics/erfc_scaled.c.  I think Fritz's re-use
> in simplification makes since given the ugly nested
> if-elseif-else constructs in the file.  As this is
> Fortran specific, I do not believe it should be moved
> up to a toplevel diretory.

I'm ambivalent about the location, though it is Fortran-specific.

> > * All included files need dependency; I do not quickly
> >   see whether that' the case.

Yes, I'd like to add a dependency for them as rebuilding after a
change is a pain. I've not yet decoded the Make-fu involved, but I
will be sure to add them for the final patch.

Thanks for the feedback, Tobias.

---
Fritz
Reply | Threaded
Open this post in threaded view
|

Re: PATCH -- Fix degree trignometric functions

gcc - fortran mailing list
On Mon, Mar 30, 2020 at 4:53 PM Tobias Burnus <[hidden email]> wrote:

>
> Hi Fritz,
>
> On 3/30/20 10:20 PM, Fritz Reese via Fortran wrote:
>
> >>> * All included files need dependency; I do not quickly
> >>>    see whether that' the case.
>
> If one looks at the build, the dependency is automatically
> obtained and tracked in …/.deps/*.Po in the build directory.
> Hence, no action needed.
>
> Cheers,
>
> Tobias

Awesome, thanks!

----
Fritz
Reply | Threaded
Open this post in threaded view
|

Re: PATCH -- Fix degree trignometric functions

gcc - fortran mailing list
In reply to this post by gcc - fortran mailing list
On Fri, Mar 27, 2020 at 7:36 PM Fritz Reese <[hidden email]> wrote:

>
> On Fri, Mar 6, 2020 at 6:18 PM Steve Kargl
> <[hidden email]> wrote:
> [...]
> > TL;DR version.
> >
> >   Fix the simplification and handling of the degree trigonometric functions.
> >   This includes fixing a number of ICEs.  See PR 93871.
>
> An updated version of the patch is attached. Regression tests (and new
> test cases) are pending.
>
Patch with regression tests attached, and rebased onto master
(013fca64fc...). The new and updated tests cover several issues
discovered and discussed in PR 93871 at
https://gcc.gnu.org/bugzilla/show_bug.cgi?id=93871 . I am
bootstrapping and testing now. If everything passes, does this look OK
to commit?

gcc/fortran/ChangeLog:

        PR fortran/93871
        * gfortran.h (GFC_ISYM_ACOSD, GFC_ISYM_ASIND, GFC_ISYM_ATAN2D,
        GFC_ISYM_ATAND, GFC_ISYM_COSD, GFC_ISYM_COTAND, GFC_ISYM_SIND,
        GFC_ISYM_TAND): New.
        * intrinsic.c (add_functions): Remove check for flag_dec_math.
        Give degree trig functions simplification and name resolution
        functions (e.g, gfc_simplify_atrigd () and gfc_resolve_atrigd ()).
        (do_simplify): Remove special casing of degree trig functions.
        * intrinsic.h (gfc_simplify_acosd, gfc_simplify_asind,
        gfc_simplify_atand, gfc_simplify_cosd, gfc_simplify_cotand,
        gfc_simplify_sind, gfc_simplify_tand, gfc_resolve_trigd2): Add new
        prototypes.
        (gfc_simplify_atrigd, gfc_simplify_trigd, gfc_resolve_cotan,
        resolve_atrigd): Remove prototypes of deleted functions.
        * iresolve.c (is_trig_resolved, copy_replace_function_shallow,
        gfc_resolve_cotan, get_radians, get_degrees, resolve_trig_call,
        gfc_resolve_atrigd, gfc_resolve_atan2d): Delete functions.
        (gfc_resolve_trigd, gfc_resolve_trigd2): Resolve to library functions.
        * simplify.c (rad2deg, deg2rad, gfc_simplify_acosd, gfc_simplify_asind,
        gfc_simplify_atand, gfc_simplify_atan2d, gfc_simplify_cosd,
        gfc_simplify_sind, gfc_simplify_tand, gfc_simplify_cotand): New
        functions.
        (gfc_simplify_atan2): Fix error message.
        (simplify_trig_call, gfc_simplify_trigd, gfc_simplify_atrigd,
        radians_f): Delete functions.
        * trans-intrinsic.c: Add LIB_FUNCTION decls for sind, cosd, tand.
        (rad2deg, gfc_conv_intrinsic_atrigd, gfc_conv_intrinsic_cotan,
        gfc_conv_intrinsic_cotand, gfc_conv_intrinsic_atan2d): New functions.
        (gfc_conv_intrinsic_function): Handle ACOSD, ASIND, ATAND, COTAN,
        COTAND, ATAN2D.
        * trigd_fe.inc: New file. Included by simplify.c to implement
        simplify_sind, simplify_cosd, simplify_tand with code common to the
        libgfortran implementation.

libgfortran/ChangeLog:

        PR fortran/93871
        * Makefile.am, Makefile.in: New make rule for intrinsics/trigd.c.
        * gfortran.map: New routines for {sind, cosd, tand}X{r4, r8, r10, r16}.
        * intrinsics/trigd.c, intrinsics/trigd_lib.inc, intrinsics/trigd.inc:
        New files. Defines native degree-valued trig functions.

gcc/testsuite/ChangeLog:

        PR fortran/93871
        * gfortran.dg/dec_math.f90: Extend coverage to real(10) and real(16).
        * gfortran.dg/dec_math_2.f90: New test.
        * gfortran.dg/dec_math_3.f90: Likewise.
        * gfortran.dg/dec_math_4.f90: Likewise.
        * gfortran.dg/dec_math_5.f90: Likewise.

trigd_v4.patch (138K) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: PATCH -- Fix degree trignometric functions

Andreas Schwab-2
FAIL: gfortran.dg/dec_math_5.f90   -O0  (test for excess errors)
Excess errors:
/opt/gcc/gcc-20200408/gcc/testsuite/gfortran.dg/dec_math_5.f90:132:9: Fatal Error: Cannot open module file 'ieee_arithmetic.mod' for reading at (1): No such file or directory
compilation terminated.

Andreas.

--
Andreas Schwab, [hidden email]
GPG Key fingerprint = 7578 EB47 D4E5 4D69 2510  2552 DF73 E780 A9DA AEC1
"And now for something completely different."
12