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
