https://gcc.gnu.org/bugzilla/show_bug.cgi?id=93871
--- Comment #33 from Thomas Henlich <thenlich at gcc dot gnu.org> ---
Going back one step, I wonder if it would be good enough to perform a correctly
rounded conversion from degrees to radians, and just use sin() as is.
...
f = sgn * sin(deg2rad(arg))
end function fcn
! Compute correctly rounded value of x * pi / 180 for x = 0 ... 90
function deg2rad(x) result(y)
real(sp) y
real(sp), intent(in) :: x
real(sp), parameter :: a = atan(1._4) / 45
real(sp), parameter :: b = atan(1._16) / 45 - a
y = a * x + b * x
end function deg2rad
(The exact values of a and b still need some work)