Ticket #15337: 15337_ulp.patch

File 15337_ulp.patch, 14.3 KB (added by Jeroen Demeyer, 9 years ago)
  • sage/libs/mpfr.pxd

    # HG changeset patch
    # User Jeroen Demeyer <jdemeyer@cage.ugent.be>
    # Date 1382968568 -3600
    # Node ID 55f5467a877620e3b808d5e2620c93dcd2662a6e
    # Parent  83ce0bc2816413611fd533af77149747fa8b152f
    Speed up ulp() method, add epsilon() method, add field arguments
    
    diff --git a/sage/libs/mpfr.pxd b/sage/libs/mpfr.pxd
    a b  
    44cdef extern from "mpfr.h":
    55    ctypedef struct __mpfr_struct:
    66        pass
    7        
     7
    88    #ctypedef __mpfr_struct mpfr_t[1]
    99    ctypedef __mpfr_struct* mpfr_t
    10     ctypedef mpfr_t mpfr_ptr 
     10    ctypedef mpfr_t mpfr_ptr
    1111    ctypedef mpfr_t mpfr_srcptr
    1212    ctypedef enum mpfr_rnd_t:
    1313        GMP_RNDN = 0
     
    1616        GMP_RNDD = 3
    1717        GMP_RND_MAX = 4
    1818        GMP_RNDNA = -1
    19    
     19
    2020    ctypedef mpfr_rnd_t mp_rnd_t
    2121    ctypedef long mp_prec_t
    2222
     
    3535    int mpfr_set (mpfr_t rop, mpfr_t op, mp_rnd_t rnd)
    3636    int mpfr_set_ui (mpfr_t rop, unsigned long int op, mp_rnd_t rnd)
    3737    int mpfr_set_si (mpfr_t rop, long int op, mp_rnd_t rnd)
    38     # int mpfr_set_uj (mpfr_t rop, uintmax_t op, mp_rnd_t rnd)
    39     # int mpfr_set_sj (mpfr_t rop, intmax_t op, mp_rnd_t rnd)
    4038    int mpfr_set_d (mpfr_t rop, double op, mp_rnd_t rnd)
    4139    int mpfr_set_ld (mpfr_t rop, long double op, mp_rnd_t rnd)
    4240    # int mpfr_set_decimal64 (mpfr_t rop, _Decimal64 op, mp_rnd_t rnd)
     
    4543    # int mpfr_set_f (mpfr_t rop, mpf_t op, mp_rnd_t rnd)
    4644    int mpfr_set_ui_2exp (mpfr_t rop, unsigned long int op, mp_exp_t e, mp_rnd_t rnd)
    4745    int mpfr_set_si_2exp (mpfr_t rop, long int op, mp_exp_t e, mp_rnd_t rnd)
    48     # int mpfr_set_uj_2exp (mpfr_t rop, uintmax_t op, intmax_t e, mp_rnd_t rnd)
    49     # int mpfr_set_sj_2exp (mpfr_t rop, intmax_t op, intmax_t e, mp_rnd_t rnd)
    5046    int mpfr_set_str (mpfr_t rop,  char *s, int base, mp_rnd_t rnd)
    5147    int mpfr_strtofr (mpfr_t rop, char *nptr, char **endptr, int base, mp_rnd_t rnd)
    5248    void mpfr_set_inf (mpfr_t x, int sign)
    5349    void mpfr_set_nan (mpfr_t x)
     50    void mpfr_set_zero (mpfr_t x, int sign)
    5451    void mpfr_swap (mpfr_t x, mpfr_t y)
    5552
    5653    # Combined Initialization and Assignment Functions
     
    7269    long double mpfr_get_ld_2exp (long *exp, mpfr_t op, mp_rnd_t rnd)
    7370    long mpfr_get_si (mpfr_t op, mp_rnd_t rnd)
    7471    unsigned long mpfr_get_ui (mpfr_t op, mp_rnd_t rnd)
    75     # intmax_t mpfr_get_sj (mpfr_t op, mp_rnd_t rnd)
    76     # uintmax_t mpfr_get_uj (mpfr_t op, mp_rnd_t rnd)
    7772    mp_exp_t mpfr_get_z_exp (mpz_t rop, mpfr_t op)
    7873    void mpfr_get_z (mpz_t rop, mpfr_t op, mp_rnd_t rnd)
    7974    # int mpfr_get_f (mpf_t rop, mpfr_t op, mp_rnd_t rnd)
     
    148143    bint mpfr_inf_p (mpfr_t op)
    149144    bint mpfr_number_p (mpfr_t op)
    150145    bint mpfr_zero_p (mpfr_t op)
     146    bint mpfr_regular_p (mpfr_t op)
    151147    int mpfr_sgn (mpfr_t op)
    152148    bint mpfr_greater_p (mpfr_t op1, mpfr_t op2)
    153149    bint mpfr_greaterequal_p (mpfr_t op1, mpfr_t op2)
     
    188184    int mpfr_log1p (mpfr_t rop, mpfr_t op, mp_rnd_t rnd)
    189185    int mpfr_expm1 (mpfr_t rop, mpfr_t op, mp_rnd_t rnd)
    190186    int mpfr_eint (mpfr_t y, mpfr_t x, mp_rnd_t rnd)
     187    int mpfr_li2 (mpfr_t y, mpfr_t x, mp_rnd_t rnd)
    191188    int mpfr_gamma (mpfr_t rop, mpfr_t op, mp_rnd_t rnd)
    192189    int mpfr_lngamma (mpfr_t rop, mpfr_t op, mp_rnd_t rnd)
    193190    int mpfr_lgamma (mpfr_t rop, int *signp, mpfr_t op, mp_rnd_t rnd)
  • sage/rings/real_double.pyx

    diff --git a/sage/rings/real_double.pyx b/sage/rings/real_double.pyx
    a b  
    626626
    627627    def ulp(self):
    628628        """
    629         Return the unit of least precision of ``self``, which is the weight of
    630         the least significant bit of ``self``. Unless ``self`` is exactly a
    631         power of two, it is gap between this number and the next closest
    632         distinct  number that can be represented.
     629        Returns the unit of least precision of ``self``, which is the
     630        weight of the least significant bit of ``self``. This is always
     631        a strictly positive number. It is also the gap between this
     632        number and the closest number with larger absolute value that
     633        can be represented.
    633634
    634635        EXAMPLES::
    635636
     
    649650            sage: b - b.ulp() == b
    650651            False
    651652
    652         Adding or subtracting something less than half an ulp never
    653         gives the same number (unless the number is exactly a power of
    654         2 and subtracting an ulp decreases the ulp)::
     653        Since the default rounding mode is round-to-nearest, adding or
     654        subtracting something less than half an ulp always gives the
     655        same number, unless the result has a smaller ulp. The latter
     656        can only happen if the input number is (up to sign) exactly a
     657        power of 2::
    655658
    656659            sage: a - a.ulp()/3 == a
    657660            True
     
    680683            4.94065645841e-324
    681684            sage: RDF(0).ulp()/2
    682685            0.0
    683            
     686
    684687        Some special values::
    685688
    686689            sage: a = RDF(1)/RDF(0); a
     
    689692            +infinity
    690693            sage: (-a).ulp()
    691694            +infinity
    692             sage: a = RR('nan')
     695            sage: a = RDF('nan')
    693696            sage: a.ulp() is a
    694697            True
     698
     699        The ulp method works correctly with small numbers::
     700
     701            sage: u = RDF(0).ulp()
     702            sage: u.ulp() == u
     703            True
     704            sage: x = u * (2^52-1)  # largest denormal number
     705            sage: x.ulp() == u
     706            True
     707            sage: x = u * 2^52  # smallest normal number
     708            sage: x.ulp() == u
     709            True
     710           
    695711        """
    696712        # First, check special values
    697713        if self._value == 0:
     
    704720        # Normal case
    705721        cdef int e
    706722        frexp(self._value, &e)
    707         return RealDoubleElement(ldexp(1.0, e-53))
     723        e -= 53
     724        # Correction for denormals
     725        if e < -1074:
     726            e = -1074
     727        return RealDoubleElement(ldexp(1.0, e))
    708728
    709729    def real(self):
    710730        """
     
    18701890
    18711891            sage: R = RealField(128)
    18721892            sage: def check_error(x):
    1873             ...     x = RDF(x)
    1874             ...     log_RDF = x.log()
    1875             ...     log_RR = R(x).log()
    1876             ...     diff = R(log_RDF) - log_RR
    1877             ...     if abs(diff) <= log_RDF.ulp():
    1878             ...         return True
    1879             ...     print "logarithm check failed for %s (diff = %s ulp)"% \
    1880             ...         (x, diff/log_RDF.ulp())
    1881             ...     return False
     1893            ....:   x = RDF(x)
     1894            ....:   log_RDF = x.log()
     1895            ....:   log_RR = R(x).log()
     1896            ....:   diff = R(log_RDF) - log_RR
     1897            ....:   if abs(diff) < log_RDF.ulp():
     1898            ....:       return True
     1899            ....:   print "logarithm check failed for %s (diff = %s ulp)"% \
     1900            ....:       (x, diff/log_RDF.ulp())
     1901            ....:   return False
    18821902            sage: all( check_error(2^x) for x in range(-100,100) )
    18831903            True
    18841904            sage: all( check_error(x) for x in sxrange(0.01, 2.00, 0.01) )
  • sage/rings/real_mpfr.pyx

    diff --git a/sage/rings/real_mpfr.pyx b/sage/rings/real_mpfr.pyx
    a b  
    23882388        """
    23892389        return self
    23902390
    2391     def ulp(self):
    2392         """
    2393         Returns the unit of least precision of ``self``, which is the weight of
    2394         the least significant bit of ``self``. Unless ``self`` is exactly a
    2395         power of two, it is gap between this number and the next closest
    2396         distinct number that can be represented.
     2391    def ulp(self, field=None):
     2392        """
     2393        Returns the unit of least precision of ``self``, which is the
     2394        weight of the least significant bit of ``self``. This is always
     2395        a strictly positive number. It is also the gap between this
     2396        number and the closest number with larger absolute value that
     2397        can be represented.
     2398
     2399        INPUT:
     2400
     2401        - ``field`` -- :class:`RealField` used as parent of the result.
     2402          If not specified, use ``parent(self)``.
     2403
     2404        .. NOTE::
     2405
     2406            The ulp of zero is defined as the smallest representable
     2407            positive number. For extremely small numbers, underflow
     2408            occurs and the output is also the smallest representable
     2409            positive number (the rounding mode is ignored, this
     2410            computation is done by rounding towards +infinity).
     2411
     2412        .. SEEALSO::
     2413
     2414            :meth:`epsilon` for a scale-invariant version of this.
    23972415
    23982416        EXAMPLES::
    23992417
    24002418            sage: a = 1.0
     2419            sage: a.ulp()
     2420            2.22044604925031e-16
     2421            sage: (-1.5).ulp()
     2422            2.22044604925031e-16
    24012423            sage: a + a.ulp() == a
    24022424            False
    24032425            sage: a + a.ulp()/2 == a
     
    24072429            sage: b = a + a.ulp()
    24082430            sage: (a+b)/2 in [a,b]
    24092431            True
     2432
     2433        The ulp of zero is the smallest non-zero number::
     2434
     2435            sage: a = RR(0).ulp()
     2436            sage: a
     2437            2.38256490488795e-323228497            # 32-bit
     2438            8.50969131174084e-1388255822130839284  # 64-bit
     2439            sage: a.fp_rank()
     2440            1
     2441
     2442        The ulp of very small numbers results in underflow, so the
     2443        smallest non-zero number is returned instead::
     2444
     2445            sage: a.ulp() == a
     2446            True
     2447
     2448        We use a different field::
     2449
     2450            sage: a = RealField(256).pi()
     2451            sage: a.ulp()
     2452            3.454467422037777850154540745120159828446400145774512554009481388067436721265e-77
     2453            sage: e = a.ulp(RealField(64))
     2454            sage: e
     2455            3.45446742203777785e-77
     2456            sage: parent(e)
     2457            Real Field with 64 bits of precision
     2458            sage: e = a.ulp(QQ)
     2459            Traceback (most recent call last):
     2460            ...
     2461            TypeError: field argument must be a RealField
     2462
     2463        For infinity and NaN, we get back positive infinity and NaN::
    24102464           
    24112465            sage: a = RR(infinity)
    24122466            sage: a.ulp()
     
    24142468            sage: (-a).ulp()
    24152469            +infinity
    24162470            sage: a = RR('nan')
    2417             sage: a.ulp() is a
    2418             True
    2419         """
    2420         cdef RealNumber x
    2421         if mpfr_nan_p(self.value):
    2422             return self
    2423         elif mpfr_inf_p(self.value):
    2424             if mpfr_sgn(self.value) == 1:
    2425                 return self
    2426             else:
    2427                 return -self
     2471            sage: a.ulp()
     2472            NaN
     2473            sage: parent(RR('nan').ulp(RealField(42)))
     2474            Real Field with 42 bits of precision
     2475        """
     2476        cdef RealField_class _parent
     2477        if field is None:
     2478            _parent = self._parent
    24282479        else:
    2429             x = self._new()
    2430             mpfr_set(x.value, self.value, (<RealField_class>self._parent).rnd)
    2431             if mpfr_sgn(self.value) == 1:
     2480            try:
     2481                _parent = field
     2482            except TypeError:
     2483                raise TypeError("field argument must be a RealField")
     2484
     2485        cdef RealNumber x = _parent._new()
     2486        cdef mp_exp_t e
     2487        if mpfr_regular_p(self.value):
     2488            # Non-zero number
     2489            e = mpfr_get_exp(self.value) - mpfr_get_prec(self.value)
     2490            # Round up in case of underflow
     2491            mpfr_set_ui_2exp(x.value, 1, e, GMP_RNDU)
     2492        else:
     2493            # Special cases: zero, infinity, NaN
     2494            if mpfr_zero_p(self.value):
     2495                mpfr_set_zero(x.value, 1)
    24322496                mpfr_nextabove(x.value)
    2433             else:
    2434                 mpfr_nextbelow(x.value)
    2435             mpfr_sub(x.value, x.value, self.value, (<RealField_class>self._parent).rnd)
    2436             mpfr_abs(x.value, x.value, (<RealField_class>self._parent).rnd)
    2437             return x
    2438            
     2497            elif mpfr_inf_p(self.value):
     2498                mpfr_set_inf(x.value, 1)
     2499            else: # NaN
     2500                mpfr_set_nan(x.value)
     2501        return x
     2502
     2503    def epsilon(self, field=None):
     2504        """
     2505        Returns ``abs(self)`` divided by `2^b` where `b` is the
     2506        precision in bits of ``self``. Equivalently, return
     2507        ``abs(self)`` multiplied by the :meth:`ulp` of 1.
     2508
     2509        This is a scale-invariant version of :meth:`ulp` and it lies
     2510        in `[u/2, u)` where `u` is ``self.ulp()`` (except in the case
     2511        of zero or underflow).
     2512
     2513        INPUT:
     2514
     2515        - ``field`` -- :class:`RealField` used as parent of the result.
     2516          If not specified, use ``parent(self)``.
     2517
     2518        OUTPUT:
     2519
     2520        ``field(self.abs() / 2^self.precision())``
     2521
     2522        EXAMPLES::
     2523
     2524            sage: RR(2^53).epsilon()
     2525            1.00000000000000
     2526            sage: RR(0).epsilon()
     2527            0.000000000000000
     2528            sage: a = RR.pi()
     2529            sage: a.epsilon()
     2530            3.48786849800863e-16
     2531            sage: a.ulp()/2, a.ulp()
     2532            (2.22044604925031e-16, 4.44089209850063e-16)
     2533            sage: a / 2^a.precision()
     2534            3.48786849800863e-16
     2535            sage: (-a).epsilon()
     2536            3.48786849800863e-16
     2537
     2538        We use a different field::
     2539
     2540            sage: a = RealField(256).pi()
     2541            sage: a.epsilon()
     2542            2.713132368784788677624750042896586252980746500631892201656843478528498954308e-77
     2543            sage: e = a.epsilon(RealField(64))
     2544            sage: e
     2545            2.71313236878478868e-77
     2546            sage: parent(e)
     2547            Real Field with 64 bits of precision
     2548            sage: e = a.epsilon(QQ)
     2549            Traceback (most recent call last):
     2550            ...
     2551            TypeError: field argument must be a RealField
     2552
     2553        Special values::
     2554
     2555            sage: RR('nan').epsilon()
     2556            NaN
     2557            sage: parent(RR('nan').epsilon(RealField(42)))
     2558            Real Field with 42 bits of precision
     2559            sage: RR('+Inf').epsilon()
     2560            +infinity
     2561            sage: RR('-Inf').epsilon()
     2562            +infinity
     2563        """
     2564        cdef RealField_class _parent
     2565        if field is None:
     2566            _parent = self._parent
     2567        else:
     2568            try:
     2569                _parent = field
     2570            except TypeError:
     2571                raise TypeError("field argument must be a RealField")
     2572
     2573        cdef RealNumber x = _parent._new()
     2574        mpfr_div_2exp(x.value, self.value, mpfr_get_prec(self.value), _parent.rnd)
     2575        mpfr_abs(x.value, x.value, _parent.rnd)
     2576        return x
     2577
    24392578
    24402579    ###################
    24412580    # Rounding etc