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

[Niels Möller]

Richard Biener <rguenther at suse.de> writes:

> #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;
> }

What's the expected value of y? I'm not that familiar with complex trig
functions.  I tried evaluating it in pari/gp, which gives

? x = 3.09126495087690770626068115234375e+8 + I*-4.045689747817175388336181640625e+8
%2 = 309126495.08769077062606811523437500000 - 404568974.78171753883361816406250000000*I
? tan(x)
  ***   at top-level: tan(x)
  ***                 ^------
  *** tan: the PARI stack overflows !
  current stack size: 8000000 (7.629 Mbytes)
  [hint] set 'parisizemax' to a non-zero value in your GPRC

which makes me think the value is somehow badly behaved.

> I see GCC never finishing constant folding of ctan () with
> backtraces like
>
> #0  0x00007ffff735666b in ?? () from /usr/lib64/libgmp.so.10
> #1  0x00007ffff7356540 in ?? () from /usr/lib64/libgmp.so.10
> #2  0x00007ffff7356540 in ?? () from /usr/lib64/libgmp.so.10
> #3  0x00007ffff7356540 in ?? () from /usr/lib64/libgmp.so.10
> #4  0x00007ffff7357157 in ?? () from /usr/lib64/libgmp.so.10
> #5  0x00007ffff7357ecf in __gmpn_mul_fft () from /usr/lib64/libgmp.so.10
> #6  0x00007ffff737a7cf in __gmpn_mulmod_bnm1 () from 
> /usr/lib64/libgmp.so.10
> #7  0x00007ffff737a553 in __gmpn_mulmod_bnm1 () from 
> /usr/lib64/libgmp.so.10
> #8  0x00007ffff737ab39 in __gmpn_mulmod_bnm1 () from 
> /usr/lib64/libgmp.so.10
> #9  0x00007ffff7358654 in __gmpn_nussbaumer_mul () from 
> /usr/lib64/libgmp.so.10

That we end up in gmp's fft multiply code also indicates a huge number.
I don't think code trying to evaluate a number to IEEE double precision
should ever do that.

> #10 0x00007ffff735829e in __gmpn_mul_n () from /usr/lib64/libgmp.so.10
> #11 0x00007ffff737e5e9 in __gmpn_preinv_mu_div_qr ()
>    from /usr/lib64/libgmp.so.10
> #12 0x00007ffff737ebf9 in ?? () from /usr/lib64/libgmp.so.10
> #13 0x00007ffff7362638 in __gmpn_tdiv_qr () from /usr/lib64/libgmp.so.10
> #14 0x00007ffff7353e24 in __gmpn_divrem () from /usr/lib64/libgmp.so.10
> #15 0x00007ffff735b64c in ?? () from /usr/lib64/libgmp.so.10
> #16 0x00007ffff735baf5 in __gmpn_sqrtrem () from /usr/lib64/libgmp.so.10
> #17 0x00007ffff7344102 in __gmpz_sqrt () from /usr/lib64/libgmp.so.10
> #18 0x00007ffff75d4c0c in ?? () from /usr/lib64/libmpfr.so.4
> #19 0x00007ffff75d50a1 in ?? () from /usr/lib64/libmpfr.so.4
> #20 0x00007ffff75d563d in mpfr_sincos_fast () from /usr/lib64/libmpfr.so.4

My best guess is that mpfr should detect the difficulty, fail and return
an inf or NaN as appropriate.

I'd imagine it's generally a bit tricky to accurately compute
trigonometric functions for large inputs, since one needs an accurate
subtraction of the right multiple of pi.

> #21 0x00007ffff75d5c8c in mpfr_sin_cos () from /usr/lib64/libmpfr.so.4
> #22 0x00007ffff7810b79 in mpc_sin_cos () from /usr/lib64/libmpc.so.3
> #23 0x00007ffff78138ab in mpc_tan () from /usr/lib64/libmpc.so.3
> #24 0x0000000000e835cf in do_mpc_arg1 (result_real=0x7fffffffcee0, 
>
> So I'm not sure which library to report the issue on (and the MPC
> bugtracker on gforge seems unused).  So I hope relevant devs follow
> the GMP list and can identify the problematic library.
>
> This has been reported as GCC bug here: 
> https://gcc.gnu.org/bugzilla/show_bug.cgi?id=88074

Regards,
/Niels

-- 
Niels Möller. PGP-encrypted email is preferred. Keyid 368C6677.
Internet email is subject to wholesale government surveillance.