Issue with large factorial

[wjamesjr]

I know that I am doing something wrong but just can't figure what. The largest factorial that I can get converted to a mpfr is 44787927. Above this number I get inf as a result. I am on gmpy2 2.15 and python 3.11.7 on an Apple silicon Mac M3.

if name == "main": n = 44787927 gmpy2.get_context().emin = gmpy2.get_emin_min() gmpy2.get_context().emax = gmpy2.get_emax_max() print("emin", gmpy2.get_context().emin) print("emax", gmpy2.get_context().emax) print("max precision", gmpy2.get_max_precision()) gmpy2.get_context().precision = 100 print("precision set to", gmpy2.get_context().precision) factn = gmpy2.factorial(n) print(factn) n = 44787928 factn = gmpy2.factorial(n) print(factn) exit(0)

I get the following result running the above code.

emin -4611686018427387903 emax 4611686018427387903 max precision 9223372036854775551 precision set to 100 1.9481086289383819889097566536014e+323228493 inf

[casevh]

gmpy2.gamma() has a similar behavior. I vaguely recollect another report of this issue and it was a limitation of the MPFR library. I haven't found the details but I'll keep looking.

[casevh]

I found this comment in MPFR's gamma.c:

 if (u < 44787929UL && bits_fac (u - 1) <= p + (rnd_mode == MPFR_RNDN))
    /* bits_fac: lower bound on the number of bits of m,
       where gamma(x) = (u-1)! = m*2^e with m odd. */
    return mpfr_fac_ui (gamma, u - 1, rnd_mode);
  /* if bits_fac(...) > p (resp. p+1 for rounding to nearest),
     then gamma(x) cannot be exact in precision p (resp. p+1).
     FIXME: remove the test u < 44787929UL after changing bits_fac
     to return a mpz_t or mpfr_t. */

[wjamesjr]

Yep, that looks like the issue. Strange number to set as a cutoff! I will fall back to a python factorial algorithm with some bit slicing. Thanks for the quick response!

[wjamesjr]

Looks like the issue runs deeper than just gamma(). The below code computing a mpz factorial fails as well, with result inf. 44787927 still works. Is there some trick to converting very large mpz values to mpfr??

n = 44787928
gmpy2.get_context().emin = gmpy2.get_emin_min()
gmpy2.get_context().emax = gmpy2.get_emax_max()
print("emin", gmpy2.get_context().emin)
print("emax", gmpy2.get_context().emax)
print("max precision", gmpy2.get_max_precision())
gmpy2.get_context().precision = 100
print("precision set to", gmpy2.get_context().precision)
factn = gmpy2.fac(n)
print(factn.num_digits())
factn2 = mpfr(factn)
print(factn2)
exit(0)

[wjamesjr]

Looks like I have a hard limit on size of mpz that will convert to mpfr. The below code works. Change 2 29 to 2 30 and it goes to inf again. Should mpfr be limited in bit size?? I thought that's what precision provided.

gmpy2.get_context().emin = gmpy2.get_emin_min()
gmpy2.get_context().emax = gmpy2.get_emax_max()
print("emin", gmpy2.get_context().emin)
print("emax", gmpy2.get_context().emax)
print("max precision", gmpy2.get_max_precision())
gmpy2.get_context().precision = 100
print("precision set to", gmpy2.get_context().precision)
t1 = mpz(1) << (2 ** 29)
t2 = mpfr(t1)
print(t2)

[casevh]

It looks like a deeper issue:

t1=mpz(1) << (2*27) mpfr(t1) mpfr('1.196380724997376356710237763087067030291123782412927478906332372851427131104587e+40403562',256) _ mpfr('1.431326839145247872477712623353078898059627334067519357500412952958814141885275e+80807124',256) mpfr('2.048696520457526277391095958728021868321933030871131210018127684872659300422929e+161614248',256) _ mpfr('inf')

and

mpfr(123456.789)63483756.2 mpfr('3.15691251903978272292764153590394795083483376196033776445078487109637376386208e+323228495',256) mpfr(123456.789)63483756.3 mpfr('1.019562878429473795214570128278084403612032953011088359304266502177487364234525e+323228496',256) mpfr(123456.789)**63483756.4 mpfr('inf')

