gsl_complex_tan with large negative imaginary number

[ohliumliu]

https://savannah.gnu.org/bugs/?47347

[ohliumliu]

1. 7141c7adb8100207ef01ff87a94d5179103fcebc This issue with tan is also applicable to tanh. In this commit, tanh was modified to have a separate branch to handle large real part. It is confirmed that this is necessary and also applicable to tan. gsl_complex_tan_2 implements this idea.
2. glibc implementation C99 has a built-in complex type. Its implement converts everything to exponential when the imagineary part is big. Is this a better implementation?
3. My testing In test_tan.c, I also checked and confirmed which intermediate value in gsl_complex_tan blows up. In general, inf is allowed in gsl, but it may not be smart enough to handle inf/inf. The idea in gsl_complex_tan_2 is to allow inf or 1/inf, but avoid inf/inf.
4. test cases: Built-in test cases of gsl doesn't seem to include this case.

TODO

1. Finished up test_tan.c and generate a patch for submission.
2. Add some test cases.
3. Study glibc implementation.

[ohliumliu]

• glibc implementation

double sinrx, cosrx;
      double den;
      const int t = (int) ((DBL_MAX_EXP - 1) * M_LN2 / 2);

      /* tan(x+iy) = (sin(2x) + i*sinh(2y))/(cos(2x) + cosh(2y))
     = (sin(x)*cos(x) + i*sinh(y)*cosh(y)/(cos(x)^2 + sinh(y)^2). */

      __sincos (__real__ x, &sinrx, &cosrx);

      if (fabs (__imag__ x) > t)
    {
      /* Avoid intermediate overflow when the real part of the
         result may be subnormal.  Ignoring negligible terms, the
         imaginary part is +/- 1, the real part is
         sin(x)*cos(x)/sinh(y)^2 = 4*sin(x)*cos(x)/exp(2y).  */
      double exp_2t = __ieee754_exp (2 * t);

      __imag__ res = __copysign (1.0, __imag__ x);
      __real__ res = 4 * sinrx * cosrx;
      __imag__ x = fabs (__imag__ x);
      __imag__ x -= t;
      __real__ res /= exp_2t;
      if (__imag__ x > t)
        {
          /* Underflow (original imaginary part of x has absolute
         value > 2t).  */
          __real__ res /= exp_2t;
        }
      else
        __real__ res /= __ieee754_exp (2 * __imag__ x);
    }
      else
    {
      double sinhix = __ieee754_sinh (__imag__ x);
      double coshix = __ieee754_cosh (__imag__ x);

      den = cosrx * cosrx + sinhix * sinhix;
      __real__ res = sinrx * cosrx / den;
      __imag__ res = sinhix * coshix / den;
    }
    }

The idea is the following.

1. In the following, we are talking about x+iy as input. In the code above, x is a complex input itself.
2. The starting point is tan(x+iy) = (sin(2x) + i*sinh(2y))/(cos(2x) + cosh(2y)) = (sin(x)*cos(x) + i*sinh(y)*cosh(y)/(cos(x)^2 + sinh(y)^2).
3. The overflow comes from sinh(y)^2, in particular exp(2y) or exp(-2y).
4. The overflow threshold for fabs(y) is t = (int) ((DBL_MAX_EXP - 1) * M_LN2 / 2). Specifically,

exp(2t) <= 2^(DBL_MAX_EXP-1)
t <= (DBL_MAX_EXP-1) * M_LN2/2

where DBL_MAX_EXP is

More precisely, this is the maximum positive integer such that value FLT_RADIX raised to this power minus 1 can be represented as a floating point number of type DBL and FLT_RADIX is 2, and M_LN2 is the natural log of 2.

When y is bigger than the threshold. The imaginary part of the result is +/-i based on the sign of y. The real part is 4*sin(x)*cos(x)/exp(2*fabs(y)). We need to process the denominator. Generally speaking, we need to writefabs(y)as the sum of a number oftandy mod t`.

In the code above, y is first replaced by fabs(y) and then by fabs(y) - t. If fabs(y) - t > t, exp(2t)*exp(2t) is used instead of exp(2*fabs(y)); if fabs(y) - t < t, exp(2t)*exp(2(fabs(y) - t)) = exp(2*fabs(y)) is used in the denominator. In the first case, it's an approximation, and in the second case, it is accurate in math. This is a simplified version of the idea in last paragraph.