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

[Vincent Lefevre]

On 2018-11-19 11:14:49 +0100, Richard Biener wrote:
[...]
> 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;
>         }

A zero value suggests an underflow somewhere, which implies that...

> 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?

increasing the precision will generally not solve the underflow
issue, which is due to the limited exponent range. The algorithm
should be modified to do some form of scaling (when possible).

MPFR has internally some very partial support for "unbounded floats"
(UBF), i.e. where the exponent is a mpz_t. They are currently used
for mpfr_fma and mpfr_fmma.

-- 
Vincent Lefèvre <vincent at vinc17.net> - Web: <https://www.vinc17.net/>
100% accessible validated (X)HTML - Blog: <https://www.vinc17.net/blog/>
Work: CR INRIA - computer arithmetic / AriC project (LIP, ENS-Lyon)