#!/usr/bin/python
# -*- coding: utf-8 -*-
"""Interval arithmetic in Python for a Mandelbrot set rendering example,
which is in intervalmand.py.

This interval arithmetic package is NOT SUITABLE FOR RELIABLE
COMPUTATION.  Reliable computation is a use of interval arithmetic to
avoid being fooled by errors introduced by the limited precision of
floating-point calculations, and implementing it requires control over
the floating-point processor's rounding mode.  I don't know how to
control floating-point rounding behavior in Python, so this module is
not useful for that.

Arithmetic on complex-number intervals with this package is about 100×
slower than with native Python complex numbers.

There are other implementations of interval arithmetic for Python,
which I haven't evaluated:

<https://code.google.com/p/pyinterval/> "An interval object consists
of a finite union of closed, possibly unbound, mathematical
intervals [on the extended real number set]."

<https://github.com/intiocean/pyinter> "a small and simple library
written in Python for performing interval and discontinous range
arithmetic."

<https://pypi.python.org/pypi/IntPy> "Maximum Accuracy Interval
Arithmetic, a mathematical tool for the solution of problems related
to numerical errors.  ...Only Real Intervals are available. No Complex
Intervals, no Interval Vectors and Matrixes, and no extensions of
basic functions. These will be our next work."  Dead since 2009.

"""

import math

class Interval(object):
    def __init__(self, min, max):
        self.min = min
        self.max = max

    def __iter__(self):
        yield self.min
        yield self.max

    def __add__(self, other):
        if not isinstance(other, Interval):
            other = interval(other)

        return Interval(self.min + other.min, self.max + other.max)

    def __iadd__(self, other):
        if not isinstance(other, Interval):
            other = interval(other)

        self.min += other.min
        self.max += other.max
        return self

    def __isub__(self, other):
        if not isinstance(other, Interval):
            other = interval(other)

        self.min -= other.max
        self.max -= other.min
        return self

    def __neg__(self):
        return Interval(-self.max, -self.min)

    def __sub__(self, other):
        return self + -interval(other)

    def __radd__(self, other):
        return self + other

    def __rsub__(self, other):
        return interval(other) - self

    def __mul__(self, other):
        if not isinstance(other, Interval):
            other = interval(other)

        a = self.min * other.min
        b = self.max * other.min
        c = self.min * other.max
        d = self.max * other.max
        return Interval(min(a, b, c, d), max(a, b, c, d))

    def __rmul__(self, other):
        return self * other

    def __div__(self, other):
        if 0 in other:
            return infinite_interval
        a = self.min / other.min
        b = self.max / other.min
        c = self.min / other.max
        d = self.max / other.max
        return Interval(min(a, b, c, d), max(a, b, c, d))

    def __rdiv__(self, other):
        return interval(other) / self

    def __contains__(self, other):
        return self.min <= other <= self.max

    def __abs__(self):
        if 0 in self:
            return Interval(0, max(abs(self.min), abs(self.max)))
        elif self.max < 0:
            return -self
        else:
            return self

    def log(self):
        return Interval(math.log(self.min) if self.min != 0 else -inf,
                        math.log(self.max) if self.max != 0 else -inf)

    def exp(self):
        try:
            max = math.exp(self.max)
        except OverflowError:
            print "overflow in exp", self.max
            max = inf
        return Interval(math.exp(self.min), max)

    def __pow__(self, other):
        "Handle some special cases of exponentiation that work on the real line."
        other = interval(other)
        if self.min >= 0:
            return (self.log() * other).exp()
        elif other.min != other.max or int(other.min) != other.min:
            # Some special cases can be handled on the real line, but
            # taking a negative number to a fractional power is not
            # among them.
            raise ValueError(self)

        try:
            a = self.min ** other.min
        except OverflowError:
            a = inf
        try:
            b = self.max ** other.min
        except OverflowError:
            b = inf

        if other.min % 2 != 0:
            # Odd integral powers increase monotonically, so this is easy:
            return Interval(a, b)
        elif 0 in self:
            # Even powers have separate positive, negative, and
            # zero-spanning cases; the positive case is handled above.
            # Zero-spanning includes zero:
            return Interval(0, max(a, b))
        else:
            # Finally, the negative case; decreasing, so backwards:
            return Interval(b, a)

    def span(self):
        return self.max - self.min

    def copy(self):
        return Interval(self.min, self.max)

inf = 1e1000
infinite_interval = Interval(-inf, inf)

