Changeset 4122:dcc3a21b3bbc

Timestamp:

04/26/07 22:04:37 (6 years ago)

Author:

William Stein <wstein@…>

Branch:

default

Parents:

4120:cb21d372affb (diff), 4121:1058fc175758 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.

Message:

merge

Files:

• sage/rings/real_mpfr.pyx (modified) (1 diff)
• sage/rings/real_mpfr.pyx (modified) (17 diffs)

sage/rings/real_mpfr.pyx

@@ -1578 +1578 @@
             return intv.simplest_rational(low_open = True)
         
+    def nearby_rational(self, max_error=None, max_denominator=None):
+        """
+        Find a rational near to self.  Exactly one of max_error
+        or max_denominator must be specified.  If max_error is specified,
+        then this returns the simplest rational in the range
+        [self-max_error .. self+max_error].  If max_denominator is
+        specified, then this returns the rational closest to self with
+        denominator at most max_denominator.  (In case of ties, we
+        pick the simpler rational.)
+
+        EXAMPLES:
+            sage: (0.333).nearby_rational(max_error=0.001)
+            1/3
+            sage: (0.333).nearby_rational(max_error=1)
+            0
+            sage: (-0.333).nearby_rational(max_error=0.0001)
+            -257/772
+
+            sage: (0.333).nearby_rational(max_denominator=100)
+            1/3
+            sage: RR(1/3 + 1/1000000).nearby_rational(max_denominator=2999999)
+            777780/2333333
+            sage: RR(1/3 + 1/1000000).nearby_rational(max_denominator=3000000)
+            1000003/3000000
+            sage: (-0.333).nearby_rational(max_denominator=1000)
+            -333/1000
+            sage: RR(3/4).nearby_rational(max_denominator=2)
+            1
+            sage: RR(pi).nearby_rational(max_denominator=120)
+            355/113
+            sage: RR(pi).nearby_rational(max_denominator=10000)
+            355/113
+            sage: RR(pi).nearby_rational(max_denominator=100000)
+            312689/99532
+            sage: RR(pi).nearby_rational(max_denominator=1)
+            3
+            sage: RR(-3.5).nearby_rational(max_denominator=1)
+            -3
+        """
+        if ((max_error is None and max_denominator is None) or
+            (max_error is not None and max_denominator is not None)):
+            raise ValueError, 'Must specify exactly one of max_error or max_denominator in nearby_rational()'
+
+        if max_error is not None:
+            from real_mpfi import RealIntervalField
+
+            intv_field = RealIntervalField(self.prec())
+            intv = intv_field(self - max_error, self + max_error)
+
+            return intv.simplest_rational()
+
+        cdef int sgn = mpfr_sgn(self.value)
+
+        if sgn == 0:
+            return Rational(0)
+
+        cdef Rational self_r = self.exact_rational()
+
+        cdef Integer self_d = self_r.denominator()
+
+        if self_d <= max_denominator:
+            return self_r
+
+        if sgn < 0:
+            self_r = -self_r
+
+        cdef Integer fl = self_r.floor()
+        cdef Rational target = self_r - fl
+
+        cdef int low_done = 0
+        cdef int high_done = 0
+
+        # We use the Stern-Brocot tree to find the nearest neighbors of
+        # self with denominator at most max_denominator.  However,
+        # navigating the Stern-Brocot tree in the straightforward way
+        # can be very slow; for instance, to get to 1/1000000 takes a
+        # million steps.  Instead, we perform many steps at once;
+        # this probably slows down the average case, but it drastically
+        # speeds up the worst case.
+
+        # Suppose we have a/b < c/d < e/f, where a/b and e/f are
+        # neighbors in the Stern-Brocot tree and c/d is the target.
+        # We alternate between moving the low and the high end toward
+        # the target as much as possible.  Suppose that there are
+        # k consecutive rightward moves in the Stern-Brocot tree
+        # traversal; then we end up with (a+k*e)/(b+k*f).  We have
+        # two constraints on k.  First, the result must be <= c/d;
+        # this gives us the following:
+        # (a+k*e)/(b+k*f) <= c/d
+        # d*a + k*(d*e) <= c*b + k*(c*f)
+        # k*(d*e) - k*(c*f) <= c*b - d*a
+        # k <= (c*b - d*a)/(d*e - c*f)
+        # when moving the high side, we get
+        # (k*a+e)/(k*b+f) >= c/d
+        # k*(d*a) + d*e >= k*(c*b) + c*f
+        # d*e - c*f >= k*(c*b - d*a)
+        # k <= (d*e - c*f)/(c*b - d*a)
+
+        # We also need the denominator to be <= max_denominator; this
+        # gives (b+k*f) <= max_denominator or
+        # k <= (max_denominator - b)/f
+        # or
+        # k <= (max_denominator - f)/b
+
+        # We use variables with the same names as in the math above.
+
+        cdef Integer a = Integer(0)
+        cdef Integer b = Integer(1)
+        cdef Integer c = target.numerator()
+        cdef Integer d = target.denominator()
+        cdef Integer e = Integer(1)
+        cdef Integer f = Integer(1)
+
+        cdef Integer k
+
+        while (not low_done) or (not high_done):
+            # Move the low side
+            k = (c*b - d*a) // (d*e - c*f)
+
+            if b+k*f > max_denominator:
+                k = (max_denominator - b) // f
+                low_done = True
+
+            if k == 0:
+                low_done = True
+
+            a = a + k*e
+            b = b + k*f
+
+            # Move the high side
+            k = (d*e - c*f) // (c*b - d*a)
+
+            if k*b + f >= max_denominator:
+                k = (max_denominator - f) // b
+                high_done = True
+
+            if k == 0:
+                high_done = True
+
+            e = k*a + e
+            f = k*b + f
+
+        # Now a/b and e/f are rationals surrounding c/d.  We know that
+        # neither is equal to c/d, since d > max_denominator and
+        # b and f are both <= max_denominator.  (We know that
+        # d > max_denominator because we return early (before we
+        # get here) if d <= max_denominator.)
+
+        low = a/b
+        high = e/f
+
+        cdef int compare = cmp(target - low, high - target)
+
+        if compare > 0:
+            result = high
+        elif compare < 0:
+            result = low
+        else:
+            compare = cmp(b, f)
+            if compare > 0:
+                result = high
+            elif compare < 0:
+                result = low
+            else:
+                compare = cmp(a, e)
+                if compare > 0:
+                    result = high
+                else:
+                    result = low
+
+        result = fl + result
+
+        if sgn < 0:
+            return -result
+        return result        
+
 
     ###########################################

