On Mon, 19 Nov 2018, Torbjörn Granlund wrote:
> [moved from gmp-devel]
>
>
> Richard Biener <[email protected]> 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 <[email protected]>
SUSE LINUX GmbH, GF: Felix Imendoerffer, Jane Smithard, Graham Norton, HRB
21284 (AG Nuernberg)
_______________________________________________
gmp-bugs mailing list
[email protected]
https://gmplib.org/mailman/listinfo/gmp-bugs
