
12

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


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 degreetrig 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


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 callingconvention 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_64pclinuxgnu
$ ./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/gccgit/gcc/testsuite/gfortran.dg/sind_intrinsic.f90:6
#2 0x40080f in main
at /data/foreese/src/gccgit/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/gccgit/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 $rip24,+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 floatingpoint
parameter directly byvalue in the SSE register %xmm0, while the
library routine expects a pointer to the value in %rdi. This is the IR
with fdumptreeoriginal:
$ 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 fdumptreeoriginal
$ ./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 middleend translation yet,
so it's not immediately clear to me why the calling code is passing
byvalue instead of byreference 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 byreference. The
biggest difference I see is the LIB_FUNCTION() lines in the
gfc_intrinsic_map in transintrinsic.c: why are these present for the
trigd functions and not for other libgfortran functions?
$ grep C4 'LIB_FUNCTION (SIND' transintrinsic.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


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 callingconvention 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_64pclinuxgnu
> $ ./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/gccgit/gcc/testsuite/gfortran.dg/sind_intrinsic.f90:6
> #2 0x40080f in main
> at /data/foreese/src/gccgit/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/gccgit/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 $rip24,+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 floatingpoint
> parameter directly byvalue in the SSE register %xmm0, while the
> library routine expects a pointer to the value in %rdi. This is the IR
> with fdumptreeoriginal:
>
> $ 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 fdumptreeoriginal
> $ ./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 middleend translation yet,
> so it's not immediately clear to me why the calling code is passing
> byvalue instead of byreference 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 byreference. The
> biggest difference I see is the LIB_FUNCTION() lines in the
> gfc_intrinsic_map in transintrinsic.c: why are these present for the
> trigd functions and not for other libgfortran functions?
>
> $ grep C4 'LIB_FUNCTION (SIND' transintrinsic.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 fdumptreeoriginal 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 20200316 17:44:02.756817000 0700
+++ trigd_inc.c 20200316 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.e30f;
@@ 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) = 1x*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


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 callingconvention mismatch between the frontend and
> > the library
[...]
> I may have misread the fdumptreeoriginal 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 PUREity?
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.e30f;
#elif TINY
static const volatile WTYPE tiny = 1.e300f;
#else
static const volatile WTYPE tiny = 0x1p10000L;
#endif
[...]
/* No spurious underflows!. In radians, cos(x) = 1x*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.e300f is a single precision literal and will always be zero
 this is probably meant to be 1.e300 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.p27L
#ifdef HAVE_GFC_REAL_16
[...]
#ifdef GFC_REAL_16_IS_FLOAT128 /* libquadmath. */
#define SMALL 3.4e17Q
[...]
#else
#define SMALL 3.4e17L
[...]
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
passbyvaluepatched 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


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 fdumptreeoriginal 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 PUREity?
Like you, I struggle with the trans* files. If you look in
transintrinsic.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 123125, 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.e30f;
> #elif TINY
> static const volatile WTYPE tiny = 1.e300f;
> #else
> static const volatile WTYPE tiny = 0x1p10000L;
> #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 constantfolding the
subtraction, which is needed to raise the FE_INEXACT exception.
Unfortunately, 1.e30f 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.e30f;
#else
static const volatile WTYPE tiny = 1.e300f;
#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) = 1x*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.e300f is a single precision literal and will always be zero
>  this is probably meant to be 1.e300 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.p27L
> #ifdef HAVE_GFC_REAL_16
> [...]
> #ifdef GFC_REAL_16_IS_FLOAT128 /* libquadmath. */
> #define SMALL 3.4e17Q
> [...]
> #else
> #define SMALL 3.4e17L
> [...]
>
> 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
> passbyvaluepatched 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.pN where N > P is chosen
appropriately for each type. This only works for ax < 0x1.pP where
P is precision.

steve


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
> transintrinsic.c, you'll find gfc_conv_intrinsic_function(),
> [...]
Ah, thanks!
> [...] This can be changed to
>
> #ifdef TINYF
> static const volatile WTYPE tiny = 1.e30f;
> #else
> static const volatile WTYPE tiny = 1.e300f;
> #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 glibc2.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.pN 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


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.e30f;
> > #else
> > static const volatile WTYPE tiny = 1.e300f;
Note, I kept the typo! 1.e300f should not have
the f suffice. It should be 1.e300.
> > 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 microoptimization.
> > [...] 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


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.e30f;
> > > #else
> > > static const volatile WTYPE tiny = 1.e300f;
>
> Note, I kept the typo! 1.e300f should not have
> the f suffice. It should be 1.e300.
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 microoptimization.
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


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*(1x**2/6).
Thanks for carry this to a finish line.

Steve


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
radianvalued as well as degreevalued intrinsics.
I don't think we need to sweat over this too much, since glibc's
radiantrig 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


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
> radianvalued as well as degreevalued intrinsics.
>
> I don't think we need to sweat over this too much, since glibc's
> radiantrig 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 inlined it with gfc_conv_intrinsic_cotand(). We
can change the inlining to generate 1/tand(x). Are you up for
changing transintrinsic.c or would you like me to do the change?

Steve


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 inlined it with gfc_conv_intrinsic_cotand(). We
> can change the inlining to generate 1/tand(x). Are you up for
> changing transintrinsic.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


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 frontend for simplification. The GMP/MPFR functions are
preprocessed 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.
* transintrinsic.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 degreevalued trig functions.


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 frontend for simplification. The GMP/MPFR functions are
> preprocessed 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.
> * transintrinsic.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 degreevalued trig functions.


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 reuse
in simplification makes since given the ugly nested
ifelseifelse 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


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 reuse
> in simplification makes since given the ugly nested
> ifelseifelse 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 Fortranspecific.
> > * 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 Makefu involved, but I
will be sure to add them for the final patch.
Thanks for the feedback, Tobias.

Fritz


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


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.
* transintrinsic.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 degreevalued 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.


FAIL: gfortran.dg/dec_math_5.f90 O0 (test for excess errors)
Excess errors:
/opt/gcc/gcc20200408/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
