"""bigfloat: floating-point numbers with very large range"""

# David Chiang
# 2011 July 18

import math

_LOG_2 = math.log(2.)
_LOG10_2 = math.log10(2.)

class bigfloat(object):
    __slots__ = ('sig','exp')
    def __init__(self, sig=0., exp=0):
        s, e = math.frexp(sig)
        self.sig = s
        self.exp = e+exp

    def _align(self, other):
        ss, se = self.sig, self.exp
        if isinstance(other, bigfloat):
            os, oe = other.sig, other.exp
        else:
            os, oe = other, 0
        if ss == 0.:
            se = oe
        elif os == 0.:
            pass
        elif se > oe:
            os = math.ldexp(os, oe-se)
        elif se < oe:
            ss = math.ldexp(ss, se-oe)
            se = oe
        return ss, os, se

    def __add__(self, other):
        ss, os, e = self._align(other)
        return bigfloat(ss+os, e)
    def __radd__(self, other):
        ss, os, e = self._align(other)
        return bigfloat(os+ss, e)

    def __sub__(self, other):
        ss, os, e = self._align(other)
        return bigfloat(ss-os, e)
    def __rsub__(self, other):
        ss, os, e = self._align(other)
        return bigfloat(os-ss, e)

    def __cmp__(self, other):
        ss, os, e = self._align(other)
        return cmp(ss, os)

    def __pos__(self):
        return bigfloat(self.sig, self.exp)
    def __neg__(self):
        return bigfloat(-self.sig, self.exp)
    def __abs__(self):
        return bigfloat(abs(self.sig), self.exp)

    def __mul__(self, other):
        if isinstance(other, bigfloat):
            return bigfloat(self.sig*other.sig, self.exp+other.exp)
        else:
            return bigfloat(self.sig*other, self.exp)
    def __rmul__(self, other):
        return bigfloat(self.sig*other, self.exp)

    def __div__(self, other):
        if isinstance(other, bigfloat):
            return bigfloat(self.sig/other.sig, self.exp-other.exp)
        else:
            return bigfloat(self.sig/other, self.exp)
    def __rdiv__(self, other):
        return bigfloat(other/self.sig, -self.exp)

    def __pow__(self, other):
        if isinstance(other, int) and abs(other) <= 10:
            return bigfloat(pow(self.sig, other), self.exp*other)
        else:
            return exp2(other*log2(self))
    def __rpow__(self, other):
        return exp2(float(self)*math.log(other)/_LOG_2)

    def __str__(self):
        if abs(self.exp) <= 10:
            return str(float(self))
        else:
            l = log10(abs(self))
            e10 = int(math.floor(l))
            s10 = math.pow(10.,l-e10)
            if self.sig < 0.: s10 = -s10
            return "%fe%+d" % (s10,e10)
    def __repr__(self):
        return "bigfloat(%s,%s)" % (self.sig, self.exp)
    def __float__(self):
        return math.ldexp(self.sig, self.exp)

def log2(b):
    return math.log(b.sig)/_LOG_2 + b.exp
def log(b):
    return math.log(b.sig) + b.exp * _LOG_2
def log10(b):
    return math.log10(b.sig) + b.exp * _LOG10_2
def exp2(f):
    e = int(math.floor(f))
    return bigfloat(math.pow(2.,f-e), e)
def exp(f):
    return exp2(f/_LOG_2)
def exp10(f):
    return exp2(f/_LOG10_2)

if __name__ == "__main__":
    b1 = bigfloat(3.,1)
    b2 = bigfloat(3.,-1)
    print b1*2
    print 2*b1
    print b1/2
    print 2/b1
    print b1+1
    print 1+b1
    print b1-1
    print 1-b1
    print b1**10
    print b1**100
    print -b1**1000
    print 3**b1, 3**float(b1)
    print log(b1), math.log(float(b1))
    print log10(b1), math.log10(float(b1))
    print exp2(-3)
    print exp2(3.5)