nan from expression with very large numbers

[GoogleCodeExporter]

a = 3341925068
b = 1954900845
(mpfr(3) ** a - 1) / ((mpfr(2) ** a - 1) * 2 ** b)

What is the expected output? A number just less than one
What do you see instead? nan

Python 3.4.2 (v3.4.2:ab2c023a9432, Oct  6 2014, 22:16:31) [MSC v.1600 64 bit 
(AMD64)] on win32

'MPFR 3.1.1'

Is this how calculations are supposed to fail when numbers get too big?

Original issue reported on code.google.com by short.jo...@gmail.com on 13 Mar 2015 at 7:25

[GoogleCodeExporter]

mpfr(3)**a exceeds the default maximum exponent and return 'inf'. The same is 
true for mpfr(2)**a. 'inf'/'inf' returns 'nan', as expected.

The MPFR library stores the exponent in a C long which is 32 bits long on 
32-bit platforms and also on 64-bit Windows. On 64-bit Linux/Unix/OSX systems, 
a C long is 64 bits in length and the exponent range can be increased on those 
systems.

>>> gmpy2.get_context()
context(precision=53, real_prec=Default, imag_prec=Default,
        round=RoundToNearest, real_round=Default, imag_round=Default,
        emax=1073741823, emin=-1073741823,
        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,
        guard_bits=0,
        convert_exact=False,
        mpfr_divmod_exact=False)
>>> gmpy2.get_context().emax=gmpy2.get_emax_max()
>>> gmpy2.get_context().emin=gmpy2.get_emin_min()
>>> gmpy2.get_context()
context(precision=53, 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,
        guard_bits=0,
        convert_exact=False,
        mpfr_divmod_exact=False)
>>> 

The following is from a Linux system:

>>> import gmpy2
>>> from gmpy2 import mpfr,mpz
>>> gmpy2.get_context().emax=gmpy2.get_emax_max()
>>> gmpy2.get_context().emin=gmpy2.get_emin_min()
>>> a = 3341925068
>>> b = 1954900845
>>> (mpfr(3) ** a - 1) / ((mpfr(2) ** a - 1) * mpz(2) ** b)
mpfr('0.99999994506351741')

For all reasonable precision values, the "- 1" terms can be ignored and result 
can be calculated using logarithms. The precision does need to be increased 
slightly.

>>> gmpy2.get_context().precision=200
>>> from gmpy2 import exp, log
>>> exp(a * log(3) - ((a+b) * log(2)))
mpfr('0.9999999450635174326119685396058203532977346729182603488114124',200)
>>> 

Let me know if you still have questions.

Thanks,
casevh

Original comment by casevh on 14 Mar 2015 at 5:58

• Changed state: Started

[GoogleCodeExporter]

Many thanks, Case, for that clear explanation and the suggestion.

Original comment by short.jo...@gmail.com on 14 Mar 2015 at 6:18

[GoogleCodeExporter]

You're welcome. I hope you enjoy using gmpy2.

casevh

Original comment by casevh on 14 Mar 2015 at 9:07

• Changed state: Done