CoCalc Public Filestmp / 2015-03-19-092125.sagewsOpen in with one click!
Authors: Harald Schilly, ℏal Snyder, William A. Stein
1/5 in RIF
False
a = RIF(1/5)
a
0.2000000000000000?
a == 1/5
False
RIF(1/5) == RIF(1/5)
False
a
0.2000000000000000?
a == a
False
a is a
True
a.lower() == a.lower()
True
a.upper() == a.upper()
True
a??
File: /usr/local/sage/sage-6.5/src/sage/rings/real_mpfi.pyx Source: cdef class RealIntervalFieldElement(sage.structure.element.RingElement): """ A real number interval. """ cdef RealIntervalFieldElement _new(self): """ Return a new real interval with same parent as ``self``. """ cdef RealIntervalFieldElement x x = PY_NEW(RealIntervalFieldElement) x._parent = self._parent mpfi_init2(x.value, (<RealIntervalField_class>self._parent).__prec) x.init = 1 return x def __init__(self, RealIntervalField_class parent, x=0, int base=10): """ Create a real interval element. Should be called by first creating a :class:`RealIntervalField`, as illustrated in the examples. EXAMPLES:: sage: R = RealIntervalField() sage: R('1.2456') 1.245600000000000? sage: R = RealIntervalField(3) sage: R('1.2456').str(style='brackets') '[1.0 .. 1.3]' EXAMPLES: Rounding:: sage: w = RealIntervalField(3)(5/2) sage: RealIntervalField(2)(w).str(2, style='brackets') '[10. .. 11.]' TESTS:: sage: a = RealIntervalField(428)(factorial(100)/exp(2)); a 1.26303298005073195998439505058085204028142920134742241494671502106333548593576383141666758300089860337889002385197008191910406895?e157 sage: a.diameter() 4.7046373946079775711568954992429894854882556641460240333441655212438503516287848720594584761250430179569094634219773739322602945e-129 Type: ``RealIntervalField?`` for many more examples. """ import sage.rings.qqbar self.init = 0 if parent is None: raise TypeError self._parent = parent mpfi_init2(self.value, parent.__prec) self.init = 1 if x is None: return cdef RealIntervalFieldElement _x, n, d cdef RealNumber rn, rn1 cdef Rational rat, rat1 cdef Integer integ, integ1 cdef RealDoubleElement dx, dx1 cdef int ix, ix1 if PY_TYPE_CHECK(x, RealIntervalFieldElement): _x = x # so we can get at x.value mpfi_set(self.value, _x.value) elif PY_TYPE_CHECK(x, RealNumber): rn = x mpfi_set_fr(self.value, rn.value) elif PY_TYPE_CHECK(x, Rational): rat = x mpfi_set_q(self.value, rat.value) elif PY_TYPE_CHECK(x, Integer): integ = x mpfi_set_z(self.value, integ.value) elif PY_TYPE_CHECK(x, int): ix = x mpfi_set_si(self.value, ix) elif isinstance(x, tuple): try: a, b = x except ValueError: raise TypeError("tuple defining an interval must have length 2") if PY_TYPE_CHECK(a, RealNumber) and PY_TYPE_CHECK(b, RealNumber): rn = a rn1 = b mpfi_interv_fr(self.value, rn.value, rn1.value) elif PY_TYPE_CHECK(a, RealDoubleElement) and PY_TYPE_CHECK(b, RealDoubleElement): dx = a dx1 = b mpfi_interv_d(self.value, dx._value, dx1._value) elif PY_TYPE_CHECK(a, Rational) and PY_TYPE_CHECK(b, Rational): rat = a rat1 = b mpfi_interv_q(self.value, rat.value, rat1.value) elif PY_TYPE_CHECK(a, Integer) and PY_TYPE_CHECK(b, Integer): integ = a integ1 = b mpfi_interv_z(self.value, integ.value, integ1.value) elif PY_TYPE_CHECK(a, int) and PY_TYPE_CHECK(b, int): ix = a ix1 = b mpfi_interv_si(self.value, ix, ix1) else: # generic fallback rn = self._parent(a).lower() rn1 = self._parent(b).upper() mpfi_interv_fr(self.value, rn.value, rn1.value) elif isinstance(x, sage.rings.qqbar.AlgebraicReal): d = x.interval(self._parent) mpfi_set(self.value, d.value) elif hasattr(x, '_real_mpfi_'): d = x._real_mpfi_(self._parent) mpfi_set(self.value, d.value) else: from sage.symbolic.expression import Expression if isinstance(x, Expression): d = x._real_mpfi_(self._parent) mpfi_set(self.value, d.value) elif isinstance(x, str): # string s = str(x).replace('..', ',').replace(' ','').replace('+infinity', '@[email protected]').replace('-infinity','[email protected]@') if mpfi_set_str(self.value, s, base): raise TypeError("Unable to convert number to real interval.") else: # try coercing to real try: rn = self._parent._lower_field()(x) rn1 = self._parent._upper_field()(x) except TypeError: raise TypeError("Unable to convert number to real interval.") mpfi_interv_fr(self.value, rn.value, rn1.value) def __reduce__(self): """ Pickling support. EXAMPLES:: sage: a = RIF(5,5.5) sage: cmp(loads(dumps(a)), a) 0 sage: R = RealIntervalField(sci_not=1, prec=200) sage: b = R('393.39203845902384098234098230948209384028340') sage: cmp(loads(dumps(b)), b) 0 sage: b = R(1)/R(0); b # R(0) has no particular sign, thus 1/R(0) covers the whole reals [-infinity .. +infinity] sage: c = loads(dumps(b)) sage: (c.lower(), c.upper()) == (b.lower(), b.upper()) True sage: b = R(-1)/R(0); b # same as above [-infinity .. +infinity] sage: c = loads(dumps(b)) sage: (c.lower(), c.upper()) == (b.lower(), b.upper()) True sage: b = R('[2 .. 3]'); b.str(error_digits=1) '2.5?5e0' sage: cmp(loads(dumps(b)), b) 0 sage: R = RealIntervalField(4000) sage: s = 1/R(3) sage: t = loads(dumps(s)) sage: (t.upper(), t.lower()) == (s.upper(), s.lower()) True sage: loads(dumps(1/RIF(0,1))) [1.0000000000000000 .. +infinity] """ return (__create__RealIntervalFieldElement_version1, (self._parent, self.upper(), self.lower())) def __dealloc__(self): """ Deallocate ``self``. EXAMPLES:: sage: R = RealIntervalField() sage: del R # indirect doctest """ if self.init: mpfi_clear(self.value) def __repr__(self): """ Return a string representation of ``self``. EXAMPLES:: sage: R = RealIntervalField() sage: R(1.2456) # indirect doctest 1.2456000000000001? """ return self.str(10) def _latex_(self): """ Return a latex represention of ``self``. EXAMPLES:: sage: latex(RIF(1.5)) # indirect doctest 1.5000000000000000? sage: latex(RIF(3e200)) # indirect doctest 3.0000000000000000? \times 10^{200} """ import re return re.sub(r"e(-?\d+)", r" \\times 10^{\1}", str(self)) def _magma_init_(self, magma): r""" Return a string representation of ``self`` in the Magma language. EXAMPLES:: sage: t = RIF(10, 10.5); t 11.? sage: magma(t) # optional - magma # indirect doctest 10.2500000000000 """ return "%s!%s" % (self.parent()._magma_init_(magma), self.center()) def _interface_init_(self, I=None): """ Raise a ``TypeError``. This function would return the string representation of ``self`` that makes sense as a default representation of a real interval in other computer algebra systems. But, most other computer algebra systems do not support interval arithmetic, so instead we just raise a ``TypeError``. Define the appropriate ``_cas_init_`` function if there is a computer algebra system you would like to support. EXAMPLES:: sage: n = RIF(1.3939494594) sage: n._interface_init_() Traceback (most recent call last): ... TypeError Here's a conversion to Maxima happens, which results in a type error:: sage: a = RealInterval('2.3') sage: maxima(a) Traceback (most recent call last): ... TypeError """ raise TypeError def _sage_input_(self, sib, coerce): r""" Produce an expression which will reproduce this value when evaluated. EXAMPLES:: sage: sage_input(RIF(e, pi), verify=True) # Verified RIF(RR(2.7182818284590451), RR(3.1415926535897936)) sage: sage_input(RealIntervalField(64)(sqrt(2)), preparse=False, verify=True) # Verified RR64 = RealField(64) RealIntervalField(64)(RR64('1.41421356237309504876'), RR64('1.41421356237309504887')) sage: sage_input(RealIntervalField(2)(12), verify=True) # Verified RealIntervalField(2)(RealField(2)(12.)) sage: sage_input(RealIntervalField(2)(13), verify=True) # Verified RR2 = RealField(2) RealIntervalField(2)(RR2(12.), RR2(16.)) sage: from sage.misc.sage_input import SageInputBuilder sage: sib = SageInputBuilder() sage: RIF(-sqrt(3), -sqrt(2))._sage_input_(sib, False) {call: {atomic:RIF}({unop:- {call: {atomic:RR}({atomic:1.7320508075688774})}}, {unop:- {call: {atomic:RR}({atomic:1.4142135623730949})}})} """ # Interval printing could often be much prettier, # but I'm feeling lazy :) if self.is_exact(): return sib(self.parent())(sib(self.lower(rnd='RNDN'))) else: # The following line would also be correct, but even though it # uses coerced=2, that doesn't help because RealNumber doesn't # print pretty for directed-rounding fields. #return sib(self.parent())(sib(self.lower(), 2), sib(self.upper(), 2)) return sib(self.parent())(sib(self.lower(rnd='RNDN')), sib(self.upper(rnd='RNDN'))) def __hash__(self): """ Return a hash value of ``self``. EXAMPLES:: sage: hash(RIF(e)) == hash(RIF(e)) # indirect doctest True """ return hash(self.str(16)) def _im_gens_(self, codomain, im_gens): """ Return the image of ``self`` under the homomorphism from the rational field to ``codomain``. This always just returns ``self`` coerced into the ``codomain``. EXAMPLES:: sage: RIF(2.1)._im_gens_(CIF, [CIF(1)]) 2.1000000000000001? sage: R = RealIntervalField(20) sage: RIF(2.1)._im_gens_(R, [R(1)]) 2.10000? """ return codomain(self) # since 1 |--> 1 def real(self): """ Return the real part of ``self``. (Since ``self`` is a real number, this simply returns ``self``.) EXAMPLES:: sage: RIF(1.2465).real() == RIF(1.2465) True """ return self # MPFR had an elaborate "truncation" scheme to avoid printing # inaccurate-looking results; this has been removed for MPFI, # because I think it's less confusing to get such results than to # see [0.333 .. 0.333] for an interval with unequal left and right # sides. def str(self, int base=10, style=None, no_sci=None, e=None, error_digits=None): r""" Return a string representation of ``self``. INPUT: - ``base`` -- base for output - ``style`` -- The printing style; either ``'brackets'`` or ``'question'`` (or ``None``, to use the current default). - ``no_sci`` -- if ``True`` do not print using scientific notation; if ``False`` print with scientific notation; if ``None`` (the default), print how the parent prints. - ``e`` -- symbol used in scientific notation - ``error_digits`` -- The number of digits of error to print, in ``'question'`` style. We support two different styles of printing; ``'question'`` style and ``'brackets'`` style. In question style (the default), we print the "known correct" part of the number, followed by a question mark:: sage: RIF(pi).str() '3.141592653589794?' sage: RIF(pi, 22/7).str() '3.142?' sage: RIF(pi, 22/7).str(style='question') '3.142?' However, if the interval is precisely equal to some integer that's not too large, we just return that integer:: sage: RIF(-42).str() '-42' sage: RIF(0).str() '0' sage: RIF(12^5).str(base=3) '110122100000' Very large integers, however, revert to the normal question-style printing:: sage: RIF(3^7).str() '2187' sage: RIF(3^7 * 2^256).str() '2.5323729916201052?e80' In brackets style, we print the lower and upper bounds of the interval within brackets:: sage: RIF(237/16).str(style='brackets') '[14.812500000000000 .. 14.812500000000000]' Note that the lower bound is rounded down, and the upper bound is rounded up. So even if the lower and upper bounds are equal, they may print differently. (This is done so that the printed representation of the interval contains all the numbers in the internal binary interval.) For instance, we find the best 10-bit floating point representation of ``1/3``:: sage: RR10 = RealField(10) sage: RR(RR10(1/3)) 0.333496093750000 And we see that the point interval containing only this floating-point number prints as a wider decimal interval, that does contain the number:: sage: RIF10 = RealIntervalField(10) sage: RIF10(RR10(1/3)).str(style='brackets') '[0.33349 .. 0.33350]' We always use brackets style for ``NaN`` and infinities:: sage: RIF(pi, infinity) [3.1415926535897931 .. +infinity] sage: RIF(NaN) [.. NaN ..] Let's take a closer, formal look at the question style. In its full generality, a number printed in the question style looks like: MANTISSA ?ERROR eEXPONENT (without the spaces). The "eEXPONENT" part is optional; if it is missing, then the exponent is 0. (If the base is greater than 10, then the exponent separator is "@" instead of "e".) The "ERROR" is optional; if it is missing, then the error is 1. The mantissa is printed in base `b`, and always contains a decimal point (also known as a radix point, in bases other than 10). (The error and exponent are always printed in base 10.) We define the "precision" of a floating-point printed representation to be the positional value of the last digit of the mantissa. For instance, in ``2.7?e5``, the precision is `10^4`; in ``8.?``, the precision is `10^0`; and in ``9.35?`` the precision is `10^{-2}`. This precision will always be `10^k` for some `k` (or, for an arbitrary base `b`, `b^k`). Then the interval is contained in the interval: .. MATH:: \text{mantissa} \cdot b^{\text{exponent}} - \text{error} \cdot b^k .. \text{mantissa} \cdot b^{\text{exponent}} + \text{error} \cdot b^k To control the printing, we can specify a maximum number of error digits. The default is 0, which means that we do not print an error at all (so that the error is always the default, 1). Now, consider the precisions needed to represent the endpoints (this is the precision that would be produced by ``v.lower().str(no_sci=False, truncate=False)``). Our result is no more precise than the less precise endpoint, and is sufficiently imprecise that the error can be represented with the given number of decimal digits. Our result is the most precise possible result, given these restrictions. When there are two possible results of equal precision and with the same error width, then we pick the one which is farther from zero. (For instance, ``RIF(0, 123)`` with two error digits could print as ``61.?62`` or ``62.?62``. We prefer the latter because it makes it clear that the interval is known not to be negative.) EXAMPLES:: sage: a = RIF(59/27); a 2.185185185185186? sage: a.str() '2.185185185185186?' sage: a.str(style='brackets') '[2.1851851851851851 .. 2.1851851851851856]' sage: a.str(16) '2.2f684bda12f69?' sage: a.str(no_sci=False) '2.185185185185186?e0' sage: pi_appr = RIF(pi, 22/7) sage: pi_appr.str(style='brackets') '[3.1415926535897931 .. 3.1428571428571433]' sage: pi_appr.str() '3.142?' sage: pi_appr.str(error_digits=1) '3.1422?7' sage: pi_appr.str(error_digits=2) '3.14223?64' sage: pi_appr.str(base=36) '3.6?' sage: RIF(NaN) [.. NaN ..] sage: RIF(pi, infinity) [3.1415926535897931 .. +infinity] sage: RIF(-infinity, pi) [-infinity .. 3.1415926535897936] sage: RealIntervalField(210)(3).sqrt() 1.732050807568877293527446341505872366942805253810380628055806980? sage: RealIntervalField(210)(RIF(3).sqrt()) 1.732050807568878? sage: RIF(3).sqrt() 1.732050807568878? sage: RIF(0, 3^-150) 1.?e-71 TESTS: Check that :trac:`13634` is fixed:: sage: RIF(0.025) 0.025000000000000002? sage: RIF.scientific_notation(True) sage: RIF(0.025) 2.5000000000000002?e-2 sage: RIF.scientific_notation(False) sage: RIF(0.025) 0.025000000000000002? """ if base < 2 or base > 36: raise ValueError, "the base (=%s) must be between 2 and 36"%base # If self is a NaN, always use brackets style. if mpfi_nan_p(self.value): if base >= 24: return "[.. @[email protected] ..]" else: return "[.. NaN ..]" if e is None: if base > 10: e = '@' else: e = 'e' if style is None: style = printing_style if error_digits is None: error_digits = printing_error_digits if mpfr_inf_p(&self.value.left) or mpfr_inf_p(&self.value.right): style = 'brackets' if style == 'brackets': t1 = self.lower().str(base=base, no_sci=no_sci, e=e, truncate=False) t2 = self.upper().str(base=base, no_sci=no_sci, e=e, truncate=False) return "[%s .. %s]"%(t1, t2) elif style == 'question': if no_sci is None: prefer_sci = self.parent().scientific_notation() elif not no_sci: prefer_sci = True else: prefer_sci = False return self._str_question_style(base, error_digits, e, prefer_sci) else: raise ValueError, 'Illegal interval printing style %s'%printing_style cpdef _str_question_style(self, int base, int error_digits, e, bint prefer_sci): r""" Compute the "question-style print representation" of this value, with the given base and error_digits. See the documentation for the str method for the definition of the question style. INPUTS: - ``base`` - base for output - ``error_digits`` - maximum number of decimal digits for error - ``e`` - symbol for exponent (typically ``'e'`` for base less than or equal to 10, ``'@'`` for larger base) - ``prefer_sci`` - ``True`` to always print in scientific notation; ``False`` to prefer non-scientific notation when possible EXAMPLES:: sage: v = RIF(e, pi) sage: v.str(2, style='brackets') '[10.101101111110000101010001011000101000101011101101001 .. 11.001001000011111101101010100010001000010110100011001]' sage: v._str_question_style(2, 0, 'e', False) '11.0?' sage: v._str_question_style(2, 3, '@', False) '10.111011100001?867' sage: v._str_question_style(2, 30, 'e', False) '10.111011100001000001011101111101011000100001001000001?476605618580184' sage: v.str(5, style='brackets') '[2.32434303404423034041024 .. 3.03232214303343241124132]' sage: v._str_question_style(5, 3, '@', False) '2.43111?662' sage: (Integer('232434', 5) + 2*662).str(5) '303233' sage: v.str(style='brackets') '[2.7182818284590450 .. 3.1415926535897936]' sage: v._str_question_style(10, 3, '@', False) '2.930?212' sage: v._str_question_style(10, 3, '@', True) '[email protected]' sage: v.str(16, style='brackets') '[2.b7e151628aed2 .. 3.243f6a8885a32]' sage: v._str_question_style(16, 3, '@', False) '2.ee1?867' sage: (Integer('2b7e', 16) + 2*867).str(16) '3244' sage: v.str(36, style='brackets') '[2.puw5nggjf8f .. 3.53i5ab8p5gz]' sage: v._str_question_style(36, 3, '@', False) '2.xh?275' sage: (Integer('2pu', 36) + 2*275).str(36) '354' TESTS:: sage: RIF(0, infinity)._str_question_style(10, 0, 'e', False) Traceback (most recent call last): ... ValueError: _str_question_style on NaN or infinity sage: for i in range(1, 9): ... print RIF(-10^i, 12345).str(style='brackets') ... print RIF(-10^i, 12345)._str_question_style(10, 0, 'e', False) ... print RIF(-10^i, 12345)._str_question_style(10, 3, 'e', False) [-10.000000000000000 .. 12345.000000000000] 0.?e5 6.17?618e3 [-100.00000000000000 .. 12345.000000000000] 0.?e5 6.13?623e3 [-1000.0000000000000 .. 12345.000000000000] 0.?e5 5.68?668e3 [-10000.000000000000 .. 12345.000000000000] 0.?e5 1.2?112e3 [-100000.00000000000 .. 12345.000000000000] 0.?e5 -4.38?562e4 [-1.0000000000000000e6 .. 12345.000000000000] 0.?e6 -4.94?507e5 [-1.0000000000000000e7 .. 12345.000000000000] 0.?e7 -4.99?501e6 [-1.0000000000000000e8 .. 12345.000000000000] 0.?e8 -5.00?501e7 sage: RIF(10^-3, 12345)._str_question_style(10, 3, 'e', False) '6.18?618e3' sage: RIF(-golden_ratio, -10^-6)._str_question_style(10, 3, 'e', False) '-0.810?810' sage: RIF(-0.85, 0.85)._str_question_style(10, 0, 'e', False) '0.?' sage: RIF(-0.85, 0.85)._str_question_style(10, 1, 'e', False) '0.0?9' sage: RIF(-0.85, 0.85)._str_question_style(10, 2, 'e', False) '0.00?85' sage: RIF(-8.5, 8.5)._str_question_style(10, 0, 'e', False) '0.?e1' sage: RIF(-85, 85)._str_question_style(10, 0, 'e', False) '0.?e2' sage: for i in range(-6, 7): ... v = RIF(2).sqrt() * 10^i ... print v._str_question_style(10, 0, 'e', False) 1.414213562373095?e-6 0.00001414213562373095? 0.0001414213562373095? 0.001414213562373095? 0.01414213562373095? 0.1414213562373095? 1.414213562373095? 14.14213562373095? 141.4213562373095? 1414.213562373095? 14142.13562373095? 141421.3562373095? 1.414213562373095?e6 sage: RIF(3^33) 5559060566555523 sage: RIF(3^33 * 2) 1.1118121133111046?e16 sage: RIF(-pi^-512, 0) -1.?e-254 sage: RealIntervalField(2)(3 * 2^18) 786432 sage: RealIntervalField(2)(3 * 2^19) 1.6?e6 sage: v = RIF(AA(2-sqrt(3))) sage: v 0.2679491924311227? sage: v.str(error_digits=1) '0.26794919243112273?4' sage: v.str(style='brackets') '[0.26794919243112269 .. 0.26794919243112276]' sage: -v -0.2679491924311227? Check that :trac:`15166` is fixed:: sage: RIF(1.84e13).exp() [2.0985787164673874e323228496 .. +infinity] # 32-bit 6.817557048799520?e7991018467019 # 64-bit sage: from sage.rings.real_mpfr import mpfr_get_exp_min, mpfr_get_exp_max sage: v = RIF(1.0 << (mpfr_get_exp_max() - 1)); v 1.0492893582336939?e323228496 # 32-bit 2.9378268945557938?e1388255822130839282 # 64-bit sage: -v -1.0492893582336939?e323228496 # 32-bit -2.9378268945557938?e1388255822130839282 # 64-bit sage: v = RIF(1.0 >> -mpfr_get_exp_min()+1); v 2.3825649048879511?e-323228497 # 32-bit 8.5096913117408362?e-1388255822130839284 # 64-bit """ if not(mpfr_number_p(&self.value.left) and mpfr_number_p(&self.value.right)): raise ValueError, "_str_question_style on NaN or infinity" if base < 2 or base > 36: raise ValueError, "the base (=%s) must be between 2 and 36"%base if error_digits < 0 or error_digits > 1000: # The restriction to 1000 is not essential. The reason to have # a restriction is that this code is not efficient for # large error_digits values (for instance, we always construct # a number with that many digits); a very large error_digits # could run out of memory, etc. # I think that 1000 is a pretty "safe" limit. The whole # purpose of question_style is to be human-readable, and # the human-readability will go way down after about 6 # error digits; 1000 error digits is just silly. raise ValueError, "error_digits (=%s) must be between 0 and 1000"%error_digits cdef mp_exp_t self_exp cdef mpz_t self_zz cdef int prec = (<RealIntervalField_class>self._parent).__prec cdef char *zz_str cdef size_t zz_str_maxlen if mpfr_equal_p(&self.value.left, &self.value.right) \ and mpfr_integer_p(&self.value.left): # This might be suitable for integer printing, but not if it's # too big. (We can represent 2^3000000 exactly in RIF, but we # don't want to print this 903090 digit number; we'd rather # just print 9.7049196389007116?e903089 .) # Represent self as m*2^k, where m is an integer with # self.prec() bits and k is an integer. (So RIF(1) would have # m = 2^52 and k=-52.) Then, as a simple heuristic, we print # as an integer if k<=0. (As a special dispensation for tiny # precisions, we also print as an integer if the number is # less than a million (actually, less than 2^20); this # will never affect "normal" uses, but it makes tiny examples # with RealIntervalField(2) prettier.) self_exp = mpfr_get_exp(&self.value.left) if mpfr_zero_p(&self.value.left) or self_exp <= prec or self_exp <= 20: mpz_init(self_zz) mpfr_get_z(self_zz, &self.value.left, GMP_RNDN) zz_str_maxlen = mpz_sizeinbase(self_zz, base) + 2 zz_str = <char *>PyMem_Malloc(zz_str_maxlen) if zz_str == NULL: mpz_clear(self_zz) raise MemoryError, "Unable to allocate memory for integer representation of interval" sig_on() mpz_get_str(zz_str, base, self_zz) sig_off() v = PyString_FromString(zz_str) PyMem_Free(zz_str) return v # We want the endpoints represented as an integer mantissa # and an exponent, using the given base. MPFR will do that for # us in mpfr_get_str, so we end up converting from MPFR to strings # to MPZ to strings... ouch. We could avoid this overhead by # copying most of the guts of mpfr_get_str to give something like # mpfr_get_mpz_exp (to give the mantissa as an mpz and the exponent); # we could also write the body of this method using the string # values, so that the second and third conversions are always # on small values (of around the size of error_digits). # We don't do any of that. cdef char *lower_s cdef char *upper_s cdef mp_exp_t lower_expo cdef mp_exp_t upper_expo cdef mpz_t lower_mpz cdef mpz_t upper_mpz sig_on() lower_s = mpfr_get_str(<char*>0, &lower_expo, base, 0, &self.value.left, GMP_RNDD) upper_s = mpfr_get_str(<char*>0, &upper_expo, base, 0, &self.value.right, GMP_RNDU) sig_off() if lower_s == <char*> 0: raise RuntimeError, "Unable to convert interval lower bound to a string" if upper_s == <char*> 0: raise RuntimeError, "Unable to convert interval upper bound to a string" # MPFR returns an exponent assuming that the implicit radix point # is to the left of the first mantissa digit. We'll be doing # arithmetic on mantissas that might not preserve the number of # digits; we adjust the exponent so that the radix point is to # the right of the last mantissa digit (that is, so the number # is mantissa*base^exponent, if you interpret mantissa as an integer). cdef long digits digits = strlen(lower_s) if lower_s[0] == '-': digits -= 1 lower_expo -= digits digits = strlen(upper_s) if upper_s[0] == '-': digits -= 1 upper_expo -= digits sig_on() mpz_init_set_str(lower_mpz, lower_s, base) mpz_init_set_str(upper_mpz, upper_s, base) mpfr_free_str(lower_s) mpfr_free_str(upper_s) sig_off() cdef mpz_t tmp mpz_init(tmp) # At several places in the function, we divide by a power of # base. This could be sped up by shifting instead, when base # is a power of 2. (I'm not bothering right now because I # expect the not-base-10 case to be quite rare, especially # when combined with high precision and question style.) # First we normalize so that both mantissas are at the same # precision. This just involves dividing the more-precise # endpoint by the appropriate power of base. However, there's # one potentially very slow case we want to avoid: if the two # exponents are sufficiently different, the simple code would # involve computing a huge power of base, and then dividing # by it to get -1, 0, or 1 (depending on the rounding). # There's one complication first: if one of the endpoints is zero, # we want to treat it as infinitely precise (otherwise, it defaults # to a precision of 2^-self.prec(), so that RIF(0, 2^-1000) # would print as 1.?e-17). (If both endpoints are zero, then # we can't get here; we already returned '0' in the integer # case above.) if mpfr_zero_p(&self.value.left): lower_expo = upper_expo if mpfr_zero_p(&self.value.right): upper_expo = lower_expo cdef mp_exp_t expo_delta if lower_expo < upper_expo: expo_delta = upper_expo - lower_expo if mpz_sizeinbase(lower_mpz, base) < expo_delta: # abs(lower) < base^expo_delta, so # floor(lower/base^expo_delta) is either -1 or 0 # (depending on the sign of lower). if mpz_sgn(lower_mpz) < 0: mpz_set_si(lower_mpz, -1) else: mpz_set_ui(lower_mpz, 0) else: mpz_ui_pow_ui(tmp, base, expo_delta) mpz_fdiv_q(lower_mpz, lower_mpz, tmp) lower_expo = upper_expo elif upper_expo < lower_expo: expo_delta = lower_expo - upper_expo if mpz_sizeinbase(upper_mpz, base) < expo_delta: # abs(upper) < base^expo_delta, so # ceiling(upper/base^expo_delta) is either 0 or 1 # (depending on the sign of upper). if mpz_sgn(upper_mpz) > 0: mpz_set_ui(upper_mpz, 1) else: mpz_set_ui(upper_mpz, 0) else: mpz_ui_pow_ui(tmp, base, expo_delta) mpz_cdiv_q(upper_mpz, upper_mpz, tmp) upper_expo = lower_expo cdef mp_exp_t expo = lower_expo # Now the basic plan is to repeat # lower = floor(lower/base); upper = ceiling(upper/base) # until the error upper-lower is sufficiently small. However, # if this loop executes many times, that's pretty inefficient # (quadratic). We do a single pre-check to see approximately # how many times we would have to go through that loop, # and do most of them with a single division. cdef mpz_t cur_error mpz_init(cur_error) mpz_sub(cur_error, upper_mpz, lower_mpz) cdef mpz_t max_error mpz_init(max_error) if error_digits == 0: mpz_set_ui(max_error, 2) else: mpz_ui_pow_ui(max_error, 10, error_digits) mpz_sub_ui(max_error, max_error, 1) mpz_mul_2exp(max_error, max_error, 1) # Now we want to compute k as large as possible such that # ceiling(upper/base^(k-1)) - floor(lower/base^(k-1)) > max_error. # We start by noting that # ceiling(upper/base^(k-1)) - floor(lower/base^(k-1)) >= # cur_error/base^(k-1), # so it suffices if # cur_error/base^(k-1) > max_error, or # cur_error/max_error > base^(k-1). # We would like to take logarithms and subtract, but taking # logarithms is expensive. Instead we use mpz_sizeinbase # as an approximate logarithm. (This could probably be # improved, either by assuming undocumented knowledge of # the internals of mpz_sizeinbase, or by writing our own # approximate logarithm.) cdef long cur_error_digits = mpz_sizeinbase(cur_error, base) cdef long max_error_digits = mpz_sizeinbase(max_error, base) # The GMP documentation claims that mpz_sizeinbase will be either # the true number of digits, or one too high. So # cur_error might have as few as cur_error_digits-1 digits, # so it might be as small as base^(cur_error_digits-2). # max_error might have as many as max_error_digits digits, so it # might be almost (but not quite) as large as base^max_error_digits. # Then their quotient will be at least slightly larger than # base^(cur_error_digits-2-max_error_digits). So we can take # k-1 = cur_error_digits-2-max_error_digits, and # k = cur_error_digits-1-max_error_digits. cdef long k = cur_error_digits - 1 - max_error_digits if k > 0: mpz_ui_pow_ui(tmp, base, k) mpz_fdiv_q(lower_mpz, lower_mpz, tmp) mpz_cdiv_q(upper_mpz, upper_mpz, tmp) expo += k mpz_sub(cur_error, upper_mpz, lower_mpz) # OK, we've almost divided enough times to fit within max_error. # (In fact, maybe we already have.) Now we just loop a few more # times until we're done. while mpz_cmp(cur_error, max_error) > 0: mpz_fdiv_q_ui(lower_mpz, lower_mpz, base) mpz_cdiv_q_ui(upper_mpz, upper_mpz, base) expo += 1 mpz_sub(cur_error, upper_mpz, lower_mpz) # Almost done. Now we need to print out a floating-point number # with a mantissa halfway between lower_mpz and upper_mpz, # an error of half of cur_error (rounded up), and an exponent # based on expo (shifted by the location of the decimal point # within the mantissa). # We briefly re-purpose lower_mpz to hold the final mantissa: mpz_add(lower_mpz, lower_mpz, upper_mpz) # According to our spec, we're supposed to divide lower_mpz # by 2, rounding away from 0. if mpz_sgn(lower_mpz) >= 0: mpz_cdiv_q_2exp(lower_mpz, lower_mpz, 1) else: mpz_fdiv_q_2exp(lower_mpz, lower_mpz, 1) # and cur_error to hold the error: mpz_cdiv_q_2exp(cur_error, cur_error, 1) cdef char *tmp_cstr tmp_cstr = <char *>PyMem_Malloc(mpz_sizeinbase(lower_mpz, base) + 2) if tmp_cstr == NULL: raise MemoryError("Unable to allocate memory for the mantissa of an interval") mpz_get_str(tmp_cstr, base, lower_mpz) digits = strlen(tmp_cstr) if tmp_cstr[0] == '-': digits -= 1 mant_string = <object> PyString_FromString(tmp_cstr+1) sign_string = '-' else: mant_string = <object> PyString_FromString(tmp_cstr) sign_string = '' PyMem_Free(tmp_cstr) if error_digits == 0: error_string = '' else: tmp_cstr = <char *>PyMem_Malloc(mpz_sizeinbase(cur_error, 10) + 2) if tmp_cstr == NULL: raise MemoryError("Unable to allocate memory for the error of an interval") mpz_get_str(tmp_cstr, 10, cur_error) error_string = <object> PyString_FromString(tmp_cstr) PyMem_Free(tmp_cstr) mpz_clear(lower_mpz) mpz_clear(upper_mpz) mpz_clear(tmp) mpz_clear(cur_error) mpz_clear(max_error) cdef bint scientific = prefer_sci # If the exponent is >0, we must use scientific notation. For # instance, RIF(10, 30) gives 2.?e1; we couldn't write that # number in the question syntax (with no error digits) without # scientific notation. if expo > 0: scientific = True # If we use scientific notation, we put the radix point to the # right of the first digit; that would give us an exponent of: cdef mp_exp_t sci_expo = expo + digits - 1 if abs(sci_expo) >= 6: scientific = True if scientific: return '%s%s.%s?%s%s%s'%(sign_string, mant_string[0], mant_string[1:], error_string, e, sci_expo) if expo + digits <= 0: return '%s0.%s%s?%s'%(sign_string, '0' * -(expo + digits), mant_string, error_string) return '%s%s.%s?%s'%(sign_string, mant_string[:expo+digits], mant_string[expo+digits:], error_string) def __copy__(self): """ Return copy of ``self`` - since ``self`` is immutable, we just return ``self`` again. EXAMPLES:: sage: a = RIF(3.5) sage: copy(a) is a True """ return self # Interval-specific functions def lower(self, rnd=None): """ Return the lower bound of this interval INPUT: - ``rnd`` -- (string) the rounding mode - ``'RNDN'`` -- round to nearest - ``'RNDD'`` -- (default) round towards minus infinity - ``'RNDZ'`` -- round towards zero - ``'RNDU'`` -- round towards plus infinity The rounding mode does not affect the value returned as a floating-point number, but it does control which variety of ``RealField`` the returned number is in, which affects printing and subsequent operations. EXAMPLES:: sage: R = RealIntervalField(13) sage: R.pi().lower().str(truncate=False) '3.1411' :: sage: x = R(1.2,1.3); x.str(style='brackets') '[1.1999 .. 1.3001]' sage: x.lower() 1.19 sage: x.lower('RNDU') 1.20 sage: x.lower('RNDN') 1.20 sage: x.lower('RNDZ') 1.19 sage: x.lower().parent() Real Field with 13 bits of precision and rounding RNDD sage: x.lower('RNDU').parent() Real Field with 13 bits of precision and rounding RNDU sage: x.lower() == x.lower('RNDU') True """ cdef RealNumber x if rnd is None: x = (<RealIntervalField_class>self._parent).__lower_field._new() else: x = (<RealField_class>(self._parent._real_field(rnd)))._new() mpfi_get_left(x.value, self.value) return x def upper(self, rnd=None): """ Return the upper bound of ``self`` INPUT: - ``rnd`` -- (string) the rounding mode - ``'RNDN'`` -- round to nearest - ``'RNDD'`` -- (default) round towards minus infinity - ``'RNDZ'`` -- round towards zero - ``'RNDU'`` -- round towards plus infinity The rounding mode does not affect the value returned as a floating-point number, but it does control which variety of ``RealField`` the returned number is in, which affects printing and subsequent operations. EXAMPLES:: sage: R = RealIntervalField(13) sage: R.pi().upper().str(truncate=False) '3.1417' :: sage: R = RealIntervalField(13) sage: x = R(1.2,1.3); x.str(style='brackets') '[1.1999 .. 1.3001]' sage: x.upper() 1.31 sage: x.upper('RNDU') 1.31 sage: x.upper('RNDN') 1.30 sage: x.upper('RNDD') 1.30 sage: x.upper('RNDZ') 1.30 sage: x.upper().parent() Real Field with 13 bits of precision and rounding RNDU sage: x.upper('RNDD').parent() Real Field with 13 bits of precision and rounding RNDD sage: x.upper() == x.upper('RNDD') True """ cdef RealNumber x if rnd is None: x = (<RealIntervalField_class>self._parent).__upper_field._new() else: x = ((<RealField_class>self._parent._real_field(rnd)))._new() mpfi_get_right(x.value, self.value) return x def endpoints(self, rnd=None): """ Return the lower and upper endpoints of ``self``. EXAMPLES:: sage: RIF(1,2).endpoints() (1.00000000000000, 2.00000000000000) sage: RIF(pi).endpoints() (3.14159265358979, 3.14159265358980) sage: a = CIF(RIF(1,2), RIF(3,4)) sage: a.real().endpoints() (1.00000000000000, 2.00000000000000) As with ``lower()`` and ``upper()``, a rounding mode is accepted:: sage: RIF(1,2).endpoints('RNDD')[0].parent() Real Field with 53 bits of precision and rounding RNDD """ return self.lower(rnd), self.upper(rnd) def absolute_diameter(self): """ The diameter of this interval (for `[a .. b]`, this is `b-a`), rounded upward, as a :class:`RealNumber`. EXAMPLES:: sage: RIF(1, pi).absolute_diameter() 2.14159265358979 """ cdef RealNumber x x = (<RealIntervalField_class>self._parent).__middle_field._new() mpfi_diam_abs(x.value, self.value) return x def relative_diameter(self): """ The relative diameter of this interval (for `[a .. b]`, this is `(b-a)/((a+b)/2)`), rounded upward, as a :class:`RealNumber`. EXAMPLES:: sage: RIF(1, pi).relative_diameter() 1.03418797197910 """ cdef RealNumber x x = (<RealIntervalField_class>self._parent).__middle_field._new() mpfi_diam_rel(x.value, self.value) return x def diameter(self): """ If 0 is in ``self``, then return :meth:`absolute_diameter()`, otherwise return :meth:`relative_diameter()`. EXAMPLES:: sage: RIF(1, 2).diameter() 0.666666666666667 sage: RIF(1, 2).absolute_diameter() 1.00000000000000 sage: RIF(1, 2).relative_diameter() 0.666666666666667 sage: RIF(pi).diameter() 1.41357985842823e-16 sage: RIF(pi).absolute_diameter() 4.44089209850063e-16 sage: RIF(pi).relative_diameter() 1.41357985842823e-16 sage: (RIF(pi) - RIF(3, 22/7)).diameter() 0.142857142857144 sage: (RIF(pi) - RIF(3, 22/7)).absolute_diameter() 0.142857142857144 sage: (RIF(pi) - RIF(3, 22/7)).relative_diameter() 2.03604377705518 """ cdef RealNumber x x = (<RealIntervalField_class>self._parent).__middle_field._new() mpfi_diam(x.value, self.value) return x def fp_rank_diameter(self): r""" Computes the diameter of this interval in terms of the "floating-point rank". The floating-point rank is the number of floating-point numbers (of the current precision) contained in the given interval, minus one. An ``fp_rank_diameter`` of 0 means that the interval is exact; an ``fp_rank_diameter`` of 1 means that the interval is as tight as possible, unless the number you're trying to represent is actually exactly representable as a floating-point number. EXAMPLES:: sage: RIF(pi).fp_rank_diameter() 1 sage: RIF(12345).fp_rank_diameter() 0 sage: RIF(-sqrt(2)).fp_rank_diameter() 1 sage: RIF(5/8).fp_rank_diameter() 0 sage: RIF(5/7).fp_rank_diameter() 1 sage: a = RIF(pi)^12345; a 2.06622879260?e6137 sage: a.fp_rank_diameter() 30524 sage: (RIF(sqrt(2)) - RIF(sqrt(2))).fp_rank_diameter() 9671406088542672151117826 # 32-bit 41538374868278620559869609387229186 # 64-bit Just because we have the best possible interval, doesn't mean the interval is actually small:: sage: a = RIF(pi)^12345678901234567890; a [2.0985787164673874e323228496 .. +infinity] # 32-bit [5.8756537891115869e1388255822130839282 .. +infinity] # 64-bit sage: a.fp_rank_diameter() 1 """ return self.lower().fp_rank_delta(self.upper()) def is_exact(self): """ Return whether this real interval is exact (i.e. contains exactly one real value). EXAMPLES:: sage: RIF(3).is_exact() True sage: RIF(2*pi).is_exact() False """ return mpfr_equal_p(&self.value.left, &self.value.right) def magnitude(self): """ The largest absolute value of the elements of the interval. EXAMPLES:: sage: RIF(-2, 1).magnitude() 2.00000000000000 sage: RIF(-1, 2).magnitude() 2.00000000000000 """ cdef RealNumber x x = (<RealIntervalField_class>self._parent).__middle_field._new() mpfi_mag(x.value, self.value) return x def mignitude(self): """ The smallest absolute value of the elements of the interval. EXAMPLES:: sage: RIF(-2, 1).mignitude() 0.000000000000000 sage: RIF(-2, -1).mignitude() 1.00000000000000 sage: RIF(3, 4).mignitude() 3.00000000000000 """ cdef RealNumber x x = (<RealIntervalField_class>self._parent).__middle_field._new() mpfi_mig(x.value, self.value) return x def center(self): """ Compute the center of the interval `[a .. b]` which is `(a+b) / 2`. EXAMPLES:: sage: RIF(1, 2).center() 1.50000000000000 """ cdef RealNumber x x = (<RealIntervalField_class>self._parent).__middle_field._new() mpfi_mid(x.value, self.value) return x def bisection(self): """ Returns the bisection of ``self`` into two intervals of half the size whose union is ``self`` and intersection is :meth:`center()`. EXAMPLES:: sage: a, b = RIF(1,2).bisection() sage: a.lower(), a.upper() (1.00000000000000, 1.50000000000000) sage: b.lower(), b.upper() (1.50000000000000, 2.00000000000000) sage: I = RIF(e, pi) sage: a, b = I.bisection() sage: a.intersection(b) == I.center() True sage: a.union(b).endpoints() == I.endpoints() True """ cdef RealIntervalFieldElement left = self._new() cdef RealIntervalFieldElement right = self._new() mpfr_set(&left.value.left, &self.value.left, GMP_RNDN) mpfi_mid(&left.value.right, self.value) mpfi_interv_fr(right.value, &left.value.right, &self.value.right) return left, right def alea(self): """ Return a floating-point number picked at random from the interval. EXAMPLES:: sage: RIF(1, 2).alea() # random 1.34696133696137 """ cdef RealNumber x x = (<RealIntervalField_class>self._parent).__middle_field._new() mpfi_alea(x.value, self.value) return x # def integer_part(self): # """ # If in decimal this number is written n.defg, returns n. # OUTPUT: # -- a Sage Integer # EXAMPLE: # sage: a = 119.41212 # sage: a.integer_part() # 119 # """ # s = self.str(base=32, no_sci=True) # i = s.find(".") # return Integer(s[:i], base=32) ######################## # Basic Arithmetic ######################## cpdef ModuleElement _add_(self, ModuleElement other): """ Add two real intervals with the same parent. EXAMPLES:: sage: R = RealIntervalField() sage: R(-1.5) + R(2.5) # indirect doctest 1 sage: R('-1.3') + R('2.3') 1.000000000000000? sage: (R(1, 2) + R(3, 4)).str(style='brackets') '[4.0000000000000000 .. 6.0000000000000000]' """ cdef RealIntervalFieldElement x x = self._new() mpfi_add(x.value, self.value, (<RealIntervalFieldElement>other).value) return x def __invert__(self): """ Return the multiplicative "inverse" of this interval. (Technically, non-precise intervals don't have multiplicative inverses.) EXAMPLES:: sage: v = RIF(2); v 2 sage: ~v 0.50000000000000000? sage: v * ~v 1 sage: v = RIF(1.5, 2.5); v.str(style='brackets') '[1.5000000000000000 .. 2.5000000000000000]' sage: (~v).str(style='brackets') '[0.39999999999999996 .. 0.66666666666666675]' sage: (v * ~v).str(style='brackets') '[0.59999999999999986 .. 1.6666666666666670]' sage: ~RIF(-1, 1) [-infinity .. +infinity] """ cdef RealIntervalFieldElement x x = self._new() mpfi_inv(x.value, self.value) return x cpdef ModuleElement _sub_(self, ModuleElement right): """ Subtract two real intervals with the same parent. EXAMPLES:: sage: R = RealIntervalField() sage: R(-1.5) - R(2.5) # indirect doctest -4 sage: R('-1.3') - R('2.7') -4.000000000000000? sage: (R(1, 2) - R(3, 4)).str(style='brackets') '[-3.0000000000000000 .. -1.0000000000000000]' """ cdef RealIntervalFieldElement x x = self._new() mpfi_sub(x.value, self.value, (<RealIntervalFieldElement>right).value) return x cpdef RingElement _mul_(self, RingElement right): """ Multiply two real intervals with the same parent. EXAMPLES:: sage: R = RealIntervalField() sage: R(-1.5) * R(2.5) # indirect doctest -3.7500000000000000? sage: R('-1.3') * R('2.3') -2.990000000000000? sage: (R(1, 2) * R(3, 4)).str(style='brackets') '[3.0000000000000000 .. 8.0000000000000000]' If two elements have different precision, arithmetic operations are performed after coercing to the lower precision:: sage: R20 = RealIntervalField(20) sage: R100 = RealIntervalField(100) sage: a = R20('393.3902834028345') sage: b = R100('393.3902834028345') sage: a 393.391? sage: b 393.390283402834500000000000000? sage: a*b 154756.? sage: b*a 154756.? sage: parent(b*a) Real Interval Field with 20 bits of precision """ cdef RealIntervalFieldElement x x = self._new() mpfi_mul(x.value, self.value, (<RealIntervalFieldElement>right).value) return x cpdef RingElement _div_(self, RingElement right): """ Divide ``self`` by ``right``, where both are real intervals with the same parent. EXAMPLES:: sage: R = RealIntervalField() sage: R(1)/R(3) # indirect doctest 0.3333333333333334? sage: R(1)/R(0) # since R(0) has no sign, gives the whole reals [-infinity .. +infinity] sage: R(1)/R(-1, 1) [-infinity .. +infinity] :: sage: R(-1.5) / R(2.5) -0.6000000000000000? sage: (R(1, 2) / R(3, 4)).str(style='brackets') '[0.25000000000000000 .. 0.66666666666666675]' """ cdef RealIntervalFieldElement x x = self._new() mpfi_div((<RealIntervalFieldElement>x).value, self.value, (<RealIntervalFieldElement>right).value) return x cpdef ModuleElement _neg_(self): """ Return the additive "inverse" of this interval. (Technically, non-precise intervals don't have additive inverses.) EXAMPLES:: sage: v = RIF(2); v 2 sage: -v # indirect doctest -2 sage: v + -v 0 sage: v = RIF(1.5, 2.5); v.str(error_digits=3) '2.000?500' sage: (-v).str(style='brackets') '[-2.5000000000000000 .. -1.5000000000000000]' sage: (v + -v).str(style='brackets') '[-1.0000000000000000 .. 1.0000000000000000]' """ cdef RealIntervalFieldElement x x = self._new() mpfi_neg(x.value, self.value) return x def __abs__(self): """ Return the absolute value of ``self``. EXAMPLES:: sage: RIF(2).__abs__() 2 sage: RIF(2.1).__abs__() 2.1000000000000001? sage: RIF(-2.1).__abs__() 2.1000000000000001? """ return self.abs() cdef RealIntervalFieldElement abs(RealIntervalFieldElement self): """ Return the absolute value of ``self``. EXAMPLES:: sage: RIF(2).abs() 2 sage: RIF(2.1).abs() 2.1000000000000001? sage: RIF(-2.1).abs() 2.1000000000000001? """ cdef RealIntervalFieldElement x x = self._new() mpfi_abs(x.value, self.value) return x def square(self): """ Return the square of ``self``. .. NOTE:: Squaring an interval is different than multiplying it by itself, because the square can never be negative. EXAMPLES:: sage: RIF(1, 2).square().str(style='brackets') '[1.0000000000000000 .. 4.0000000000000000]' sage: RIF(-1, 1).square().str(style='brackets') '[0.00000000000000000 .. 1.0000000000000000]' sage: (RIF(-1, 1) * RIF(-1, 1)).str(style='brackets') '[-1.0000000000000000 .. 1.0000000000000000]' """ cdef RealIntervalFieldElement x x = self._new() mpfi_sqr(x.value, self.value) return x # Bit shifting def _lshift_(RealIntervalFieldElement self, n): """ Return ``self*(2^n)`` for an integer ``n``. EXAMPLES:: sage: RIF(1.0)._lshift_(32) 4294967296 sage: RIF(1.5)._lshift_(2) 6 """ cdef RealIntervalFieldElement x if n > sys.maxsize: raise OverflowError, "n (=%s) must be <= %s"%(n, sys.maxsize) x = self._new() mpfi_mul_2exp(x.value, self.value, n) return x def __lshift__(x, y): """ Returns `x * 2^y`, for `y` an integer. Much faster than an ordinary multiplication. EXAMPLES:: sage: RIF(1.0) << 32 4294967296 """ if isinstance(x, RealIntervalFieldElement) and isinstance(y, (int,long, Integer)): return x._lshift_(y) return sage.structure.element.bin_op(x, y, operator.lshift) def _rshift_(RealIntervalFieldElement self, n): """ Return ``self/(2^n)`` for an integer ``n``. EXAMPLES:: sage: RIF(1.0)._rshift_(32) 2.3283064365386963?e-10 sage: RIF(1.5)._rshift_(2) 0.37500000000000000? """ if n > sys.maxsize: raise OverflowError("n (=%s) must be <= %s" % (n, sys.maxsize)) cdef RealIntervalFieldElement x x = self._new() mpfi_div_2exp(x.value, self.value, n) return x def __rshift__(x, y): """ Returns `x / 2^y`, for `y` an integer. Much faster than an ordinary division. EXAMPLES:: sage: RIF(1024.0) >> 14 0.062500000000000000? """ if isinstance(x, RealIntervalFieldElement) and \ isinstance(y, (int,long,Integer)): return x._rshift_(y) return sage.structure.element.bin_op(x, y, operator.rshift) def multiplicative_order(self): r""" Return `n` such that ``self^n == 1``. Only `\pm 1` have finite multiplicative order. EXAMPLES:: sage: RIF(1).multiplicative_order() 1 sage: RIF(-1).multiplicative_order() 2 sage: RIF(3).multiplicative_order() +Infinity """ if self == 1: return 1 elif self == -1: return 2 return sage.rings.infinity.infinity # def sign(self): # return mpfr_sgn(self.value) def precision(self): """ Returns the precision of ``self``. EXAMPLES:: sage: RIF(2.1).precision() 53 sage: RealIntervalField(200)(2.1).precision() 200 """ return (<RealIntervalField_class>self._parent).__prec prec = precision ################### # Rounding etc ################### # Not implemented on intervals (for good reason!) # def round(self): # """ # Rounds self to the nearest real number. There are 4 # rounding modes. They are # EXAMPLES: # RNDN -- round to nearest: # sage: R = RealIntervalField(20,False,'RNDN') # sage: R(22.454) # 22.454 # sage: R(22.491) # 22.490 # RNDZ -- round towards zero: # sage: R = RealIntervalField(20,False,'RNDZ') # sage: R(22.454) # 22.453 # sage: R(22.491) # 22.490 # RNDU -- round towards plus infinity: # sage: R = RealIntervalField(20,False,'RNDU') # sage: R(22.454) # 22.454 # sage: R(22.491) # 22.491 # RNDU -- round towards minus infinity: # sage: R = RealIntervalField(20,False,'RNDD') # sage: R(22.454) # 22.453 # sage: R(22.491) # 22.490 # """ # cdef RealIntervalFieldElement x # x = self._new() # mpfr_round(x.value, self.value) # return x def floor(self): """ Return the floor of ``self``. EXAMPLES:: sage: R = RealIntervalField() sage: (2.99).floor() 2 sage: (2.00).floor() 2 sage: floor(RR(-5/2)) -3 sage: R = RealIntervalField(100) sage: a = R(9.5, 11.3); a.str(style='brackets') '[9.5000000000000000000000000000000 .. 11.300000000000000710542735760101]' sage: floor(a).str(style='brackets') '[9.0000000000000000000000000000000 .. 11.000000000000000000000000000000]' sage: a.floor() 10.? sage: ceil(a) 11.? sage: a.ceil().str(style='brackets') '[10.000000000000000000000000000000 .. 12.000000000000000000000000000000]' """ return self.parent()(self.lower().floor(), self.upper().floor()) def ceil(self): """ Return the ceiling of ``self``. OUTPUT: integer EXAMPLES:: sage: (2.99).ceil() 3 sage: (2.00).ceil() 2 sage: (2.01).ceil() 3 sage: R = RealIntervalField(30) sage: a = R(-9.5, -11.3); a.str(style='brackets') '[-11.300000012 .. -9.5000000000]' sage: a.floor().str(style='brackets') '[-12.000000000 .. -10.000000000]' sage: a.ceil() -10.? sage: ceil(a).str(style='brackets') '[-11.000000000 .. -9.0000000000]' """ return self.parent()(self.lower().ceil(), self.upper().ceil()) ceiling = ceil # def trunc(self): # """ # Truncates this number # EXAMPLES: # sage: (2.99).trunc() # 2.00000000000000 # sage: (-0.00).trunc() # -0.000000000000000 # sage: (0.00).trunc() # 0.000000000000000 # """ # cdef RealIntervalFieldElement x # x = self._new() # mpfr_trunc(x.value, self.value) # return x # def frac(self): # """ # frac returns a real number > -1 and < 1. it satisfies the # relation: # x = x.trunc() + x.frac() # EXAMPLES: # sage: (2.99).frac() # 0.990000000000000 # sage: (2.50).frac() # 0.500000000000000 # sage: (-2.79).frac() # -0.790000000000000 # """ # cdef RealIntervalFieldElement x # x = self._new() # mpfr_frac(x.value, self.value, (<RealIntervalField>self._parent).rnd) # return x ########################################### # Conversions ########################################### # def __float__(self): # return mpfr_get_d(self.value, (<RealIntervalField>self._parent).rnd) # def __int__(self): # """ # Returns integer truncation of this real number. # """ # s = self.str(32) # i = s.find('.') # return int(s[:i], 32) # def __long__(self): # """ # Returns long integer truncation of this real number. # """ # s = self.str(32) # i = s.find('.') # return long(s[:i], 32) # def __complex__(self): # return complex(float(self)) # def _complex_number_(self): # return sage.rings.complex_field.ComplexField(self.prec())(self) # def _pari_(self): # return sage.libs.pari.all.pari.new_with_bits_prec(str(self), (<RealIntervalField>self._parent).__prec) def unique_floor(self): """ Returns the unique floor of this interval, if it is well defined, otherwise raises a ``ValueError``. OUTPUT: - an integer. EXAMPLES:: sage: RIF(pi).unique_floor() 3 sage: RIF(100*pi).unique_floor() 314 sage: RIF(100, 200).unique_floor() Traceback (most recent call last): ... ValueError: interval does not have a unique floor """ a, b = self.lower().floor(), self.upper().floor() if a == b: return a else: raise ValueError("interval does not have a unique floor") def unique_ceil(self): """ Returns the unique ceiling of this interval, if it is well defined, otherwise raises a ``ValueError``. OUTPUT: - an integer. EXAMPLES:: sage: RIF(pi).unique_ceil() 4 sage: RIF(100*pi).unique_ceil() 315 sage: RIF(100, 200).unique_ceil() Traceback (most recent call last): ... ValueError: interval does not have a unique ceil """ a, b = self.lower().ceil(), self.upper().ceil() if a == b: return a else: raise ValueError("interval does not have a unique ceil") def unique_round(self): """ Returns the unique round (nearest integer) of this interval, if it is well defined, otherwise raises a ``ValueError``. OUTPUT: - an integer. EXAMPLES:: sage: RIF(pi).unique_round() 3 sage: RIF(1000*pi).unique_round() 3142 sage: RIF(100, 200).unique_round() Traceback (most recent call last): ... ValueError: interval does not have a unique round (nearest integer) sage: RIF(1.2, 1.7).unique_round() Traceback (most recent call last): ... ValueError: interval does not have a unique round (nearest integer) sage: RIF(0.7, 1.2).unique_round() 1 sage: RIF(-pi).unique_round() -3 sage: (RIF(4.5).unique_round(), RIF(-4.5).unique_round()) (5, -5) TESTS:: sage: RIF(-1/2, -1/3).unique_round() Traceback (most recent call last): ... ValueError: interval does not have a unique round (nearest integer) sage: RIF(-1/2, 1/3).unique_round() Traceback (most recent call last): ... ValueError: interval does not have a unique round (nearest integer) sage: RIF(-1/3, 1/3).unique_round() 0 sage: RIF(-1/2, 0).unique_round() Traceback (most recent call last): ... ValueError: interval does not have a unique round (nearest integer) sage: RIF(1/2).unique_round() 1 sage: RIF(-1/2).unique_round() -1 sage: RIF(0).unique_round() 0 """ a, b = self.lower().round(), self.upper().round() if a == b: return a else: raise ValueError("interval does not have a unique round (nearest integer)") def unique_integer(self): """ Return the unique integer in this interval, if there is exactly one, otherwise raises a ``ValueError``. EXAMPLES:: sage: RIF(pi).unique_integer() Traceback (most recent call last): ... ValueError: interval contains no integer sage: RIF(pi, pi+1).unique_integer() 4 sage: RIF(pi, pi+2).unique_integer() Traceback (most recent call last): ... ValueError: interval contains more than one integer sage: RIF(100).unique_integer() 100 """ a, b = self.lower().ceil(), self.upper().floor() if a == b: return a elif a < b: raise ValueError, "interval contains more than one integer" else: raise ValueError, "interval contains no integer" def simplest_rational(self, low_open=False, high_open=False): """ Return the simplest rational in this interval. Given rationals `a / b` and `c / d` (both in lowest terms), the former is simpler if `b<d` or if `b = d` and `|a| < |c|`. If optional parameters ``low_open`` or ``high_open`` are ``True``, then treat this as an open interval on that end. EXAMPLES:: sage: RealIntervalField(10)(pi).simplest_rational() 22/7 sage: RealIntervalField(20)(pi).simplest_rational() 355/113 sage: RIF(0.123, 0.567).simplest_rational() 1/2 sage: RIF(RR(1/3).nextabove(), RR(3/7)).simplest_rational() 2/5 sage: RIF(1234/567).simplest_rational() 1234/567 sage: RIF(-8765/432).simplest_rational() -8765/432 sage: RIF(-1.234, 0.003).simplest_rational() 0 sage: RIF(RR(1/3)).simplest_rational() 6004799503160661/18014398509481984 sage: RIF(RR(1/3)).simplest_rational(high_open=True) Traceback (most recent call last): ... ValueError: simplest_rational() on open, empty interval sage: RIF(1/3, 1/2).simplest_rational() 1/2 sage: RIF(1/3, 1/2).simplest_rational(high_open=True) 1/3 sage: phi = ((RealIntervalField(500)(5).sqrt() + 1)/2) sage: phi.simplest_rational() == fibonacci(362)/fibonacci(361) True """ if mpfr_equal_p(&self.value.left, &self.value.right): if low_open or high_open: raise ValueError, 'simplest_rational() on open, empty interval' return self.lower().exact_rational() if mpfi_has_zero(self.value): return Rational(0) if mpfi_is_neg(self.value): return -(self._neg_().simplest_rational(low_open=high_open, high_open=low_open)) low = self.lower() high = self.upper() # First, we try using approximate arithmetic of slightly higher # precision. cdef RealIntervalFieldElement highprec highprec = RealIntervalField(int(self.prec() * 1.2))(self) cdef Rational try1 = highprec._simplest_rational_helper() # Note that to compute "try1 >= low", Sage converts try1 to a # floating-point number rounding down, and "try1 <= high" # rounds up (since "low" and "high" are in downward-rounding # and upward-rounding fields, respectively). if try1 >= low and try1 <= high: ok = True if low_open and (try1 == low.exact_rational()): ok = False if high_open and (try1 == high.exact_rational()): ok = False if ok: return try1 # We could try again with higher precision; instead, we # go directly to using exact arithmetic. return _simplest_rational_exact(low.exact_rational(), high.exact_rational(), low_open, high_open) cdef Rational _simplest_rational_helper(self): """ Returns the simplest rational in an interval which is either equal to or slightly larger than ``self``. We assume that both endpoints of ``self`` are nonnegative. """ low = self.lower() cdef RealIntervalFieldElement new_elt if low <= 1: if low == 0: return Rational(0) if self.upper() >= 1: return Rational(1) new_elt = ~self return ~(new_elt._simplest_rational_helper()) fl = low.floor() new_elt = self - fl return fl + new_elt._simplest_rational_helper() ########################################### # Comparisons: ==, !=, <, <=, >, >= ########################################### def is_NaN(self): """ Check to see if ``self`` is Not-a-Number ``NaN``. EXAMPLES:: sage: a = RIF(0) / RIF(0.0,0.00); a [.. NaN ..] sage: a.is_NaN() True """ return mpfi_nan_p(self.value) def __richcmp__(left, right, int op): """ Rich comparison between ``left`` and ``right``. For more information, see :mod:`sage.rings.real_mpfi`. EXAMPLES:: sage: RIF(0) < RIF(2) True sage: RIF(0, 1) < RIF(2, 3) True Because these are possible ranges, they are only equal if they are exact and follow inequality if the intervals are disjoint:: sage: RIF(2) == RIF(2) True sage: RIF(0, 1) == RIF(0, 1) False sage: RIF(0, 2) < RIF(2, 3) False sage: RIF(0, 2) > RIF(2, 3) False """ return (<Element>left)._richcmp(right, op) cdef _richcmp_c_impl(left, Element right, int op): """ Implements comparisons between intervals. (See the file header comment for more information on interval comparison.) EXAMPLES:: sage: 0 < RIF(1, 3) True sage: 1 < RIF(1, 3) False sage: 2 < RIF(1, 3) False sage: 4 < RIF(1, 3) False sage: RIF(0, 1/2) < RIF(1, 3) True sage: RIF(0, 1) < RIF(1, 3) False sage: RIF(0, 2) < RIF(1, 3) False sage: RIF(1, 2) < RIF(1, 3) False sage: RIF(1, 3) < 4 True sage: RIF(1, 3) < 3 False sage: RIF(1, 3) < 2 False sage: RIF(1, 3) < 0 False sage: 0 <= RIF(1, 3) True sage: 1 <= RIF(1, 3) True sage: 2 <= RIF(1, 3) False sage: 4 <= RIF(1, 3) False sage: RIF(0, 1/2) <= RIF(1, 3) True sage: RIF(0, 1) <= RIF(1, 3) True sage: RIF(0, 2) <= RIF(1, 3) False sage: RIF(1, 2) <= RIF(1, 3) False sage: RIF(1, 3) <= 4 True sage: RIF(1, 3) <= 3 True sage: RIF(1, 3) <= 2 False sage: RIF(1, 3) <= 0 False sage: RIF(1, 3) > 0 True sage: RIF(1, 3) > 1 False sage: RIF(1, 3) > 2 False sage: RIF(1, 3) > 4 False sage: RIF(1, 3) > RIF(0, 1/2) True sage: RIF(1, 3) > RIF(0, 1) False sage: RIF(1, 3) > RIF(0, 2) False sage: RIF(1, 3) > RIF(1, 2) False sage: 4 > RIF(1, 3) True sage: 3 > RIF(1, 3) False sage: 2 > RIF(1, 3) False sage: 0 > RIF(1, 3) False sage: RIF(1, 3) >= 0 True sage: RIF(1, 3) >= 1 True sage: RIF(1, 3) >= 2 False sage: RIF(1, 3) >= 4 False sage: RIF(1, 3) >= RIF(0, 1/2) True sage: RIF(1, 3) >= RIF(0, 1) True sage: RIF(1, 3) >= RIF(0, 2) False sage: RIF(1, 3) >= RIF(1, 2) False sage: 4 >= RIF(1, 3) True sage: 3 >= RIF(1, 3) True sage: 2 >= RIF(1, 3) False sage: 0 >= RIF(1, 3) False sage: 0 == RIF(0) True sage: 0 == RIF(1) False sage: 1 == RIF(0) False sage: 0 == RIF(0, 1) False sage: 1 == RIF(0, 1) False sage: RIF(0, 1) == RIF(0, 1) False sage: RIF(1) == 0 False sage: RIF(1) == 1 True sage: RIF(0) == RIF(0) True sage: RIF(pi) == RIF(pi) False sage: RIF(0, 1) == RIF(1, 2) False sage: RIF(1, 2) == RIF(0, 1) False sage: 0 != RIF(0) False sage: 0 != RIF(1) True sage: 1 != RIF(0) True sage: 0 != RIF(0, 1) False sage: 1 != RIF(0, 1) False sage: RIF(0, 1) != RIF(0, 1) False sage: RIF(1) != 0 True sage: RIF(1) != 1 False sage: RIF(0) != RIF(0) False sage: RIF(pi) != RIF(pi) False sage: RIF(0, 1) != RIF(1, 2) False sage: RIF(1, 2) != RIF(0, 1) False """ cdef RealIntervalFieldElement lt, rt lt = left rt = right if op == 0: #< return mpfr_less_p(&lt.value.right, &rt.value.left) elif op == 2: #== # a == b iff a<=b and b <= a # (this gives a result with two comparisons, where the # obvious approach would use three) return mpfr_lessequal_p(&lt.value.right, &rt.value.left) \ and mpfr_lessequal_p(&rt.value.right, &lt.value.left) elif op == 4: #> return mpfr_less_p(&rt.value.right, &lt.value.left) elif op == 1: #<= return mpfr_lessequal_p(&lt.value.right, &rt.value.left) elif op == 3: #!= return mpfr_less_p(&lt.value.right, &rt.value.left) \ or mpfr_less_p(&rt.value.right, &lt.value.left) elif op == 5: #>= return mpfr_lessequal_p(&rt.value.right, &lt.value.left) def __nonzero__(self): """ Return ``True`` if ``self`` is not known to be exactly zero. EXAMPLES:: sage: RIF(0).__nonzero__() False sage: RIF(1).__nonzero__() True sage: RIF(1, 2).__nonzero__() True sage: RIF(0, 1).__nonzero__() True sage: RIF(-1, 1).__nonzero__() True """ return not (mpfr_zero_p(&self.value.left) and mpfr_zero_p(&self.value.right)) def __cmp__(left, right): """ Compare two intervals lexicographically. Return 0 if they are the same interval, -1 if the second is larger, or 1 if the first is larger. EXAMPLES:: sage: cmp(RIF(0), RIF(1)) -1 sage: cmp(RIF(0, 1), RIF(1)) -1 sage: cmp(RIF(0, 1), RIF(1, 2)) -1 sage: cmp(RIF(0, 0.99999), RIF(1, 2)) -1 sage: cmp(RIF(1, 2), RIF(0, 1)) 1 sage: cmp(RIF(1, 2), RIF(0)) 1 sage: cmp(RIF(0, 1), RIF(0, 2)) -1 sage: cmp(RIF(0, 1), RIF(0, 1)) 0 sage: cmp(RIF(0, 1), RIF(0, 1/2)) 1 """ return (<Element>left)._cmp(right) cdef int _cmp_c_impl(left, Element right) except -2: """ Implements the lexicographic total order on intervals. """ cdef RealIntervalFieldElement lt, rt lt = left rt = right cdef int i i = mpfr_cmp(&lt.value.left, &rt.value.left) if i < 0: return -1 elif i > 0: return 1 i = mpfr_cmp(&lt.value.right, &rt.value.right) if i < 0: return -1 elif i > 0: return 1 else: return 0 def __contains__(self, other): """ Test whether one interval (or real number) is totally contained in another. EXAMPLES:: sage: RIF(0, 2) in RIF(1, 3) False sage: RIF(0, 2) in RIF(0, 2) True sage: RIF(1, 2) in RIF(0, 3) True sage: 1.0 in RIF(0, 2) True sage: pi in RIF(3.1415, 3.1416) True sage: 22/7 in RIF(3.1415, 3.1416) False """ cdef RealIntervalFieldElement other_intv cdef RealNumber other_rn if PY_TYPE_CHECK(other, RealIntervalFieldElement): other_intv = other return mpfi_is_inside(other_intv.value, self.value) elif PY_TYPE_CHECK(other, RealNumber): other_rn = other return mpfi_is_inside_fr(other_rn.value, self.value) try: other_intv = self._parent(other) return mpfi_is_inside(other_intv.value, self.value) except TypeError as msg: return False def contains_zero(self): """ Return ``True`` if ``self`` is an interval containing zero. EXAMPLES:: sage: RIF(0).contains_zero() True sage: RIF(1, 2).contains_zero() False sage: RIF(-1, 1).contains_zero() True sage: RIF(-1, 0).contains_zero() True """ return mpfi_has_zero(self.value) def overlaps(self, RealIntervalFieldElement other): """ Return ``True`` if ``self`` and other are intervals with at least one value in common. For intervals ``a`` and ``b``, we have ``a.overlaps(b)`` iff ``not(a!=b)``. EXAMPLES:: sage: RIF(0, 1).overlaps(RIF(1, 2)) True sage: RIF(1, 2).overlaps(RIF(0, 1)) True sage: RIF(0, 1).overlaps(RIF(2, 3)) False sage: RIF(2, 3).overlaps(RIF(0, 1)) False sage: RIF(0, 3).overlaps(RIF(1, 2)) True sage: RIF(0, 2).overlaps(RIF(1, 3)) True """ return mpfr_greaterequal_p(&self.value.right, &other.value.left) \ and mpfr_greaterequal_p(&other.value.right, &self.value.left) def intersection(self, other): """ Return the intersection of two intervals. If the intervals do not overlap, raises a ``ValueError``. EXAMPLES:: sage: RIF(1, 2).intersection(RIF(1.5, 3)).str(style='brackets') '[1.5000000000000000 .. 2.0000000000000000]' sage: RIF(1, 2).intersection(RIF(4/3, 5/3)).str(style='brackets') '[1.3333333333333332 .. 1.6666666666666668]' sage: RIF(1, 2).intersection(RIF(3, 4)) Traceback (most recent call last): ... ValueError: intersection of non-overlapping intervals """ cdef RealIntervalFieldElement x x = self._new() cdef RealIntervalFieldElement other_intv if PY_TYPE_CHECK(other, RealIntervalFieldElement): other_intv = other else: # Let type errors from _coerce_ propagate... other_intv = self._parent(other) mpfi_intersect(x.value, self.value, other_intv.value) if mpfr_less_p(&x.value.right, &x.value.left): raise ValueError, "intersection of non-overlapping intervals" return x def union(self, other): """ Return the union of two intervals, or of an interval and a real number (more precisely, the convex hull). EXAMPLES:: sage: RIF(1, 2).union(RIF(pi, 22/7)).str(style='brackets') '[1.0000000000000000 .. 3.1428571428571433]' sage: RIF(1, 2).union(pi).str(style='brackets') '[1.0000000000000000 .. 3.1415926535897936]' sage: RIF(1).union(RIF(0, 2)).str(style='brackets') '[0.00000000000000000 .. 2.0000000000000000]' sage: RIF(1).union(RIF(-1)).str(style='brackets') '[-1.0000000000000000 .. 1.0000000000000000]' """ cdef RealIntervalFieldElement x x = self._new() cdef RealIntervalFieldElement other_intv cdef RealNumber other_rn if PY_TYPE_CHECK(other, RealIntervalFieldElement): other_intv = other mpfi_union(x.value, self.value, other_intv.value) elif PY_TYPE_CHECK(other, RealNumber): other_rn = other mpfi_set(x.value, self.value) mpfi_put_fr(x.value, other_rn.value) else: # Let type errors from _coerce_ propagate... other_intv = self._parent(other) mpfi_union(x.value, self.value, other_intv.value) return x def min(self, _other): """ Return an interval containing the minimum of ``self`` and ``other``. EXAMPLES:: sage: a=RIF(-1, 1).min(0).endpoints() sage: a[0] == -1.0 and a[1].abs() == 0.0 # in MPFI, the sign of 0.0 is not specified True sage: RIF(-1, 1).min(pi).endpoints() (-1.00000000000000, 1.00000000000000) sage: RIF(-1, 1).min(RIF(-100, 100)).endpoints() (-100.000000000000, 1.00000000000000) sage: RIF(-1, 1).min(RIF(-100, 0)).endpoints() (-100.000000000000, 0.000000000000000) The generic min does not always do the right thing:: sage: min(0, RIF(-1, 1)) 0 sage: min(RIF(-1, 1), RIF(-100, 100)).endpoints() (-1.00000000000000, 1.00000000000000) """ cdef RealIntervalFieldElement other if isinstance(_other, RealIntervalFieldElement): other = <RealIntervalFieldElement>_other else: other = self._parent(_other) if mpfr_cmp(&self.value.right, &other.value.left) <= 0: return self elif mpfr_cmp(&other.value.right, &self.value.left) <= 0: return other cdef RealIntervalFieldElement x = self._new() mpfr_min(&x.value.left, &self.value.left, &other.value.left, GMP_RNDD) mpfr_min(&x.value.right, &self.value.right, &other.value.right, GMP_RNDU) return x def max(self, _other): """ Return an interval containing the maximum of ``self`` and ``other``. EXAMPLES:: sage: RIF(-1, 1).max(0).endpoints() (0.000000000000000, 1.00000000000000) sage: RIF(-1, 1).max(RIF(2, 3)).endpoints() (2.00000000000000, 3.00000000000000) sage: RIF(-1, 1).max(RIF(-100, 100)).endpoints() (-1.00000000000000, 100.000000000000) The generic max does not always do the right thing:: sage: max(0, RIF(-1, 1)) 0 sage: max(RIF(-1, 1), RIF(-100, 100)).endpoints() (-1.00000000000000, 1.00000000000000) """ cdef RealIntervalFieldElement other if isinstance(_other, RealIntervalFieldElement): other = <RealIntervalFieldElement>_other else: other = self._parent(_other) if mpfr_cmp(&self.value.right, &other.value.left) <= 0: return other elif mpfr_cmp(&other.value.right, &self.value.left) <= 0: return self cdef RealIntervalFieldElement x = self._new() mpfr_max(&x.value.left, &self.value.left, &other.value.left, GMP_RNDD) mpfr_max(&x.value.right, &self.value.right, &other.value.right, GMP_RNDU) return x ############################ # Special Functions ############################ def sqrt(self): """ Return a square root of ``self``. Raises an error if ``self`` is nonpositive. If you use :meth:`square_root()` then an interval will always be returned (though it will be ``NaN`` if self is nonpositive). EXAMPLES:: sage: r = RIF(4.0) sage: r.sqrt() 2 sage: r.sqrt()^2 == r True :: sage: r = RIF(4344) sage: r.sqrt() 65.90902821313633? sage: r.sqrt()^2 == r False sage: r in r.sqrt()^2 True sage: r.sqrt()^2 - r 0.?e-11 sage: (r.sqrt()^2 - r).str(style='brackets') '[-9.0949470177292824e-13 .. 1.8189894035458565e-12]' :: sage: r = RIF(-2.0) sage: r.sqrt() Traceback (most recent call last): ... ValueError: self (=-2) is not >= 0 :: sage: r = RIF(-2, 2) sage: r.sqrt() Traceback (most recent call last): ... ValueError: self (=0.?e1) is not >= 0 """ if self.lower() < 0: raise ValueError, "self (=%s) is not >= 0"%self return self.square_root() def square_root(self): """ Return a square root of ``self``. An interval will always be returned (though it will be ``NaN`` if self is nonpositive). EXAMPLES:: sage: r = RIF(-2.0) sage: r.square_root() [.. NaN ..] sage: r.sqrt() Traceback (most recent call last): ... ValueError: self (=-2) is not >= 0 """ cdef RealIntervalFieldElement x x = self._new() sig_on() mpfi_sqrt(x.value, self.value) sig_off() return x # MPFI does not have cbrt. # def cube_root(self): # """ # Return the cubic root (defined over the real numbers) of self. # EXAMPLES: # sage: r = 125.0; r.cube_root() # 5.00000000000000 # sage: r = -119.0 # sage: r.cube_root()^3 - r # illustrates precision loss # -0.0000000000000142108547152020 # """ # cdef RealIntervalFieldElement x # x = self._new() # sig_on() # mpfr_cbrt(x.value, self.value, (<RealIntervalField>self._parent).rnd) # sig_off() # return x # MPFI does not have pow. # def __pow(self, RealIntervalFieldElement exponent): # cdef RealIntervalFieldElement x # x = self._new() # sig_on() # mpfr_pow(x.value, self.value, exponent.value, (<RealIntervalField>self._parent).rnd) # sig_off() # if mpfr_nan_p(x.value): # return self._complex_number_()**exponent._complex_number_() # return x # def __pow__(self, exponent, modulus): # """ # Compute self raised to the power of exponent, rounded in # the direction specified by the parent of self. # If the result is not a real number, self and the exponent are # both coerced to complex numbers (with sufficient precision), # then the exponentiation is computed in the complex numbers. # Thus this function can return either a real or complex number. # EXAMPLES: # sage: R = RealIntervalField(30) # sage: a = R('1.23456') # sage: a^20 # 67.646297 # sage: a^a # 1.2971114 # sage: b = R(-1) # sage: b^(1/2) # 1.0000000*I # 32-bit # -0.00000000000000000010842021 + 0.99999999*I # 64-bit # """ # cdef RealIntervalFieldElement x # if not PY_TYPE_CHECK(self, RealIntervalFieldElement): # return self.__pow__(float(exponent)) # if not PY_TYPE_CHECK(exponent, RealIntervalFieldElement): # x = self # exponent = x._parent(exponent) # return self.__pow(exponent) def __pow__(self, exponent, modulus): """ Raise ``self`` to ``exponent``. EXAMPLES:: sage: R = RealIntervalField(17) sage: x = R((-e,pi)) sage: x2 = x^2; x2.lower(), x2.upper() (0.0000, 9.870) sage: x3 = x^3; x3.lower(), x3.upper() (-26.83, 31.01) """ if exponent == 2: return self.square() if isinstance(exponent, (int, long, Integer)): q, r = divmod (exponent, 2) if r == 0: # x^(2q) = (x^q)^2 xq = sage.rings.ring_element.RingElement.__pow__(self, q) return xq.abs().square() else: return sage.rings.ring_element.RingElement.__pow__(self, exponent) return (self.log() * exponent).exp() def log(self, base='e'): """ Return the logarithm of ``self`` to the given ``base``. EXAMPLES:: sage: R = RealIntervalField() sage: r = R(2); r.log() 0.6931471805599453? sage: r = R(-2); r.log() 0.6931471805599453? + 3.141592653589794?*I """ cdef RealIntervalFieldElement x if self < 0: return self.parent().complex_field()(self).log(base) if base == 'e': x = self._new() sig_on() mpfi_log(x.value, self.value) sig_off() return x elif base == 10: return self.log10() elif base == 2: return self.log2() else: return self.log() / (self.parent()(base)).log() def log2(self): """ Return ``log`` to the base 2 of ``self``. EXAMPLES:: sage: r = RIF(16.0) sage: r.log2() 4 :: sage: r = RIF(31.9); r.log2() 4.995484518877507? :: sage: r = RIF(0.0, 2.0) sage: r.log2() [-infinity .. 1.0000000000000000] """ cdef RealIntervalFieldElement x if self < 0: return self.parent().complex_field()(self).log(2) x = self._new() sig_on() mpfi_log2(x.value, self.value) sig_off() return x def log10(self): """ Return log to the base 10 of ``self``. EXAMPLES:: sage: r = RIF(16.0); r.log10() 1.204119982655925? sage: r.log() / log(10.0) 1.204119982655925? :: sage: r = RIF(39.9); r.log10() 1.600972895686749? :: sage: r = RIF(0.0) sage: r.log10() [-infinity .. -infinity] :: sage: r = RIF(-1.0) sage: r.log10() 1.364376353841841?*I """ cdef RealIntervalFieldElement x if self < 0: return self.parent().complex_field()(self).log(10) x = self._new() sig_on() mpfi_log10(x.value, self.value) sig_off() return x def exp(self): r""" Returns `e^\mathtt{self}` EXAMPLES:: sage: r = RIF(0.0) sage: r.exp() 1 :: sage: r = RIF(32.3) sage: a = r.exp(); a 1.065888472748645?e14 sage: a.log() 32.30000000000000? :: sage: r = RIF(-32.3) sage: r.exp() 9.38184458849869?e-15 """ cdef RealIntervalFieldElement x x = self._new() sig_on() mpfi_exp(x.value, self.value) sig_off() return x def exp2(self): """ Returns `2^\mathtt{self}` EXAMPLES:: sage: r = RIF(0.0) sage: r.exp2() 1 :: sage: r = RIF(32.0) sage: r.exp2() 4294967296 :: sage: r = RIF(-32.3) sage: r.exp2() 1.891172482530207?e-10 """ cdef RealIntervalFieldElement x x = self._new() sig_on() mpfi_exp2(x.value, self.value) sig_off() return x def is_int(self): r""" Checks to see whether this interval includes exactly one integer. OUTPUT: If this contains exactly one integer, it returns the tuple ``(True, n)``, where ``n`` is that integer; otherwise, this returns ``(False, None)``. EXAMPLES:: sage: a = RIF(0.8,1.5) sage: a.is_int() (True, 1) sage: a = RIF(1.1,1.5) sage: a.is_int() (False, None) sage: a = RIF(1,2) sage: a.is_int() (False, None) sage: a = RIF(-1.1, -0.9) sage: a.is_int() (True, -1) sage: a = RIF(0.1, 1.9) sage: a.is_int() (True, 1) sage: RIF(+infinity,+infinity).is_int() (False, None) """ a = self.lower() if a.is_NaN() or a.is_infinity(): return False, None a = a.ceil() b = self.upper() if b.is_NaN() or b.is_infinity(): return False, None b = b.floor() if a == b: return True, a else: return False, None # MPFI does not have exp10. (Could easily be synthesized if an[...]