# Re: PATCH -- Fix degree trignometric functions Classic List Threaded 21 messages 12
Open this post in threaded view
|

## Re: PATCH -- Fix degree trignometric functions

 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
Open this post in threaded view
|

## Re: PATCH -- Fix degree trignometric functions

 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
Open this post in threaded view
|

## Re: PATCH -- Fix degree trignometric functions

 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 :    sub    \$0x220,%rsp    0x000000000040075d :    mov    0x200904(%rip),%rax  # 0x601068    0x0000000000400764 :    movq   %rax,%xmm0    0x0000000000400769 :    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
Open this post in threaded view
|

## Re: PATCH -- Fix degree trignometric functions

Open this post in threaded view
|

## Re: PATCH -- Fix degree trignometric functions

 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
Open this post in threaded view
|

## Re: PATCH -- Fix degree trignometric functions

Open this post in threaded view
|

## Re: PATCH -- Fix degree trignometric functions

 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
Open this post in threaded view
|

## Re: PATCH -- Fix degree trignometric functions

 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
Open this post in threaded view
|

## Re: PATCH -- Fix degree trignometric functions

 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
Open this post in threaded view
|

## Re: PATCH -- Fix degree trignometric functions

 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
Open this post in threaded view
|

## Re: PATCH -- Fix degree trignometric functions

 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
Open this post in threaded view
|

## Re: PATCH -- Fix degree trignometric functions

 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
Open this post in threaded view
|

## Re: PATCH -- Fix degree trignometric functions

 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
Open this post in threaded view
|

## Re: PATCH -- Fix degree trignometric functions

 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
Open this post in threaded view
|

## Re: PATCH -- Fix degree trignometric 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 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.
Open this post in threaded view
|

## Re: PATCH -- Fix degree trignometric 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 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
Open this post in threaded view
|

## Re: PATCH -- Fix degree trignometric functions

 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
Open this post in threaded view
|

## Re: PATCH -- Fix degree trignometric functions

 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
Open this post in threaded view
|

## Re: PATCH -- Fix degree trignometric functions

 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