def interval(object):
    "Coerce an object to an interval."
    if isinstance(object, Interval):
        return object
    else:
        return Interval(object, object)

class ComplexInterval(object):
    def __init__(self, re, im):
        self.re = re
        self.im = im

    def __iter__(self):
        yield self.re
        yield self.im

    def __add__(self, other):
        other = complex_interval(other)
        return ComplexInterval(self.re + other.re, self.im + other.im)

    def __neg__(self):
        return ComplexInterval(-self.re, -self.im)

    def __sub__(self, other):
        return self + -complex_interval(other)

    def __radd__(self, other):
        return self + other

    def __rsub__(self, other):
        return complex_interval(other) - self

    def __mul__(self, other):
        if not isinstance(other, ComplexInterval):
            other = complex_interval(other)

        return ComplexInterval(self.re * other.re - self.im * other.im,
                               self.im * other.re + self.re * other.im)

    def __imul__(self, other):
        if not isinstance(other, ComplexInterval):
            other = complex_interval(other)
        other_re = other.re
        other_im = other.im
        self_re = self.re
        self_im = self.im

        # self.re, self.im = (self_re * other_re - self_im * other_im,
        #                     self_im * other_re + self_re * other_im)

        imsneg = self_im * other_im
        reim = self_re * other_im
        self.re *= other_re
        self.re -= imsneg
        self.im *= other_re
        self.im += reim

        return self

    def __iadd__(self, other):
        if not isinstance(other, ComplexInterval):
            other = complex_interval(other)
        
        self.re += other.re
        self.im += other.im

        return self

    def __rmul__(self, other):
        return self * other

    # I forget how complex division works but I'm not using it right now,
    # so I'm leaving it unimplemented.
    def __rdiv__(self, other):
        return complex_interval(other) / self

    def __contains__(self, other):
        other = complex(other)
        return other.real in self.re and other.imag in self.im

    def __abs__(self):
        rv = self.abs_sq()**0.5
        assert type(rv) is Interval, rv
        return rv

    def abs_sq(self):
        return self.re*self.re + self.im*self.im

    def copy(self):
        return ComplexInterval(self.re.copy(), self.im.copy())

def complex_interval(object):
    "Coerce an object to a ComplexInterval."
    if isinstance(object, ComplexInterval):
        return object
    elif isinstance(object, Interval):
        return ComplexInterval(object, interval(0))
    elif isinstance(object, complex):
        return ComplexInterval(interval(object.real), interval(object.imag))
    elif isinstance(object, int) or isinstance(object, long) or isinstance(object, float):
        return ComplexInterval(interval(object), interval(0))
    else:
        raise TypeError(object)

def mandelbrot_escape_time(z0, maxiter=63):
    """Compute the escape time for the Mandelbrot set for an interval.

    This computes a conservative approximation to the number of
    iterations needed for points in the given complex interval to
    escape given the orbit of the Mandelbrot set.  It's guaranteed
    that the (real) interval it returns will be a superset of the
    number of iterations any number inside the (complex) interval
    given as argument z0 needs to escape the usual cutoff radius of 2
    from the origin.  However, it may be overbroad.

    """
    z = complex_interval(z0)
    c = z.copy()
    minimum = None
    for ii in range(maxiter):
        abs_sq_z = z.abs_sq()
        # print "absz, minimum", abs_z, minimum
        if minimum is None and abs_sq_z.max > 4:
            minimum = ii
        if abs_sq_z.min > 4:
            return Interval(minimum, ii)
        z *= z
        z += c

    # If we have reached maxiter, we might either be entirely inside
    # the (approximate) Mandelbrot set, or just partly.  If we're
    # entirely inside, minimum is still None.
    return Interval(minimum if minimum is not None else maxiter,
                    maxiter)

def main():
    "Simple test function."
    # print 0, interval(0)
    # print "log 0", interval(0).log()
    # print "2 log 0", 2 * interval(0).log()
    # print "exp(2 log 0)", (2 * interval(0).log()).exp()
    print "0²", interval(0)**2
    print abs(-Interval(3, 4) - Interval(0, 0.1))**2
    print abs(complex_interval(1j) * complex_interval(1j))
    print mandelbrot_escape_time(ComplexInterval(Interval(-2., 2.), Interval(-2., 2.)))
    print mandelbrot_escape_time(ComplexInterval(Interval(0., .25), Interval(0., .25)))

if __name__ == "__main__":
    import cgitb
    cgitb.enable(format='text')
    main()