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

[Vincent Lefevre]

On 2018-11-19 15:26:00 +0100, Richard Biener wrote:
> On Mon, 19 Nov 2018, Vincent Lefevre wrote:
> > AFAIK, GCC already changes the exponent range to the binary64 one
> > (the lack of subnormals is not an issue, because in case of underflow,
> > GCC ignores the computations done with MPFR).
> 
> GCC does ignore results in case of overflow or underflow but it
> doesn't currently limit mpfrs exponent range.  Or rather most
> languages don't - I do see mpfr_set_e{min,max} calls being
> done in the Fortran frontend only, but only temporarily and
> seemingly only to call mpfr_subnormalize.

Thanks for the clarification. FYI, mpfr_subnormalize was added
for gfortran.

> > But for its internal computations, I assume that mpc extends the
> > exponent range (if it didn't you would have the same kind of problem,
> > but much more often). Here it seems that the issue is that the maximum
> > supported exponent range is not large enough for the chosen algorithm
> > (which may not be the right one for this kind of input).
> 
> It happens that limiting the range makes the CPOW computation fast.

Yes, because for some functions, underflow or overflow can be detected
before the actual computation.

-- 
Vincent Lefèvre  - Web: 
100% accessible validated (X)HTML - Blog: 
Work: CR INRIA - computer arithmetic / AriC project (LIP, ENS-Lyon)
___
gmp-bugs mailing list
gmp-bugs@gmplib.org
https://gmplib.org/mailman/listinfo/gmp-bugs

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

[Dennis Clarke]

On 11/19/18 5:40 AM, Torbjörn Granlund wrote:

Richard Biener  writes:

   > __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 
   #include 

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


using that test I see 100% cpu and slowly growing RSS with a lot of time 
being spent inside mpn_fft_fft and then mpn_fft_mul_2exp_modF :


Program received signal SIGINT, Interrupt.
0x0008002dcd04 in __gmpn_copyi () from 
/home/dclarke/local/lib/libgmp.so.10

(gdb) where
#0  0x0008002dcd04 in __gmpn_copyi ()
   from /home/dclarke/local/lib/libgmp.so.10
#1  0x00080029a021 in mpn_fft_mul_2exp_modF (r=0x800e47610, 
a=0x80150bcb0,

d=1600, n=40) at mul_fft.c:267
#2  0x00080029a3f0 in mpn_fft_fft (Ap=0x800e8e130, K=32, 
ll=0x800e00098,

omega=160, n=40, inc=16, tp=0x800e47610) at mul_fft.c:391
#3  0x00080029a351 in mpn_fft_fft (Ap=0x800e8dc30, K=64, 
ll=0x800e000a0,

omega=80, n=40, inc=8, tp=0x800e47610) at mul_fft.c:383
#4  0x00080029a3b2 in mpn_fft_fft (Ap=0x800e8dc10, K=128, 
ll=0x800e000a8,

omega=40, n=40, inc=4, tp=0x800e47610) at mul_fft.c:384
#5  0x00080029a351 in mpn_fft_fft (Ap=0x800e8dc10, K=256, 
ll=0x800e000b0,

omega=20, n=40, inc=2, tp=0x800e47610) at mul_fft.c:383
#6  0x00080029a351 in mpn_fft_fft (Ap=0x800e8dc10, K=512, 
ll=0x800e000b8,

omega=10, n=40, inc=1, tp=0x800e47610) at mul_fft.c:383
#7  0x0008002995ff in mpn_mul_fft_internal (op=0x8014d4150, pl=9728, 
k=9,
Ap=0x800e8dc10, Bp=0x800e8b410, A=0x8014fd610, B=0x801060950, 
nprime=40,

l=19, Mp=5, fft_l=0x800e00070, T=0x800e47610, sqr=1) at mul_fft.c:738
#8  0x000800298ef8 in __gmpn_mul_fft (op=0x8014d4150, pl=9728,
n=0x801336980, nl=9475, m=0x801336980, ml=9475, k=9) at mul_fft.c:904
#9  0x0008002cc6d6 in __gmpn_sqrmod_bnm1 (rp=0x8014a8500, rn=19456,
ap=0x801336980, an=9475, tp=0x8014d4150) at sqrmod_bnm1.c:194
#10 0x00080029c568 in __gmpn_nussbaumer_mul (pp=0x8014a8500,
ap=0x801336980, an=9475, bp=0x801336980, bn=9475) at 
nussbaumer_mul.c:61