sage/rings/real_mpfr.pyx

@@ -48 +48 @@
 include '../ext/interrupt.pxi'
 include "../ext/stdsage.pxi"
+include "../ext/random.pxi"
 
 cimport sage.rings.ring
@@ -216 +217 @@
         """
         if hasattr(x, '_mpfr_'):
-            # This design with the hasattr is very annoying.
-            # The only thing that uses it right now is symbolic constants
-            # and symbolic function evaluation.
-            # Getting rid of this would speed things up.
             return x._mpfr_(self)
         cdef RealNumber z
@@ -246 +243 @@
         elif self.__prec <= 53 and is_RealDoubleElement(x):
             return self(x)
-        import sage.functions.constants
-        return self._coerce_try(x, [sage.functions.constants.ConstantRing])
+        raise TypeError
+        #import sage.functions.constants
+        #return self._coerce_try(x, [sage.functions.constants.ConstantRing])
 
     def __cmp__(self, other):
@@ -411 +409 @@
             0.69314718055994530941723212146
         """
-        cdef RealNumber x
-        x = self._new()        
+        cdef RealNumber x = self._new()
         mpfr_const_log2(x.value, self.rnd)
         return x
+        
+    def random_element(self, min=-1, max=1, distribution=None):
+        """
+        Returns a uniformly distributed random number between 
+        min and max (default -1 to 1). 
+        
+        EXAMPLES:
+            sage: RealField(100).random_element(-5, 10)
+            4.2682457657074627882421620493
+            sage: RealField(10).random_element()
+            .27
+        """
+        cdef RealNumber x = self._new()
+        mpfr_urandomb(x.value, state)
+        if min == 0 and max == 1:
+            return x
+        else:
+            return (max-min)*x + min
 
     def factorial(self, int n):
@@ -663 +678 @@
 
     def _latex_(self):
-        return str(self)
+        s = self.str()
+        parts = s.split('e')
+        if len(parts) > 1:
+            # scientific notation
+            if parts[1][0] == '+':
+                parts[1] = parts[1][1:]
+            s = "%s \\times 10^{%s}" % (parts[0], parts[1])
+        return s
 
     def _interface_init_(self):
