[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 #32 from Steve Kargl <sgk at troutmask dot apl.washington.edu> ---
On Thu, Feb 27, 2020 at 05:53:38PM +0000, thenlich at gcc dot gnu.org wrote:
> https://gcc.gnu.org/bugzilla/show_bug.cgi?id=93871
> 
> --- Comment #31 from Thomas Henlich <thenlich at gcc dot gnu.org> ---
> I wonder, if some "correct" rounding could further increase
> accuracy: We know the sign and "real" magnitude of the
> difference deg2rad-π/180 and can round the result of sin()
> accordingly.

I think we can exploit a table-driven method.  Once x is 
folded into [0,90], x is split by x = n + dx where
n = 0, 1, ..., 90 and dx is in [0.,1.).

sin(x) = sin(n+d) = sin(n) * cos(d) + cos(n) * sin(d)

The table contains a high,low decomposition of {sin(1),cos(1)},
{sin(3),cos(3)}, ..., {sin(90),cos(90)}.  Denote the decomposition
with shi(n), slo(n), chi(n), clo(n), and both C(d) = cos(d) and
S(d) * sin(d) are minimax polynomials. So, we have

if (x == n)  /* Integral value of x.  */
   sin(x) = shi(n) + slo(n)   ! Is addition needed.
else
   sin(x) = shi(n)*C(d) + slo(n)*C(d) + chi(n)*S(d) + clo(n)*S(d) 

The else-branch summation can be done with Kahan summation.