--Type  for more, q to quit, c to continue without paging--
#11 0x00080029b80f in __gmpn_sqr (p=0x8014a8500, a=0x801336980, n=9475)
at sqr.c:97
#12 0x000800371726 in mpfr_sqrhigh_n (rp=0x801483500, np=0x801324180,
n=18947) at mulders.c:180
#13 0x00080031c0ab in mpfr_mul (a=0x7fffde50, b=0x7fffde50,
c=0x7fffde50, rnd_mode=MPFR_RNDD) at mul.c:968
#14 0x00080036b2c0 in mpfr_sqr (a=0x7fffde50, b=0x7fffde50,
rnd_mode=MPFR_RNDD) at sqr.c:561
#15 0x000800312bd0 in mpfr_exp_3 (y=0x7fffe350, x=0x7fffe380,
rnd_mode=MPFR_RNDD) at exp3.c:232
#16 0x00080031487f in mpfr_exp (y=0x7fffe350, x=0x7fffe380,
rnd_mode=MPFR_RNDD) at exp.c:176
#17 0x00080033a525 in mpfr_sinh_cosh (sh=0x7fffe550,
ch=0x7fffe530, xt=0x7fffe880, rnd_mode=MPFR_RNDN)
at sinh_cosh.c:108
#18 0x0008003bf37b in mpc_sin_cos (rop_sin=0x7fffe7e0,
rop_cos=0x7fffe7a0, op=0x7fffe860, rnd_sin=17, rnd_cos=17)
at sin_cos.c:390
#19 0x0008003c3756 in mpc_tan (rop=0x7fffe860, op=0x7fffe860,
rnd=0) at tan.c:196
#20 0x002013c0 in main () at test_case.c:19
(gdb) quit
A debugging session is active.

Inferior 1 [process 58829] will be killed.

Quit anyway? (y or n) y
titan$

This is gmp 6.1.2 and mpfr-4.0.1-p13 and mpc-1.1.0 on FreeBSD 12.0RC1
 and I also see similar behavior on Solaris SPARC.



I simplified your test case using 1 as the real part and an integer for
the imaginary part.  With the real part set to 1, the computation
finishes in about 10 seconds, with every doubling of it the runtime
almost triples.  (If my observation generalises to an extended real part
range, the compilation for your test case should terminate in about
10*3^log2(4.045689747817175388336181640625e+8/1)/3600/24/365 years.)



Well the test case provided merely packs up and goes away with no
indication that it will return anytime soon.

I am just going to watch this thread closely.

Dennis

___
gmp-bugs mailing list
gmp-bugs@gmplib.org
https://gmplib.org/mailman/listinfo/gmp-bugs

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

[Richard Biener]

On Mon, 19 Nov 2018, Vincent Lefevre wrote:

> On 2018-11-19 14:17:13 +0100, Niels Möller wrote:
> > You might want to add an entry for "Exponent range" to the concept index
> > in the manual.
> 
> The concept index already has "Exponent", and the corresponding
> paragraph is:
> 
>An “exponent” is a component of a regular floating-point number.  Its
> C data type is ‘mpfr_exp_t’.  Valid exponents are restricted to a subset
> of this type, and the exponent range can be changed globally as
> described in *note Exception Related Functions::.  Special values do not
> have an exponent.
> 
> I think that this is quite explicit.
> 
> > I guess gcc should set things up similarly to the example on that page,
> > "This is an example of how to emulate binary double IEEE 754
> > arithmetic", with the caveat "This emulates a double IEEE 754 arithmetic
> > with correct rounding in the subnormal range, which may not be the case
> > for your hardware".
> 
> AFAIK, GCC already changes the exponent range to the binary64 one
> (the lack of subnormals is not an issue, because in case of underflow,
> GCC ignores the computations done with MPFR).

GCC does ignore results in case of overflow or underflow but it
doesn't currently limit mpfrs exponent range.  Or rather most
languages don't - I do see mpfr_set_e{min,max} calls being
done in the Fortran frontend only, but only temporarily and
seemingly only to call mpfr_subnormalize.

> But for its internal computations, I assume that mpc extends the
> exponent range (if it didn't you would have the same kind of problem,
> but much more often). Here it seems that the issue is that the maximum
> supported exponent range is not large enough for the chosen algorithm
> (which may not be the right one for this kind of input).

It happens that limiting the range makes the CPOW computation fast.

Richard.

-- 
Richard Biener 
SUSE LINUX GmbH, GF: Felix Imendoerffer, Jane Smithard, Graham Norton, HRB 
21284 (AG Nuernberg)
___
gmp-bugs mailing list
gmp-bugs@gmplib.org
https://gmplib.org/mailman/listinfo/gmp-bugs

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

[Richard Biener]

On Mon, 19 Nov 2018, Niels Möller wrote:

> Vincent Lefevre  writes:
> 
> > MPFR has a bounded (user configurable) exponent range, with underflow
> > and overflow exceptions for out-of-range results, like in IEEE 754.
> 
> Good! I didn't find that in the manual until I knew the phrase to search
> for. It's in
> https://www.mpfr.org/mpfr-current/mpfr.html#Exception-Related-Functions,
> and relevant functions are mpfr_set_emin, mpfr_set_emax and mpfr_subnormalize.
> 
> You might want to add an entry for "Exponent range" to the concept index
> in the manual.

Indeed I didn't find these functions either ...

> I guess gcc should set things up similarly to the example on that page,
> "This is an example of how to emulate binary double IEEE 754
> arithmetic", with the caveat "This emulates a double IEEE 754 arithmetic
> with correct rounding in the subnormal range, which may not be the case
> for your hardware".

OK, GCC doesn't use those, for the testcase if I set emax/emin to the
range from 128bit IEEE then the folding finishes quickly.