@@ -1263 +1285 @@
         Returns integer truncation of this real number.
         """
-        s = self.str(32)
-        i = s.find('.')
-        return int(s[:i], 32)
+        cdef Integer z = Integer()
+        mpfr_get_z(z.value, self.value, GMP_RNDZ)
+        return z.__int__()
 
     def __long__(self):
@@ -1271 +1293 @@
         Returns long integer truncation of this real number.
         """
-        s = self.str(32)
-        i = s.find('.')
-        return long(s[:i], 32)
+        cdef Integer z = Integer()
+        mpfr_get_z(z.value, self.value, GMP_RNDZ)
+        return z.__long__()
 
     def __complex__(self):
@@ -1816 +1838 @@
 
     def sqrt(self):
-        """
+        r"""
         Return a square root of self.
 
         If self is negative a complex number is returned.
 
-        If you use self.square_root() then a real number will always
-        be returned (though it will be NaN if self is negative).
+        If you use \code{self.sqrt_approx()} then a real number will
+        always be returned (though it will be NaN if self is
+        negative).
 
         EXAMPLES:
@@ -1833 +1856 @@
             sage: r = 4344
             sage: r.sqrt()
-            65.9090282131363
+            2*sqrt(1086)
+            
+            sage: r = 4344.0
             sage: r.sqrt()^2 == r
             True
             sage: r.sqrt()^2 - r            
-             0.000000000000000
+            0.000000000000000
 
             sage: r = -2.0
@@ -1844 +1869 @@
             """
         if self >= 0:
-            return self.square_root()
+            return self.sqrt_approx()
         return self._complex_number_().sqrt()
-        
-
-    def square_root(self):
+
+    def sqrt_approx(self):
         """
         Return a square root of self.  A real number will always be
@@ -1857 +1881 @@
         EXAMPLES:
             sage: r = -2.0
-            sage: r.square_root()
+            sage: r.sqrt_approx()
             NaN
             sage: r.sqrt()
@@ -1918 +1942 @@
             1.0000000*I                    # 32-bit
             -1.0842022e-19 + 1.0000000*I   # 64-bit
+
+        We raise a real number to a symbolic object:
+            sage: 1.5^x
+            1.50000000000000^x
+            sage: -2.3^(x+y^3+sin(x))
+            -2.30000000000000^(y^3 + sin(x) + x)
         """
         cdef RealNumber x
@@ -1923 +1953 @@
             return self.__pow__(float(exponent))
         if not PY_TYPE_CHECK(exponent, RealNumber):
-            x = self
-            exponent = x._parent(exponent)
+            try:
+                x = self
+                exponent = x._parent(exponent)
+            except TypeError:
+                try:
+                    return exponent.parent()(self)**exponent
+                except AttributeError:
+                    raise TypeError
         return self.__pow(exponent)
 
@@ -2259 +2295 @@
         return x
 
+    def coth(self):
+        """
+        EXAMPLES:
+            sage: RealField(100)(2).coth()
+            1.0373147207275480958778097648
+        """
+        return 1/self.tanh()
+
+    def csch(self):
+        """
+        EXAMPLES:
+            sage: RealField(100)(2).csch()
+            0.27572056477178320775835148216
+        """
+        return 1/self.sinh()
+
+    def sech(self):
+        """
+        EXAMPLES:
+            sage: RealField(100)(2).sech()
+            0.26580222883407969212086273982        
+        """
+        return 1/self.cosh()
+
     def acosh(self):
         """
@@ -2419 +2479 @@
 
         EXAMPLE:
-             sage: r = sqrt(2); r
+             sage: r = sqrt(2.0); r
              1.41421356237310
              sage: r.algdep(5)
@@ -2436 +2496 @@
 
          EXAMPLE:
-              sage: r = sqrt(2); r
+              sage: r = sqrt(2.0); r
               1.41421356237310
               sage: r.algdep(5)
@@ -2536 +2596 @@
 
 def is_RealField(x):
-    return PY_TYPE_CHECK(x, RealField)
+    return bool(PY_TYPE_CHECK(x, RealField))
 
 def is_RealNumber(x):
-    return PY_TYPE_CHECK(x, RealNumber)
+    return bool(PY_TYPE_CHECK(x, RealNumber))
 
 def __create__RealField_version0(prec, sci_not, rnd):