Bug in gmp/mpfr/mpc, never ending ctan computation

[Richard Biener]

On Mon, 19 Nov 2018, Torbjörn Granlund wrote:

> [moved from gmp-devel]
> 
> 
> Richard Biener <rguenther at suse.de> writes:
> 
>   For
> >
>   #include <complex.h>
> >
>   int main()
>   {
>     _Complex double x;
>     __real x = 3.09126495087690770626068115234375e+8;
>     __imag x = -4.045689747817175388336181640625e+8;
>     volatile _Complex double y = ctan (x);
>     return 0;
>   }
> 
> If we get a test case somewhat closer to GMP, then it is likely
> somebody will look at it.
> 
> I expect the maintainers of library called by gcc (mpc?) would be helped
> by seeing the actual call to their library.

OK, I can reproduce the issue with mpc 1.1.0, mpfr 3.1.4 and gmp 6.1.0
using

#include <mpc.h>
#include <mpfr.h>

int main()
{
  mpc_t m;
  mpc_init2 (m, 53);
  mpfr_set_str (mpc_realref (m), "3.09126495087690770626068115234375e+8", 
10, GMP_RNDN);
  mpfr_set_str (mpc_imagref (m), "-4.045689747817175388336181640625e+8", 
10, GMP_RNDN);
  mpfr_clear_flags ();
  mpc_tan (m, m, 0);
  return 0;
}

breaking at tan.c:189 and printing prec and then err I see

Breakpoint 1, mpc_tan (rop=0x7fffffffdd50, op=0x7fffffffdd50, rnd=0)
    at /space/rguenther/src/svn/trunk2/mpc/src/tan.c:189
189           prec += mpc_ceil_log2 (prec) + err;
$3 = 53
$4 = 7

Breakpoint 1, mpc_tan (rop=0x7fffffffdd50, op=0x7fffffffdd50, rnd=0)
    at /space/rguenther/src/svn/trunk2/mpc/src/tan.c:189
189           prec += mpc_ceil_log2 (prec) + err;
$5 = 66
$6 = 66

Breakpoint 1, mpc_tan (rop=0x7fffffffdd50, op=0x7fffffffdd50, rnd=0)
    at /space/rguenther/src/svn/trunk2/mpc/src/tan.c:189
189           prec += mpc_ceil_log2 (prec) + err;
$7 = 139
$8 = 139

...
Breakpoint 1, mpc_tan (rop=0x7fffffffdd50, op=0x7fffffffdd50, rnd=0)
    at /space/rguenther/src/svn/trunk2/mpc/src/tan.c:189
189           prec += mpc_ceil_log2 (prec) + err;
$29 = 303084
$30 = 303084
...

So somehow err ends up equal to prec and ever increasing?  Which means
we always end up in

      /* OP is no pure real nor pure imaginary, so in theory the real and
         imaginary parts of its tangent cannot be null. However due to
         rounding errors this might happen. Consider for example
         tan(1+14*I) = 1.26e-10 + 1.00*I. For small precision sin(op) and
         cos(op) differ only by a factor I, thus after mpc_div x = I and
         its real part is zero. */
      if (mpfr_zero_p (mpc_realref (x)) || mpfr_zero_p (mpc_imagref (x)))
        {
          err = prec; /* double precision */
          continue;
        }

As I questioned in the first followup I wonder if GCC can tell mpc
to not exceed 53bits of precision for the result?  Even with denormals
a precision of 303084 bits isn't possible?

Richard.

-- 
Richard Biener <rguenther at suse.de>
SUSE LINUX GmbH, GF: Felix Imendoerffer, Jane Smithard, Graham Norton, HRB 21284 (AG Nuernberg)