Ticket #13134: 13134_ulp.patch

File 13134_ulp.patch, 3.8 KB (added by jdemeyer, 10 years ago)
  • sage/rings/real_double.pyx

    # HG changeset patch
    # User Jeroen Demeyer <jdemeyer@cage.ugent.be>
    # Date 1340100027 -7200
    # Node ID 0a193a60c44d25f1c8e52a09f49752b3f9af5d1f
    # Parent  840bb929a987235cfb6ccfd90b290a76f18bfefd
    Fix ulp() in RealDoubleField
    
    diff --git a/sage/rings/real_double.pyx b/sage/rings/real_double.pyx
    a b  
    599599       
    600600        EXAMPLES::
    601601 
    602             sage: a = RDF(1)
     602            sage: a = RDF(pi)
     603            sage: a.ulp()
     604            4.4408920985e-16
     605            sage: b = a + a.ulp()
     606
     607        Adding or subtracting an ulp always gives a different number::
     608
     609            sage: a + a.ulp() == a
     610            False
    603611            sage: a - a.ulp() == a
    604612            False
    605             sage: a - a.ulp()/2 == a
     613            sage: b + b.ulp() == b
     614            False
     615            sage: b - b.ulp() == b
     616            False
     617
     618        Adding or subtracting something less than half an ulp never
     619        gives the same number (unless the number is exactly a power of
     620        2 and subtracting an ulp decreases the ulp)::
     621
     622            sage: a - a.ulp()/3 == a
    606623            True
     624            sage: a + a.ulp()/3 == a
     625            True
     626            sage: b - b.ulp()/3 == b
     627            True
     628            sage: b + b.ulp()/3 == b
     629            True
     630            sage: c = RDF(1)
     631            sage: c - c.ulp()/3 == c
     632            False
     633            sage: c.ulp()
     634            2.22044604925e-16
     635            sage: (c - c.ulp()).ulp()
     636            1.11022302463e-16
    607637           
    608             sage: a = RDF.pi()
    609             sage: b = a + a.ulp()
    610             sage: (a+b)/2 in [a,b]
    611             True
     638        The ulp is always positive::
     639
     640            sage: RDF(-1).ulp()
     641            2.22044604925e-16
     642
     643        The ulp of zero is the smallest positive number in RDF::
     644
     645            sage: RDF(0).ulp()
     646            4.94065645841e-324
     647            sage: RDF(0).ulp()/2
     648            0.0
    612649           
     650        Some special values::
     651
    613652            sage: a = RDF(1)/RDF(0); a
    614653            +infinity
    615654            sage: a.ulp()
     
    620659            sage: a.ulp() is a
    621660            True
    622661        """
    623         cdef int e, v = gsl_isinf(self._value)
     662        # First, check special values
     663        if self._value == 0:
     664            return RealDoubleElement(ldexp(1.0, -1074))
    624665        if gsl_isnan(self._value):
    625666            return self
    626         elif self._value == 0:
    627             return RealDoubleElement(ldexp(1.0, e-1082))
    628         elif v == 1:
    629             return self
    630         elif v == -1:
    631             return -self
    632         else:
    633             frexp(self._value, &e)
    634             return RealDoubleElement(ldexp(1.0, e-54))
     667        if gsl_isinf(self._value):
     668            return self.abs()
     669
     670        # Normal case
     671        cdef int e
     672        frexp(self._value, &e)
     673        return RealDoubleElement(ldexp(1.0, e-53))
    635674
    636675    def real(self):
    637676        """
     
    17311770        Make sure that we can take the log of small numbers accurately
    17321771        and the fix doesn't break preexisting values (:trac:`12557`)::
    17331772
     1773            sage: R = RealField(128)
    17341774            sage: def check_error(x):
    17351775            ...     x = RDF(x)
    17361776            ...     log_RDF = x.log()
    1737             ...     log_exact = log(x.n(1000))
    1738             ...     return abs(log_RDF-log_exact) <= log_RDF.ulp()
     1777            ...     log_RR = R(x).log()
     1778            ...     diff = R(log_RDF) - log_RR
     1779            ...     if abs(diff) <= log_RDF.ulp():
     1780            ...         return True
     1781            ...     print "logarithm check failed for %s (diff = %s ulp)"% \
     1782            ...         (x, diff/log_RDF.ulp())
     1783            ...     return False
    17391784            sage: all( check_error(2^x) for x in range(-100,100) )
    17401785            True
    17411786            sage: all( check_error(x) for x in sxrange(0.01, 2.00, 0.01) )