It looks like emin/emax is global and not per mpfr_t so for GCC we
have to jump through some additional hoops.  I'll see to adjust
GCC to set the limits based on the maximum range ever needed
(largest FP mode supported by the target) and look for fallout.

Richard.

> Regards,
> /Niels
> 
> 

-- 
Richard Biener 
SUSE LINUX GmbH, GF: Felix Imendoerffer, Jane Smithard, Graham Norton, HRB 
21284 (AG Nuernberg)
___
gmp-bugs mailing list
gmp-bugs@gmplib.org
https://gmplib.org/mailman/listinfo/gmp-bugs

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

[Vincent Lefevre]

On 2018-11-19 12:49:12 +0100, Richard Biener wrote:
> I see.  So I guess one issue might be that mpfr numbers have a
> precision (mantissa bits?) but no restriction on the exponent?

MPFR has a bounded (user configurable) exponent range, with underflow
and overflow exceptions for out-of-range results, like in IEEE 754.

-- 
Vincent Lefèvre  - Web: 
100% accessible validated (X)HTML - Blog: 
Work: CR INRIA - computer arithmetic / AriC project (LIP, ENS-Lyon)
___
gmp-bugs mailing list
gmp-bugs@gmplib.org
https://gmplib.org/mailman/listinfo/gmp-bugs

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

[Richard Biener]

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

> Richard Biener  writes:
> 
>   :/  But then with GCC we want deterministic behavior -- ISL
>   added some "computation budget" and error reporting if that was
>   taken up.  Maybe gmp/mpfr/mpc could consider adding such a thing
>   (reporting budget overruns via a signal could be acceptable if
>   the API churn would be too bad for another way).
> 
> I'll leave this up to the mpc maintainers.
> 
>   > I simplified your test case using 1 as the real part and an integer for
>   > the imaginary part.  With the real part set to 1, the computation
>   > finishes in about 10 seconds, with every doubling of it the runtime
>   > almost triples.  (If my observation generalises to an extended real part
>   > range, the compilation for your test case should terminate in about
>   > 10*3^log2(4.045689747817175388336181640625e+8/1)/3600/24/365 years.)
> 
>   Ick.  I wonder if the actual ctan result is infinite or we just
>   fail to compute it in an efficient manner...
> 
> The result should be extremely small, so returning Inf would be
> quite inaccurate.  :-)

I see.  So I guess one issue might be that mpfr numbers have a
precision (mantissa bits?) but no restriction on the exponent?
That is, for ctan eventually zero _is_ a validly rounded result
and the computation should simply terminate...

Richard.

-- 
Richard Biener 
SUSE LINUX GmbH, GF: Felix Imendoerffer, Jane Smithard, Graham Norton, HRB 
21284 (AG Nuernberg)
___
gmp-bugs mailing list
gmp-bugs@gmplib.org
https://gmplib.org/mailman/listinfo/gmp-bugs

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

[Torbjörn Granlund]

Richard Biener  writes:

  :/  But then with GCC we want deterministic behavior -- ISL
  added some "computation budget" and error reporting if that was
  taken up.  Maybe gmp/mpfr/mpc could consider adding such a thing
  (reporting budget overruns via a signal could be acceptable if
  the API churn would be too bad for another way).

I'll leave this up to the mpc maintainers.

  > I simplified your test case using 1 as the real part and an integer for
  > the imaginary part.  With the real part set to 1, the computation
  > finishes in about 10 seconds, with every doubling of it the runtime
  > almost triples.  (If my observation generalises to an extended real part
  > range, the compilation for your test case should terminate in about
  > 10*3^log2(4.045689747817175388336181640625e+8/1)/3600/24/365 years.)

  Ick.  I wonder if the actual ctan result is infinite or we just
  fail to compute it in an efficient manner...

The result should be extremely small, so returning Inf would be
quite inaccurate.  :-)

-- 
Torbjörn
Please encrypt, key id 0xC8601622
___
gmp-bugs mailing list
gmp-bugs@gmplib.org
https://gmplib.org/mailman/listinfo/gmp-bugs

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

[Torbjörn Granlund]

t...@gmplib.org (Torbjörn Granlund) writes:

  I simplified your test case using 1 as the real part and an integer for
  the imaginary part.  With the real part set to 1, the computation
  finishes in about 10 seconds, with every doubling of it the runtime
  almost triples.  (If my observation generalises to an extended real part
  range, the compilation for your test case should terminate in about
  10*3^log2(4.045689747817175388336181640625e+8/1)/3600/24/365 years.)

I meant to write:

I simplified your test case using 1 as the real part and an integer for
the imaginary part.  With the imaginary part set to 1, the computation
finishes in about 10 seconds, with every doubling of it the runtime almost
triples.  (If my observation generalises to an extended imaginary part
range, the compilation for your test case should terminate in about
10*3^log2(4.045689747817175388336181640625e+8/1)/3600/24/365 years.)

-- 
Torbjörn
Please encrypt, key id 0xC8601622
___
gmp-bugs mailing list
gmp-bugs@gmplib.org
https://gmplib.org/mailman/listinfo/gmp-bugs