For reference

gmpy2.get_context() context(precision=256, real_prec=Default, imag_prec=Default, round=RoundToNearest, real_round=Default, imag_round=Default, emax=4611686018427387903, emin=-4611686018427387903, subnormalize=False, trap_underflow=False, underflow=False, trap_overflow=False, overflow=True, trap_inexact=False, inexact=True, trap_invalid=False, invalid=False, trap_erange=False, erange=False, trap_divzero=False, divzero=False, allow_complex=False, rational_division=False, allow_release_gil=False)

Time to dig into the MPFR source code.

On Fri, Jan 5, 2024 at 8:11 AM Bill James @.***> wrote:

Looks like I have a hard limit on size of mpz that will convert to mpfr. The below code works. Change 2 29 to 2 30 and it goes to inf again. Should mpfr be limited in bit size?? I thought that's what precision provided.

gmpy2.get_context().emin = gmpy2.get_emin_min() gmpy2.get_context().emax = gmpy2.get_emax_max() print("emin", gmpy2.get_context().emin) print("emax", gmpy2.get_context().emax) print("max precision", gmpy2.get_max_precision()) gmpy2.get_context().precision = 100 print("precision set to", gmpy2.get_context().precision) t1 = mpz(1) << (2 ** 29) t2 = mpfr(t1) print(t2)

— Reply to this email directly, view it on GitHub https://github.com/aleaxit/gmpy/issues/464#issuecomment-1878912852, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAMR23YSAMWNKZSKGSIAL2LYNAQ2BAVCNFSM6AAAAABBN54TPKVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTQNZYHEYTEOBVGI . You are receiving this because you commented.Message ID: @.***>

[wjamesjr]

Reported a bug to the MPFR development team. After some discussion they have decided to investigate. They say that they do not test for large bit precisions due to the amount of time required. But in this case they have decided to look into it to find if it is some int32 overflow or something like large multiplication issues and such.

[skirpichev]

Could we kept here a link to MPFR issue or maillist thread?

[wjamesjr]

https://sympa.inria.fr/sympa/arc/mpfr/2024-10/msg00002.html

Above is mail list thread hyperlink where I reported the bug.

[casevh]

I finally remembered this issue from a few years ago. It is a limitation of the MPFR on all platforms. This is from MPFR's gamma.c:

  if (u < 44787929UL && bits_fac (u - 1) <= p + (rnd_mode == MPFR_RNDN))
    /* bits_fac: lower bound on the number of bits of m,
       where gamma(x) = (u-1)! = m*2^e with m odd. */
    return mpfr_fac_ui (gamma, u - 1, rnd_mode);
  /* if bits_fac(...) > p (resp. p+1 for rounding to nearest),
     then gamma(x) cannot be exact in precision p (resp. p+1).
     FIXME: remove the test u < 44787929UL after changing bits_fac
     to return a mpz_t or mpfr_t. */

On Fri, Oct 11, 2024 at 7:49 PM Bill James @.***> wrote:

https://sympa.inria.fr/sympa/arc/mpfr/2024-10/msg00002.html

Above is mail list thread hyperlink where I reported the bug.

— Reply to this email directly, view it on GitHub https://github.com/aleaxit/gmpy/issues/464#issuecomment-2408318165, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAMR232SBYFJEGM7YIB4DL3Z3CEZ3AVCNFSM6AAAAABPZ2VIICVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDIMBYGMYTQMJWGU . You are receiving this because you commented.Message ID: @.***>

[wjamesjr]

The sample code that I gave to the MPFR development team had an issue. The c++ mpfr wrapper was using an int vs long argument for a call to convert digits to bits for precision. After changing that argument to a long and increasing emax and emin values the code worked. Which means that MPFR will work for bit precisions over 2 30, at least for the code I was testing. Now, I still can't get gmpy to work for the same algorithm for above the 2 30 bit precision threshold. I will submit test code as a separate issue.