Ticket #13033: 13033_pari_to_mpfr.patch

File 13033_pari_to_mpfr.patch, 12.6 KB (added by jdemeyer, 10 years ago)

  • sage/rings/real_mpfi.pyx 

    # HG changeset patch
# User Jeroen Demeyer <jdemeyer@cage.ugent.be>
# Date 1338169914 25200
# Node ID dfd3b8044a4d904b60bba6db3a89111290cbc657
# Parent  8709b976e906b35b6174d82ea0c3776d2f36e24a
Increase exponent range of real numbers to the maximal possible range

diff --git a/sage/rings/real_mpfi.pyx b/sage/rings/real_mpfi.pyx
    a b  
    11"""
    2 Field of Arbitrary Precision Real Intervals
     2Arbitrary Precision Real Intervals
    33
    44AUTHORS:
    55
    … …  
    21222122            sage: a.fp_rank_diameter()
    21232123            30524
    21242124            sage: (RIF(sqrt(2)) - RIF(sqrt(2))).fp_rank_diameter()
    2125             9671406088542672151117826
    2126        
     2125            9671406088542672151117826            # 32-bit
     2126            41538374868278620559869609387229186  # 64-bit
     2127
    21272128        Just because we have the best possible interval, doesn't mean the
    21282129        interval is actually small::
    2129        
    2130             sage: a = RIF(pi)^1234567890; a
    2131             [2.0985787164673874e323228496 .. +infinity]
     2130
     2131            sage: a = RIF(pi)^12345678901234567890; a
     2132            [2.0985787164673874e323228496 .. +infinity]            # 32-bit
     2133            [5.8756537891115869e1388255822130839282 .. +infinity]  # 64-bit
    21322134            sage: a.fp_rank_diameter()
    21332135            1
    21342136        """
    21352137        return self.lower().fp_rank_delta(self.upper())
    21362138
    21372139    def is_exact(self):
    2138         return mpfr_equal_p(&self.value.left, &self.value.right)           
     2140        return mpfr_equal_p(&self.value.left, &self.value.right)
    21392141
    21402142    def magnitude(self):
    21412143        """

  • sage/rings/real_mpfr.pyx 

    diff --git a/sage/rings/real_mpfr.pyx b/sage/rings/real_mpfr.pyx
    a b  
    11r"""
    2 Field of Arbitrary Precision Real Numbers
     2Arbitrary Precision Real Numbers
    33
    44AUTHORS:
    55
    … …  
    2424
    2525- Robert Bradshaw (2009-09): decimal literals, optimizations
    2626
     27- Jeroen Demeyer (2012-05-27): set the MPFR exponent range to the
     28  maximal possible value (:trac:`13033`)
     29
    2730This is a binding for the MPFR arbitrary-precision floating point
    2831library.
    2932
    … …  
    3538In Sage (as in MPFR), floating-point numbers of precision
    3639`p` are of the form `s m 2^{e-p}`, where
    3740`s \in \{-1, 1\}`, `2^{p-1} \leq m < 2^p`, and
    38 `-2^{30} + 1 \leq e \leq 2^{30} - 1`; plus the special
    39 values ``+0``, ``-0``,
    40 ``+infinity``, ``-infinity``, and
    41 ``NaN`` (which stands for Not-a-Number).
     41`-2^B + 1 \leq e \leq 2^B - 1` where `B = 30` on 32-bit systems
     42and `B = 62` on 64-bit systems;
     43additionally, there are the special values ``+0``, ``-0``,
     44``+infinity``, ``-infinity`` and ``NaN`` (which stands for Not-a-Number).
    4245
    4346Operations in this module which are direct wrappers of MPFR
    4447functions are "correctly rounded"; we briefly describe what this
    … …  
    180183
    181184def mpfr_prec_min():
    182185    """
    183     Return the mpfr variable MPFR_PREC_MIN.
     186    Return the mpfr variable ``MPFR_PREC_MIN``.
     187
     188    EXAMPLES::
     189
     190        sage: from sage.rings.real_mpfr import mpfr_prec_min
     191        sage: mpfr_prec_min()
     192        2
     193        sage: R = RealField(2)
     194        sage: R(2) + R(1)
     195        3.0
     196        sage: R(4) + R(1)
     197        4.0
     198        sage: R = RealField(1)
     199        Traceback (most recent call last):
     200        ...
     201        ValueError: prec (=1) must be >= 2 and <= 2147483391
    184202    """
    185203    return MPFR_PREC_MIN
    186204
    … …  
    199217        sage: R = RealField(2^31-256)
    200218        Traceback (most recent call last):
    201219        ...
    202         ValueError: prec (=2147483392) must be >= 2 and <= 2147483391.
     220        ValueError: prec (=2147483392) must be >= 2 and <= 2147483391
    203221    """
    204222    global MY_MPFR_PREC_MAX
    205223    return MY_MPFR_PREC_MAX
    206224
     225def mpfr_get_exp_min():
     226    """
     227    Return the current minimal exponent for MPFR numbers.
     228
     229    EXAMPLES::
     230
     231        sage: from sage.rings.real_mpfr import mpfr_get_exp_min
     232        sage: mpfr_get_exp_min()
     233        -1073741823            # 32-bit
     234        -4611686018427387903   # 64-bit
     235        sage: 0.5 >> (-mpfr_get_exp_min())
     236        2.38256490488795e-323228497            # 32-bit
     237        8.50969131174084e-1388255822130839284  # 64-bit
     238        sage: 0.5 >> (-mpfr_get_exp_min()+1)
     239        0.000000000000000
     240    """
     241    return mpfr_get_emin()
     242
     243def mpfr_get_exp_max():
     244    """
     245    Return the current maximal exponent for MPFR numbers.
     246
     247    EXAMPLES::
     248
     249        sage: from sage.rings.real_mpfr import mpfr_get_exp_max
     250        sage: mpfr_get_exp_max()
     251        1073741823            # 32-bit
     252        4611686018427387903   # 64-bit
     253        sage: 0.5 << mpfr_get_exp_max()
     254        1.04928935823369e323228496            # 32-bit
     255        2.93782689455579e1388255822130839282  # 64-bit
     256        sage: 0.5 << (mpfr_get_exp_max()+1)
     257        +infinity
     258    """
     259    return mpfr_get_emax()
     260
     261def mpfr_set_exp_min(mp_exp_t e):
     262    """
     263    Set the minimal exponent for MPFR numbers.
     264
     265    EXAMPLES::
     266
     267        sage: from sage.rings.real_mpfr import mpfr_get_exp_min, mpfr_set_exp_min
     268        sage: old = mpfr_get_exp_min()
     269        sage: mpfr_set_exp_min(-1000)
     270        sage: 0.5 >> 1000
     271        4.66631809251609e-302
     272        sage: 0.5 >> 1001
     273        0.000000000000000
     274        sage: mpfr_set_exp_min(old)
     275        sage: 0.5 >> 1001
     276        2.33315904625805e-302
     277    """
     278    if mpfr_set_emin(e) != 0:
     279        raise OverflowError("bad value for mpfr_set_exp_min()")
     280
     281def mpfr_set_exp_max(mp_exp_t e):
     282    """
     283    Set the maximal exponent for MPFR numbers.
     284
     285    EXAMPLES::
     286
     287        sage: from sage.rings.real_mpfr import mpfr_get_exp_max, mpfr_set_exp_max
     288        sage: old = mpfr_get_exp_max()
     289        sage: mpfr_set_exp_max(1000)
     290        sage: 0.5 << 1000
     291        5.35754303593134e300
     292        sage: 0.5 << 1001
     293        +infinity
     294        sage: mpfr_set_exp_max(old)
     295        sage: 0.5 << 1001
     296        1.07150860718627e301
     297    """
     298    if mpfr_set_emax(e) != 0:
     299        raise OverflowError("bad value for mpfr_set_exp_max()")
     300
     301def mpfr_get_exp_min_min():
     302    """
     303    Get the minimal value allowed for :func:`mpfr_set_exp_min`.
     304
     305    EXAMPLES::
     306
     307        sage: from sage.rings.real_mpfr import mpfr_get_exp_min_min, mpfr_set_exp_min
     308        sage: mpfr_get_exp_min_min()
     309        -1073741823            # 32-bit
     310        -4611686018427387903   # 64-bit
     311
     312    This is really the minimal value allowed::
     313
     314        sage: mpfr_set_exp_min(mpfr_get_exp_min_min() - 1)
     315        Traceback (most recent call last):
     316        ...
     317        OverflowError: bad value for mpfr_set_exp_min()
     318    """
     319    return mpfr_get_emin_min()
     320
     321def mpfr_get_exp_max_max():
     322    """
     323    Get the maximal value allowed for :func:`mpfr_set_exp_max`.
     324
     325    EXAMPLES::
     326
     327        sage: from sage.rings.real_mpfr import mpfr_get_exp_max_max, mpfr_set_exp_max
     328        sage: mpfr_get_exp_max_max()
     329        1073741823            # 32-bit
     330        4611686018427387903   # 64-bit
     331
     332    This is really the maximal value allowed::
     333
     334        sage: mpfr_set_exp_max(mpfr_get_exp_max_max() + 1)
     335        Traceback (most recent call last):
     336        ...
     337        OverflowError: bad value for mpfr_set_exp_max()
     338    """
     339    return mpfr_get_emax_max()
     340
     341# On Sage startup, set the exponent range to the maximum allowed
     342mpfr_set_exp_min(mpfr_get_emin_min())
     343mpfr_set_exp_max(mpfr_get_emax_max())
     344
    207345#*****************************************************************************
    208346#
    209347#       Real Field
    210348#
    211349#*****************************************************************************
    212 # The real field is in Pyrex, so mpfr elements will have access to
     350# The real field is in Cython, so mpfr elements will have access to
    213351# their parent via direct C calls, which will be faster.
    214352
    215353_rounding_modes = ['RNDN', 'RNDZ', 'RNDU', 'RNDD']
    … …  
    295433        global MY_MPFR_PREC_MAX
    296434        cdef RealNumber rn
    297435        if prec < MPFR_PREC_MIN or prec > MY_MPFR_PREC_MAX:
    298             raise ValueError, "prec (=%s) must be >= %s and <= %s."%(
     436            raise ValueError, "prec (=%s) must be >= %s and <= %s"%(
    299437                prec, MPFR_PREC_MIN, MY_MPFR_PREC_MAX)
    300438        self.__prec = prec
    301439        if not isinstance(rnd, str):
    … …  
    12181356
    12191357        mpz_import(mantissa, lg(g) - 2, 1, wordsize/8, 0, 0, &g[2])
    12201358
    1221         cdef int exponent
     1359        cdef mp_exp_t exponent
    12221360        exponent = expo(g)
    12231361
    12241362        # Round to nearest for best results when setting a low-precision
    … …  
    17511889            sage: RR(0).nextbelow().nextbelow().fp_rank()
    17521890            -2
    17531891            sage: RR(1).fp_rank()
    1754             4835703278458516698824705
     1892            4835703278458516698824705            # 32-bit
     1893            20769187434139310514121985316880385  # 64-bit
    17551894            sage: RR(-1).fp_rank()
    1756             -4835703278458516698824705
     1895            -4835703278458516698824705            # 32-bit
     1896            -20769187434139310514121985316880385  # 64-bit
    17571897            sage: RR(1).fp_rank() - RR(1).nextbelow().fp_rank()
    17581898            1
    17591899            sage: RR(-infinity).fp_rank()
    1760             -9671406552413433770278913
     1900            -9671406552413433770278913            # 32-bit
     1901            -41538374868278621023740371006390273  # 64-bit
    17611902            sage: RR(-infinity).fp_rank() - RR(-infinity).nextabove().fp_rank()
    17621903            -1
    17631904        """
    … …  
    17661907
    17671908        cdef Integer z = PY_NEW(Integer)
    17681909
    1769         cdef mp_exp_t EXP_MIN = -(1<<30) + 1
    1770         cdef mp_exp_t EXP_MAX = (1<<30) - 1
     1910        cdef mp_exp_t EXP_MIN = mpfr_get_exp_min()
     1911        cdef mp_exp_t EXP_MAX = mpfr_get_exp_max()
    17711912        # fp_rank(0.0) = 0
    17721913        # fp_rank(m*2^e-p) = (m-2^{p-1})+(e-EXP_MIN)*2^{p-1}+1
    17731914        #                  = m+(e-EXP_MIN-1)*2^{p-1}+1
    … …  
    18371978        There are lots of floating-point numbers around 0::
    18381979       
    18391980            sage: R2(-1).fp_rank_delta(R2(1))
    1840             4294967298
    1841         """
    1842 
     1981            4294967298            # 32-bit
     1982            18446744073709551618  # 64-bit
     1983        """
    18431984        # We create the API for forward compatibility, because it can have
    18441985        # a (somewhat) more efficient implementation than this; but for now,
    18451986        # we just go with the stupid implementation.
    … …  
    23482489            sage: (1.0).nexttoward(RR('-infinity')).str(truncate=False)
    23492490            '0.99999999999999989'
    23502491            sage: RR(infinity).nexttoward(0)
    2351             2.09857871646739e323228496
     2492            2.09857871646739e323228496            # 32-bit
     2493            5.87565378911159e1388255822130839282  # 64-bit
    23522494            sage: RR(pi).str(truncate=False)
    23532495            '3.1415926535897931'
    23542496            sage: RR(pi).nexttoward(22/7).str(truncate=False)
    … …  
    23762518        EXAMPLES::
    23772519       
    23782520            sage: RR('-infinity').nextabove()
    2379             -2.09857871646739e323228496
     2521            -2.09857871646739e323228496            # 32-bit
     2522            -5.87565378911159e1388255822130839282  # 64-bit
    23802523            sage: RR(0).nextabove()
    2381             2.38256490488795e-323228497
     2524            2.38256490488795e-323228497            # 32-bit
     2525            8.50969131174084e-1388255822130839284  # 64-bit
    23822526            sage: RR('+infinity').nextabove()
    23832527            +infinity
    23842528            sage: RR(-sqrt(2)).str(truncate=False)
    … …  
    24022546            sage: RR('-infinity').nextbelow()
    24032547            -infinity
    24042548            sage: RR(0).nextbelow()
    2405             -2.38256490488795e-323228497
     2549            -2.38256490488795e-323228497            # 32-bit
     2550            -8.50969131174084e-1388255822130839284  # 64-bit
    24062551            sage: RR('+infinity').nextbelow()
    2407             2.09857871646739e323228496
     2552            2.09857871646739e323228496              # 32-bit
     2553            5.87565378911159e1388255822130839282    # 64-bit
    24082554            sage: RR(-sqrt(2)).str(truncate=False)
    24092555            '-1.4142135623730951'
    24102556            sage: RR(-sqrt(2)).nextbelow().str(truncate=False)
    … …  
    25782724            2
    25792725            sage: RealField(100)(0.0)._pari_().sizeword()
    25802726            2
     2727
     2728        Check that the largest and smallest exponents representable by
     2729        PARI convert correctly::
     2730
     2731            sage: a = pari(0.5) << (sys.maxint+1)/4
     2732            sage: RR(a) >> (sys.maxint+1)/4
     2733            0.500000000000000
     2734            sage: a = pari(0.5) >> (sys.maxint-3)/4
     2735            sage: RR(a) << (sys.maxint-3)/4
     2736            0.500000000000000
    25812737        """
    25822738        # This uses interfaces of MPFR and PARI which are documented
    25832739        # (and not marked subject-to-change).  It could be faster
    … …  
    25982754        cdef int rounded_prec
    25992755        rounded_prec = (self.prec() + wordsize - 1) & ~(wordsize - 1)
    26002756
    2601         # Yes, assigning to self works fine, even in Pyrex.
     2757        # Yes, assigning to self works fine, even in Cython.
    26022758        if rounded_prec > prec:
    26032759            self = RealField(rounded_prec)(self)
    26042760