CoCalc Public Filestmp / 2015-03-19-092125.sagews
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)
0
sage: R = RealIntervalField(sci_not=1, prec=200)
sage: b = R('393.39203845902384098234098230948209384028340')
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.lower(), c.upper()) == (b.lower(), b.upper())
True
sage: b = R(-1)/R(0); b # same as above
[-infinity .. +infinity]
sage: (c.lower(), c.upper()) == (b.lower(), b.upper())
True
sage: b = R('[2 .. 3]'); b.str(error_digits=1)
'2.5?5e0'
0
sage: R = RealIntervalField(4000)
sage: s = 1/R(3)
sage: (t.upper(), t.lower()) == (s.upper(), s.lower())
True
[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
# 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:
# 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
########################

"""
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()
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

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[...]