[Bug libfortran/93871] COTAN is slow for complex types

[sgk at troutmask dot apl.washington.edu]

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=93871

--- Comment #28 from Steve Kargl <sgk at troutmask dot apl.washington.edu> ---
On Wed, Feb 26, 2020 at 04:02:23PM +0000, sgk at troutmask dot
apl.washington.edu wrote:
> 
> > This is a best effort to still be able to use the standard library functions
> > but also get an increased accuracy expected from the degree functions.
> > 
> > Thus:
> > 1. Calculate cosd(x) = sind(90 - x)
> > 2. Calculate cotand(x) = tand(90 - x)
> > 3. Reduce range of sind() argument from (0...360) further to x-360 if it is
> > above 270, and to 180-x if it is above 90
> > 
> 

Exhaustive testing in the domain [2**(-8),720] for REAL(4) gives

troutmask:sgk[347] gfcx -o z -O2 -fdec e.f90 && ./z
Total values tested: 146014209
               ulp > 1.5   max ulp   x at max ulp 
WIP Patch[1]:   6558266   4331134.9   359.999969    
fcn below[2]:       717   1.6212178   1.79070497    

[1] Work in progress has implemented range reduction in the
    interval [0,360], but does not fold the reduced range
    as is done in fcn(x).  The reference solution in the
    ulp computation is fcn(x) converted to double precision.

[2] Only 9 exceed a max ulp of 1.6.

      ! Compute SIND(x) where x is in degrees.
      function fcn(x) result(f)
         real(sp) f
         real(sp), intent(in) :: x
         integer sgn
         real(sp) arg
         real(sp), parameter :: deg2rad = atan(1._sp) / 45
         !
         ! Required for sin(-x) = - sin(x)
         !
         sgn = 1
         if (x < 0) sgn = -1
         arg = abs(x)

         ! 0 <= arg < 360
         arg = modulo(arg, 360._sp)

         ! Fold  into [0,90] range
         if (arg <= 180) then
            if (arg > 90) arg = 180 - arg
         else
            sgn = -sgn
            arg = arg - 180
            if (arg > 90) arg = 180 - arg
         end if
         if (arg == 180) then
            f = sign(0._sp, real(sgn,sp))
         else
            f = sgn * sin(arg * deg2rad)
         end if
      end function fcn

Probably need to punt the above into libgfortran.  Note,
for x < 2**(-8) we can do sin(x) ~ x