Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
""" Real intervals with a fixed absolute precision """ from __future__ import print_function, absolute_import
from sage.ext.stdsage cimport PY_NEW
from sage.libs.gmp.mpz cimport *
from sage.structure.element cimport RingElement, ModuleElement, Element, FieldElement from sage.rings.ring cimport Field from sage.rings.integer cimport Integer
from sage.structure.parent cimport Parent from sage.structure.element cimport parent
cpdef inline Integer shift_floor(Integer x, long shift): r""" Return `x / 2^s` where `s` is the value of ``shift``, rounded towards `-\infty`. For internal use.
EXAMPLES::
sage: from sage.rings.real_interval_absolute import shift_floor sage: shift_floor(15, 2) 3 sage: shift_floor(-15, 2) -4 """
cpdef inline Integer shift_ceil(Integer x, long shift): r""" Return `x / 2^s` where `s` is the value of ``shift``, rounded towards `+\infty`. For internal use.
EXAMPLES::
sage: from sage.rings.real_interval_absolute import shift_ceil sage: shift_ceil(15, 2) 4 sage: shift_ceil(-15, 2) -3 sage: shift_ceil(32, 2) 8 sage: shift_ceil(-32, 2) -8 """
""" The only piece of data is the precision.
TESTS::
sage: from sage.rings.real_interval_absolute import RealIntervalAbsoluteField sage: RealIntervalAbsoluteField.create_key(1000) 1000 """
""" Ensures uniqueness.
TESTS::
sage: from sage.rings.real_interval_absolute import RealIntervalAbsoluteField sage: RealIntervalAbsoluteField(23) is RealIntervalAbsoluteField(23) # indirect doctest True """
cdef class RealIntervalAbsoluteField_class(Field): """ This field is similar to the :class:`RealIntervalField` except instead of truncating everything to a fixed relative precision, it maintains a fixed absolute precision.
Note that unlike the standard real interval field, elements in this field can have different size and experience coefficient blowup. On the other hand, it avoids precision loss on addition and subtraction. This is useful for, e.g., series computations for special functions.
EXAMPLES::
sage: from sage.rings.real_interval_absolute import RealIntervalAbsoluteField sage: R = RealIntervalAbsoluteField(10); R Real Interval Field with absolute precision 2^-10 sage: R(3/10) 0.300? sage: R(1000003/10) 100000.300? sage: R(1e100) + R(1) - R(1e100) 1 """
cdef long _absprec
def __init__(self, absprec): """ Initialize ``self``.
EXAMPLES::
sage: from sage.rings.real_interval_absolute import RealIntervalAbsoluteField sage: RealIntervalAbsoluteField(100) Real Interval Field with absolute precision 2^-100 sage: RealIntervalAbsoluteField(-100) Traceback (most recent call last): File "<ipython console>", line 1, in <module> File "real_interval_absolute.pyx", line 81, in sage.rings.real_interval_absolute.RealIntervalAbsoluteField.__init__ (sage/rings/real_interval_absolute.c:2463) ValueError: Absolute precision must be positive. """
def __reduce__(self): """ Used for pickling.
TESTS::
sage: from sage.rings.real_interval_absolute import RealIntervalAbsoluteField sage: R = RealIntervalAbsoluteField(100) sage: loads(dumps(R)) Real Interval Field with absolute precision 2^-100 """
def _element_constructor_(self, x): """ Construct an element with ``self`` as the parent.
TESTS::
sage: from sage.rings.real_interval_absolute import RealIntervalAbsoluteField sage: R = RealIntervalAbsoluteField(100) sage: R(1/2) # indirect doctest 0.50000000000000000000000000000000? """
cpdef _coerce_map_from_(self, R): """ Anything that coerces into the reals coerces into this ring.
TESTS::
sage: from sage.rings.real_interval_absolute import RealIntervalAbsoluteField sage: R = RealIntervalAbsoluteField(100) sage: R.has_coerce_map_from(RR) # indirect doctest True sage: R.has_coerce_map_from(QQ) True sage: R.has_coerce_map_from(Qp(5)) False
sage: R(1/2) + 100 100.5000000000000000000000000000000? sage: R(1/2) + 1.75 2.2500000000000000000000000000000?
sage: R10 = RealIntervalAbsoluteField(10) sage: R10(1/4) + R(1/4) 0.50000? """ else:
def _repr_(self): """ Return the string representation of ``self``.
TESTS::
sage: from sage.rings.real_interval_absolute import RealIntervalAbsoluteField sage: R = RealIntervalAbsoluteField(100) sage: print(R) Real Interval Field with absolute precision 2^-100 sage: R._repr_() 'Real Interval Field with absolute precision 2^-100' """
def absprec(self): """ Returns the absolute precision of self.
EXAMPLES::
sage: from sage.rings.real_interval_absolute import RealIntervalAbsoluteField sage: R = RealIntervalAbsoluteField(100) sage: R.absprec() 100 sage: RealIntervalAbsoluteField(5).absprec() 5 """
cdef inline shift_left(value, shift): """ Utility function for operands that don't support the ``<<`` operator. """ # Better than the OverflowError we would get from trying to multiply. else:
cdef class RealIntervalAbsoluteElement(FieldElement):
# This could be optimized by letting these be raw mpz_t. cdef Integer _mantissa # left endpoint cdef Integer _diameter
def __init__(self, RealIntervalAbsoluteField_class parent, value): """ Create a :class:`RealIntervalAbsoluteElement`.
EXAMPLES::
sage: from sage.rings.real_interval_absolute import RealIntervalAbsoluteField sage: R = RealIntervalAbsoluteField(50) sage: R(1) 1 sage: R(1/3) 0.333333333333334? sage: R(1.3) 1.300000000000000? sage: R(pi) 3.141592653589794? sage: R((11, 12)) 12.? sage: R((11, 11.00001)) 11.00001?
sage: R100 = RealIntervalAbsoluteField(100) sage: R(R100((5,6))) 6.? sage: R100(R((5,6))) 6.? sage: RIF(CIF(NaN)) [.. NaN ..] """
else:
else: else: except OverflowError: raise TypeError(type(value))
def __reduce__(self): """ Used for pickling.
EXAMPLES::
sage: from sage.rings.real_interval_absolute import RealIntervalAbsoluteField sage: R = RealIntervalAbsoluteField(50) sage: loads(dumps(R(1/16))) 0.06250000000000000? sage: R = RealIntervalAbsoluteField(100) sage: loads(dumps(R(1/3))) 0.333333333333333333333333333334? sage: loads(dumps(R(pi))).endpoints() == R(pi).endpoints() True """
cdef _new_c(self, Integer _mantissa, Integer _diameter): cdef RealIntervalAbsoluteElement x
cpdef lower(self): """ Return the lower bound of ``self``.
EXAMPLES::
sage: from sage.rings.real_interval_absolute import RealIntervalAbsoluteField sage: R = RealIntervalAbsoluteField(50) sage: R(1/4).lower() 1/4 """
cpdef midpoint(self): """ Return the midpoint of ``self``.
EXAMPLES::
sage: from sage.rings.real_interval_absolute import RealIntervalAbsoluteField sage: R = RealIntervalAbsoluteField(100) sage: R(1/4).midpoint() 1/4 sage: R(pi).midpoint() 7964883625991394727376702227905/2535301200456458802993406410752 sage: R(pi).midpoint().n() 3.14159265358979 """
cpdef upper(self): """ Return the upper bound of ``self``.
EXAMPLES::
sage: from sage.rings.real_interval_absolute import RealIntervalAbsoluteField sage: R = RealIntervalAbsoluteField(50) sage: R(1/4).upper() 1/4 """
cpdef absolute_diameter(self): """ Return the diameter ``self``.
EXAMPLES::
sage: from sage.rings.real_interval_absolute import RealIntervalAbsoluteField sage: R = RealIntervalAbsoluteField(10) sage: R(1/4).absolute_diameter() 0 sage: a = R(pi) sage: a.absolute_diameter() 1/1024 sage: a.upper() - a.lower() 1/1024 """
cpdef endpoints(self): """ Return the left and right endpoints of ``self``, as a tuple.
EXAMPLES::
sage: from sage.rings.real_interval_absolute import RealIntervalAbsoluteField sage: R = RealIntervalAbsoluteField(10) sage: R(1/4).endpoints() (1/4, 1/4) sage: R((1,2)).endpoints() (1, 2) """
def _real_mpfi_(self, R): """ Create a (relative) real interval out of this absolute real interval.
EXAMPLES::
sage: from sage.rings.real_interval_absolute import RealIntervalAbsoluteField sage: R = RealIntervalAbsoluteField(10) sage: R(1/2)._real_mpfi_(RIF) 0.50000000000000000?
sage: a = RealIntervalAbsoluteField(100)(1/3) sage: RIF(a) 0.3333333333333334? """
cpdef long mpfi_prec(self): """ Return the precision needed to represent this value as an mpfi interval.
EXAMPLES::
sage: from sage.rings.real_interval_absolute import RealIntervalAbsoluteField sage: R = RealIntervalAbsoluteField(10) sage: R(10).mpfi_prec() 14 sage: R(1000).mpfi_prec() 20 """
def _repr_(self): """ Leverage real interval printing.
TESTS::
sage: from sage.rings.real_interval_absolute import RealIntervalAbsoluteField sage: R = RealIntervalAbsoluteField(10) sage: R(1/3) # indirect doctest 0.334? sage: R(10^50/3) 3.3333333333333333333333333333333333333333333333333334?e49 sage: R(0) 0 """
def __hash__(self): """ Hash to the midpoint of the interval.
TESTS::
sage: from sage.rings.real_interval_absolute import RealIntervalAbsoluteField sage: R = RealIntervalAbsoluteField(10) sage: hash(R(10)) 10 sage: hash(R((11,13))) 12 sage: hash(R(1/4)) == hash(1/4) True sage: hash(R(pi)) 891658780 # 32-bit 532995478001132060 # 64-bit """
def __contains__(self, x): """ Return whether the given value lies in this interval.
EXAMPLES::
sage: from sage.rings.real_interval_absolute import RealIntervalAbsoluteField sage: R = RealIntervalAbsoluteField(50) sage: 1 in R((1,2)) True sage: 2 in R((1,2)) True sage: 3 in R((1,2)) False sage: 1.75 in R((1,2)) True """
cpdef bint is_positive(self): """ Return whether ``self`` is definitely positive.
EXAMPLES::
sage: from sage.rings.real_interval_absolute import RealIntervalAbsoluteField sage: R = RealIntervalAbsoluteField(10) sage: R(10).is_positive() True sage: R((10,11)).is_positive() True sage: R((0,11)).is_positive() False sage: R((-10,11)).is_positive() False sage: R((-10,-1)).is_positive() False sage: R(pi).is_positive() True """
cpdef bint contains_zero(self): """ Return whether ``self`` contains zero.
EXAMPLES::
sage: from sage.rings.real_interval_absolute import RealIntervalAbsoluteField sage: R = RealIntervalAbsoluteField(10) sage: R(10).contains_zero() False sage: R((10,11)).contains_zero() False sage: R((0,11)).contains_zero() True sage: R((-10,11)).contains_zero() True sage: R((-10,-1)).contains_zero() False sage: R((-10,0)).contains_zero() True sage: R(pi).contains_zero() False """
cpdef bint is_negative(self): """ Return whether ``self`` is definitely negative.
EXAMPLES::
sage: from sage.rings.real_interval_absolute import RealIntervalAbsoluteField sage: R = RealIntervalAbsoluteField(100) sage: R(10).is_negative() False sage: R((10,11)).is_negative() False sage: R((0,11)).is_negative() False sage: R((-10,11)).is_negative() False sage: R((-10,-1)).is_negative() True sage: R(pi).is_negative() False """
cdef bint is_exact(self): return not self._diameter
def __nonzero__(self): """ Return ``True`` for anything except exact zero.
EXAMPLES::
sage: from sage.rings.real_interval_absolute import RealIntervalAbsoluteField sage: R = RealIntervalAbsoluteField(10) sage: bool(R(1)) True sage: bool(R(0)) False sage: bool(R((0,1))) True """
def __neg__(self): """ TESTS::
sage: from sage.rings.real_interval_absolute import RealIntervalAbsoluteField sage: R = RealIntervalAbsoluteField(100) sage: -R(1/2) -0.50000000000000000000000000000000? sage: -R((101,102)) -102.? """
def __abs__(self): """ EXAMPLES::
sage: from sage.rings.real_interval_absolute import RealIntervalAbsoluteField sage: R = RealIntervalAbsoluteField(100) sage: abs(-R(1/4)) 0.2500000000000000000000000000000? """
cpdef abs(self): """ Return the absolute value of ``self``.
EXAMPLES::
sage: from sage.rings.real_interval_absolute import RealIntervalAbsoluteField sage: R = RealIntervalAbsoluteField(100) sage: R(1/3).abs() 0.333333333333333333333333333334? sage: R(-1/3).abs() 0.333333333333333333333333333334? sage: R((-1/3, 1/2)).abs() 1.? sage: R((-1/3, 1/2)).abs().endpoints() (0, 1/2) sage: R((-3/2, 1/2)).abs().endpoints() (0, 3/2) """ else:
cpdef _add_(self, _other): """ TESTS::
sage: from sage.rings.real_interval_absolute import RealIntervalAbsoluteField sage: R = RealIntervalAbsoluteField(10) sage: R(1) + R(2) # indirect doctest 3 sage: R(1e100) + R(0.1) + R(-1e100) 0.100? sage: (R((1,2)) + 1).endpoints() (2, 3) sage: (R((1,2)) + R((-10,0))).endpoints() (-9, 2) """
cpdef _sub_(self, _other): """ TESTS::
sage: from sage.rings.real_interval_absolute import RealIntervalAbsoluteField sage: R = RealIntervalAbsoluteField(10) sage: R(1) - R(2) # indirect doctest -1 sage: R(1e100) - R(0.1) - R(1e100) -0.100? sage: (R((1,2)) - 1).endpoints() (0, 1) sage: (R((1,2)) - R((10,100))).endpoints() (-99, -8) sage: R(pi) - R(pi) 0.000? """
cpdef _mul_(self, _other): """ TESTS::
sage: from sage.rings.real_interval_absolute import RealIntervalAbsoluteField sage: R = RealIntervalAbsoluteField(10) sage: R(2) * R(3) # indirect doctest 6 sage: R(2) * R(-3) -6 sage: elts = [R((left, right)) for left in [-2..2] for right in [left+1..2]] sage: elts [-2.?, -1.?, 0.?e1, 0.?e1, -1.?, 0.?, 0.?e1, 1.?, 1.?, 2.?] sage: for a in elts: ....: for b in elts: ....: if (a*b).lower() != (a._real_mpfi_(RIF)*b._real_mpfi_(RIF)).lower(): ....: print(a, b) ....: if (a*b).upper() != (a._real_mpfi_(RIF)*b._real_mpfi_(RIF)).upper(): ....: print(a, b) sage: R(pi) * R(pi) - R(pi^2) 0.00? """
# Break symmetry.
else: res = self._new_c(shift_floor((self._mantissa + self._diameter) * other._mantissa, absprec), shift_ceil((self._mantissa + self._diameter) * other._diameter, absprec)) else: # They both contain zero.
cpdef _acted_upon_(self, x, bint self_on_left): """ ``Absprec * relprec -> absprec`` works better than coercing both operands to absolute precision first.
EXAMPLES::
sage: from sage.rings.real_interval_absolute import RealIntervalAbsoluteField sage: R = RealIntervalAbsoluteField(4) sage: 1/4 * R(100) 25 sage: 1/100 * R(100) 1 sage: R(1/100) * R(100) 1.?e1 sage: RIF(1/100) * R(100) 1.0?
sage: R(1.5)._acted_upon_(3, True) 4.500? """ neg = self._acted_upon_(-x, self_on_left) return None if neg is None else -neg return self * RealIntervalAbsoluteElement(self._parent, (x.lower(), x.upper())) # Remember, x > 0 else: left = (self._mantissa * x.upper()).floor() right = ((self._mantissa + self._diameter) * x.lower()).ceil()
def __invert__(self): """ EXAMPLES::
sage: from sage.rings.real_interval_absolute import RealIntervalAbsoluteField sage: R = RealIntervalAbsoluteField(10) sage: ~R(2) 0.50000? sage: ~~R(2) 2 sage: ~R(3) 0.334? sage: ~~R(3) 3.00?
sage: R = RealIntervalAbsoluteField(200) sage: ~R(1e10) 1.00000000000000000000000000000000000000000000000000?e-10 sage: ~~R(1e10) 1.00000000000000000000000000000000000000000000000000?e10 sage: (~R((1,2))).endpoints() (1/2, 1) sage: (~R((1/4,8))).endpoints() (1/8, 4) sage: R(1/pi) - 1/R(pi) 0.?e-60 """ raise ZeroDivisionError("Inversion of an interval containing zero.") cdef bint negate self = -self negate = True else:
# Let our (positive) interval be [2^-B a, 2^-B b]. # Then its inverse is [2^-B 2^(2B)/b , 2^-B 2^(2B)/a].
cdef mpz_t scaling_factor # Use diameter as temp value for right endpoint... # Divide a second time to get the rounding correct... finally:
res = -res
cpdef _div_(self, _other): """ TESTS::
sage: from sage.rings.real_interval_absolute import RealIntervalAbsoluteField sage: R = RealIntervalAbsoluteField(10) sage: R(5)/R(2) # indirect doctest 2.5000? sage: a = R((5,6))/R((2,4)) sage: a.endpoints() (5/4, 3) sage: a = R((-5,6))/R((2,4)) sage: a.endpoints() (-5/2, 3) sage: R(1e100) / R(2e100) 0.50000? """ raise ZeroDivisionError("Division by an interval containing zero.")
cdef mpz_t temp
negate = True self = -self
else:
finally:
res = -res
def __lshift__(RealIntervalAbsoluteElement self, long n): """ EXAMPLES::
sage: from sage.rings.real_interval_absolute import RealIntervalAbsoluteField sage: R = RealIntervalAbsoluteField(10) sage: R(1) << 2 4 sage: R(3) << -2 0.75000? sage: (R((1/2, 5)) << 10).endpoints() (512, 5120) """
def __rshift__(RealIntervalAbsoluteElement self, long n): """ EXAMPLES::
sage: from sage.rings.real_interval_absolute import RealIntervalAbsoluteField sage: R = RealIntervalAbsoluteField(10) sage: R(1) >> 2 0.2500? sage: R(3) >> -2 12 sage: (R((1/2, 5)) >> 10).endpoints() (0, 5/1024) """
cdef shift(self, long n): else:
def __pow__(self, exponent, dummy): """ EXAMPLES::
sage: from sage.rings.real_interval_absolute import RealIntervalAbsoluteField sage: R = RealIntervalAbsoluteField(10) sage: R(10)^10 10000000000 sage: R(10)^100 10000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000 sage: R(10)^-2 0.010? sage: R(10)^-10 0.001? sage: R(100)^(1/2) 10.00? sage: (R((2,3))^2).endpoints() (4, 9) sage: (R((-5,3))^3).endpoints() (-125, 75) sage: (R((-5,3))^4).endpoints() (0, 625) """ cdef RealIntervalAbsoluteElement base except TypeError: base = exponent.parent()(self) if not exponent: raise ZeroDivisionError else: return base cdef double height exponent = exponent._real_mpfi_(RIF)
def __getattr__(self, name): """ EXAMPLES::
sage: from sage.rings.real_interval_absolute import RealIntervalAbsoluteField sage: R = RealIntervalAbsoluteField(100) sage: R(1).sin() 0.841470984807896506652502321631? sage: R(2).log() 0.693147180559945309417232121458? sage: R(1).exp().log() 1.00000000000000000000000000000?
sage: R((0,10)).sin().endpoints() (-1, 1) sage: R(1).sin().parent() Real Interval Field with absolute precision 2^-100 """ else:
def sqrt(self): """ Return the square root of ``self``.
EXAMPLES::
sage: from sage.rings.real_interval_absolute import RealIntervalAbsoluteField sage: R = RealIntervalAbsoluteField(100) sage: R(2).sqrt() 1.414213562373095048801688724210? sage: R((4,9)).sqrt().endpoints() (2, 3) """
cdef class MpfrOp: """ This class is used to endow absolute real interval field elements with all the methods of (relative) real interval field elements.
EXAMPLES::
sage: from sage.rings.real_interval_absolute import RealIntervalAbsoluteField sage: R = RealIntervalAbsoluteField(100) sage: R(1).sin() 0.841470984807896506652502321631? """ cdef object name cdef RealIntervalAbsoluteElement value
def __init__(self, value, name): """ EXAMPLES::
sage: from sage.rings.real_interval_absolute import RealIntervalAbsoluteField, MpfrOp sage: R = RealIntervalAbsoluteField(100) sage: MpfrOp(R(1), 'tan')() 1.557407724654902230506974807459? """
def __call__(self, *args): """ EXAMPLES::
sage: from sage.rings.real_interval_absolute import RealIntervalAbsoluteField, MpfrOp sage: R = RealIntervalAbsoluteField(100) sage: curried_log = MpfrOp(R(2), 'log') sage: curried_log() 0.693147180559945309417232121458? sage: curried_log(2) 1 """ absprec = min(absprec, (<RealIntervalAbsoluteField_class>(<RealIntervalAbsoluteElement>a)._parent)._absprec) relprec = max(absprec, (<RealIntervalAbsoluteElement>a).mpfi_prec()) parent = RealIntervalAbsoluteField(absprec) relprec = 53 new_args.append(a._real_mpfi_(R)) else: |