"""real.py -- Real numbers with unlimited accuracy. Version 1.13.
By Jurjen N.E. Bos (jurjen@q2c.nl).
*** Fixed 2001 Aug 23 to work with python2.1 by K M Briggs ***

A class that support "exact" calculations by using unlimited precision
with error estimation.
For a more thorough explanation, read "real.README".

The two main functions are rep() and pa().
Arguments: a function from number of digits to result or a string that
   computes to default precision.
Result:
   rep: consecutive approximations of the result
   pa: prints all digits of result (hangs if number is exactly
                representable).

Using pa or rep to compute a given digit does not take much longer
than computing only to that precision, but is more fun.

Computations are performed using class Real.

Numbers can be created using function r():

   r(n, d, p) is n/d with p digits precision
   r(n, p) is n (string/real/int/long/float) with p digits precision
   r(n, a) is n with same precision as the Real a
   r(n) is integer n with default precision
   r(f) is float f with 53 bits precision
   r(x) is the Real x
   r(s) is value of string s

(Warning: if you increase precision of a float/Real/String, new digits
get "invented"!)

Supported operators:

   +, -, abs
   +, -, *, /, **; also both ways with int
   << and >> with int
   int(), long()

Supported mathematical functions:

   floor, ceil
   sqrt, root (n-th root), sqr (square)
   exp, log, log10
   atan, asin, acos, sin, cos, tan
   sinh, cosh, tanh, asinh, acosh, atanh
   fact (for all arguments!)
   pi and e are functions taking number of digits (with default)

These functions are also available as methods for NumPy users.
The methods have the names used there (e.g., arcsin instead of asin)

Examples: (type from real import * first)

   exp(pi()*sqrt(163))     ; almost an integer!
   exp(exp(exp(3)))        ; this is quite big
   fact(r(10000000,100))   ; so is this
   sin(exp(exp(4)))        ; gives wrong result using most libraries
                           ; (but this is correct, of course!)

   pa(pi)                  ; print pi
   pa(e)                   ; print e
   pa('(lambda d:(1-fact(d))/d)(10**-r(default/2))')
                           ; print Euler's constant (!)

   pa('sqr(fact(-1/r(2)))')
                           ; pi, using fact(-1/2)=sqrt(pi)
   pa('log10(sin(10**-r(5432,100,default)))')
   pa('log10(tan(10**-r(5432,100,default)))')
                           ; note the timing, and compare results
"""

# Version history:
# 1.12:  longShift added in atan()
# 1.11:  bug fix (1.10 doesn't even compile)
# 1.10:  mixing with float and even string supported
# 1.9:   special case for **-1,**0,**1,**2,**3
#        sqr() function disappeared again: you can use **2 now
# 1.8:   removed excessive precision of exp for small arguments
# 1.7:   functions are also methods (for NumPy users)
#        exp attribute had to be renamed to expon for this
# 1.6:   added function root(power, base)
#        added function sqr(x) = x*x
#        doesn't import itself anymore
#        bugfix: fact(r(0,n)) now gives proper precision
#        better docstrings
# 1.5:   pa() and rep() have useful second argument
# 1.4:   pa() and rep() now evaluate in global context
# 1.3:   allows for extremely high precision factorials
#        added real.doc file
# 1.2:   handles regular floating point
# 1.1:   uses _relBits to compute relative precision to get a
#        consistent use of 0
#        improved Bernoulli number computation:
#        now only does 2 gcds per number
# 1.0:   first version with factorial

import real
import math, regex, marshal, string, sys
from types import IntType, LongType, FloatType
from types import StringType, InstanceType, NoneType
IntTypes = (IntType,LongType)

###################################################
# auxiliary functions and values
###################################################

default = 32   # default number of digits

_log102 = math.log10(2) # used by repr, r and pa
_log2 = math.log(2)     # used by bitsize and log
PrecError = "PrecError" # if numbers too inaccurate

# the Mac has 53 bits floats: your mileage may vary
floatPrecision = 53

# Initialization: determine our representation length.
# Used by bitsize, but also by of _pow10, exp, log, and fact
# Uncomment the "# was:" lines, and comment the next ones,
# if you have problems, and want to use bitsize instead.

# offset = length of 1L; slope=extra bits per byte
# underEst = maximum overestimation of guess
# The approximation (len(marshal.dumps(long(x)))-_offset)*_slope
# is between underEst and 1 too low.
# We just assume they exist; otherwise, bitsize gives internal error.
_offset = len(marshal.dumps(1L))
_underEst = 1
while len(marshal.dumps(1L<<_underEst)) == _offset:
   _underEst = _underEst+1
_slope = _underEst/float(len(marshal.dumps(1L<<_underEst))-_offset)

# bitsize x = number of significant bits in x
# In other words: abs(x)>>bitsize = 0; abs(x)>>bitsize-1 = 1.
# It's a bit dirty, but it is fast, and checks its output.
def bitsize(x):
   "Number of signifant bits of an int/long"
   if x == 0: return 0
   # Make lower estimate using dumps.
   coarse = int((len(marshal.dumps(long(x)))-_offset)*_slope)
   # There are at most underEst bits left;
   # assume this fits in an int.
   x = abs(int(x>>coarse))
   fine = int(math.floor(math.log(x)/_log2))
   # We assume fine is not too high, and at most one too low.
   if x>>fine: fine = fine+1
   if x>>fine or not x>>fine-1: raise "Internal error in bitsize!"
   return coarse+fine

# Check the bitsize error estimation
for i in range(1,100):
   est = int((len(marshal.dumps(long(1L<<i-1)))-_offset)*_slope)
   if not i-_underEst<=est<i:
      raise "Unexpected long representation"
del est, i

# make sure underEst is large enough to compensate rounding errors
if _underEst<10: _underEst = 10

# To avoid the "outrageous left shift count" error, use this instead of
# the builtin << operator.
# If you have a version of Python that doesn't gives this error:
# Uncomment the "# was:" lines, and comment the next ones.
def longShift(a, b=1L):
   "longShift(a,b) = b<<a without complaints"
   while a>30000: a, b = a-30000, b<<30000
   return b<<a

# Compute integer power of 10 to given precision.
def _pow10(n, p):
   "fast integer power of 10 to given precision"
   # Compute 5**exp10 first; use extra bits for squarings.
   p = p+_underEst+bitsize(n)
   mask = 1<<bitsize(n>>1)
   # was: rm, re = 1L<<p, -p
   rm, re = longShift(p), -p
   while mask:
      # Use low estimate of bitsize.
      # was: b = 2*bitsize(rm)-p
      b = 2*int((len(marshal.dumps(rm))-_offset)*_slope)-p
      rm, re = rm*rm>>b, 2*re+b
      if n&mask: rm = rm*5
      mask = mask>>1
   return rm, re+n

# compute relative precision (in bits) of x.
# Returns default value for x zero.
def _relBits(x):
   "relative precision of a Real"
   if x.mant: return bitsize(x.mant)
   else: return int(math.ceil(default/_log102))

# Make a Real. This is a separate function from Real for speed reasons.
# has zillions of input options; see real.doc

# Pattern to match a number, for string constructor of r().
_nummatch = regex.symcomp(
   "^\(<int>[+-]?[0-9]*\)"       # integer part of mantissa
   "\(<frac>\.[0-9]*\)?"         # fractional part of mantissa
   "\(<prec>+-[0-9]\)?"          # precision indicator
   "\(<exp>[eE][+-]?[0-9]+\)?$") # exponent

def r(*args):
   "make a Real"
   if len(args) == 1:
      n, = args
      d,p = 1,None
   elif len(args) == 2:
      n,p = args
      d = 1
   elif len(args) == 3:
      n,d,p = args
      if type(n) not in IntTypes or type(d) not in IntTypes:
         raise TypeError, "Bad argument type for r()"
   else: raise TypeError, "Too many arguments for r()"
   # check if there is a precision requested
   pflag = p!=None
   # compute requested precision (with default)
   if type(p)==InstanceType and p.__class__==Real: p = _relBits(p)
   elif type(p) in (IntType,FloatType,NoneType):
      p = int(math.ceil((p or default)/_log102))
   else: raise TypeError, "Bad argument type for precision of r()"
   # find Real representation
   if type(n) == InstanceType and n.__class__==Real: pass
   elif type(n) in IntTypes:
      # special output calculation for rounded division
      sa = p-bitsize(n)+bitsize(d)
      if sa<-1: return Real ((long(n)>>-sa-1)/d+1 >>1, -sa)
      # was: else: return Real( (long(n)<<sa+1)/d+1 >>1, -sa)
      else: return Real( longShift(sa+1,long(n))/d+1 >>1, -sa)
   elif type(n) == StringType:
      if _nummatch.match(n) == -1:
         raise ValueError, "String parse error"
      intPart = _nummatch.group('int') or "0"
      fracPart = _nummatch.group('frac') or "."
      precPart = _nummatch.group('prec') or "+-0"
      expPart = _nummatch.group('exp') or "e0"
      mant = string.atol(intPart+fracPart[1:])
      exp10 = string.atoi(expPart[1:])-len(fracPart)+1
      rm, re = _pow10(abs(exp10), bitsize(mant))
      # Determine number of bits to shift away for precision.
      if exp10>0:
         if precPart == "+-0": sa = bitsize(rm)-1
         else: sa = bitsize(rm*string.atoi(precPart[2:]))
         n = Real((mant*rm>>sa-1)+1>>1, re+sa)
      else:
         if precPart == "+-0": sa = bitsize(rm)+1
         else: sa = bitsize(rm/string.atoi(precPart[2:]))
         # was: n = Real((mant<<sa)/rm+1>>1, -re-sa+1)
         n = Real(longShift(sa,mant)/rm+1 >>1, -re-sa+1)
   elif type(n) == FloatType:
      rm,re = math.frexp(n)
      if not pflag: p = floatPrecision
      return Real(long(rm*(1L<<p)),re-p)
   else: raise TypeError, "Invalid type for r()"
   # we come here only if Real or String
   # force the output precision if needed
   if pflag:
      sa = p-_relBits(n)
      if sa<-1: n = Real( (n.mant>>-sa-1)+1 >>1, n.expon-sa)
      # was: n = Real( (n.mant<<sa+1)+1 >>1, n.expon-sa)
      else: n = Real( longShift(sa+1,n.mant)+1 >>1, n.expon-sa)
   return n

##################################################
# the main class
##################################################
# The value of a Real number is anywhere in range
# (mant-1)*2**expon .. (mant+1)*2**expon
# including lower but excluding upper bound.
# Examples:
# the interval [3.136, 3.145) is given by mantissa 804 and exponent -8
# the interval [-3, -2.5) is given by mantissa -11 and exponent -2.
class Real:
   "Arbitrary precision numbers"
   def __init__(self, m, e):
      self.mant = long(m)
      self.expon = int(e)

   def __repr__(self):
      # Please fasten seatbelts before figuring this out!
      # We want an error in the range 1..5;
      # generate appropriate exponent of 10.
      exp10 = int(math.ceil((self.expon+1)*_log102))-1
      error = '+-' + `int(math.ceil(10**(self.expon*_log102-exp10)))`
      # Conversion factor for mantissa is 2**self.expon/10**exp10.
      rm, re = _pow10(abs(exp10), bitsize(self.mant))
      # was: if exp10>0: m = (self.mant<<self.expon-re+1)/rm
      if exp10>0: m = longShift(self.expon-re+1,self.mant)/rm
      else: m = self.mant*rm>>-self.expon-re-1
      sign = m<0 and '-' or ''
      m = `abs(m)+1>>1`
      if 0 >= exp10 >= -len(m)-50:
         # Cutoff final "L" and insert period in proper location.
         return sign + m[:exp10-1] + '.' + \
            string.zfill(m[exp10-1:-1],-exp10) + error
      else:
         # Use exponential notation.
         return sign + m[0] + '.' + m[1:-1] + error + \
            'e' + `exp10+len(m)-2`

   def __pos__(a): return a
   def __neg__(a): return Real(-a.mant, a.expon)
   def __abs__(a): return Real(abs(a.mant), a.expon)

   def __add__(a,b):
      if type(b)==InstanceType and b.__class__==Real:
         # b is a Real: make sure a is less accurate than b.
         if a.expon < b.expon: a,b = b,a
         s = a.expon-b.expon
         m = a.mant + (b.mant>>s)
         # If m odd, we can't represent the sum: enlarge interval 1 bit.
         if m&1: return Real(m>>2, a.expon+2)
         else: return Real(m>>1, a.expon+1)
      elif type(b) in IntTypes: # Int or long have infinite precision.
         if a.expon <= 0:
            # Maximum output precision is twice the original precision.
            expmin = bitsize(b)-2*_relBits(a)
            if expmin < a.expon:
               # We have no precision problems here.
               # was: return Real(a.mant+(long(b)<<-a.expon), a.expon)
               return Real(a.mant+longShift(-a.expon,long(b)), a.expon)
            elif expmin <= 0:
               # We have to round a off.
               return Real(
                  (a.mant>>expmin-a.expon)+longShift(-expmin,long(b)),
                  expmin)
            else:
               # We have to round both off.
               return Real(
                  ((a.mant>>-a.expon)+long(b)>>expmin-1)+1>>1,
                  expmin)
         # exact right shift
         # was: if a.expon <= bitsize(b) and b>>a.expon<<a.expon == b:
         if a.expon <= bitsize(b) and \
               longShift(a.expon,b>>a.expon) == b:
            return Real(a.mant+(b>>a.expon), a.expon)
         return Real(a.mant+(b>>a.expon)+1>>1, a.expon+1)
      else: return a+r(b)

   def __mul__(a,b):
      if type(b)==InstanceType and b.__class__==Real:
         # b is a Real: make sure a is less accurate than b.
         al, bl = _relBits(a), _relBits(b)
         if al > bl: a,al,b,bl = b,bl,a,al
         # Error of product: abs(a.mant)+abs(b.mant>>bl-al) < 2*(1<<al).
         m = a.mant*(b.mant>>bl-al)>>al
         # It is now at most 2 to high, 3 too low.
         return Real(m+2>>2, a.expon+b.expon+bl+2)
      elif type(b) in IntTypes: # Int or long have infinite precision.
         if b == 0: return 0
         bs = bitsize(b)
         # was: if 1L<<bs-1 == abs(b):
         if longShift(bs-1) == abs(b):
            return Real((b>>bs-1)*a.mant, a.expon+bs-1)
         # Shifted product is at most 1 too high, 2 too low.
         return Real((a.mant*b >> bs)+1>>1, a.expon+bs+1)
      else: return a*r(b)

   def __div__(a,b):
      if type(b)==InstanceType and b.__class__==Real:
         # b is a Real
         al, bl = bitsize(a.mant), bitsize(b.mant)
         if bl<1: raise PrecError, "Division by possible zero"
         # Shift amount to get result of correct precision.
         sa = min(bl-1, 2*bl-al-2)
         # Very seldom case: a too precise; reduce precision.
         if sa<0:
            a = Real(a.mant>>-sa, a.expon-sa)
            sa = 0
         # Error of quotient:
         # (1<<sa)/b.mant + (a.mant<<sa)/(b.mant*(b.mant-1)) < 1+1;
         # this may be 3 too low or 2 too high.
         # was: return Real((a.mant<<sa)/b.mant+2>>2,
         return Real(longShift(sa,a.mant)/b.mant+2>>2,
            a.expon-b.expon-sa+2)
      elif type(b) in IntTypes: # Int or long have infinite precision.
         if b == 0: raise ValueError, "Division by zero"
         bs = bitsize(b)
         # was: if 1L<<bs-1 == abs(b):
         if longShift(bs-1) == abs(b):
            return Real((b>>bs-1)*a.mant, a.expon-bs+1)
         # Quotient is at most 2 too high, 3 too low.
         # was: return Real(((a.mant<<bs)/b+2)>>2, a.expon-bs+2)
         return Real( longShift(bs,a.mant)/b+2 >>2, a.expon-bs+2)
      else: return a/r(b)

   def __rdiv__(a,b):
      if type(b) in IntTypes: # Int or long have infinite precision.
         sa = 2*bitsize(a.mant)-2
         if sa<0: raise PrecError, "Division by possible zero"
         # Error of quotient is (b<<sa)/(a.mant*(a.mant-1)) < 1.
         # was: return Real((long(b)<<sa)/a.mant+1>>1, -sa-a.expon+1)
         return Real( longShift(sa,long(b))/a.mant+1 >>1, -sa-a.expon+1)
      else: return r(b)/a

   def __radd__(a,b): return a+b
   def __sub__(a,b): return a+-b
   def __rsub__(a,b): return -a+b
   def __rmul__(a,b): return a*b

   def __pow__(a,b):
      if b in [-1,0,1,2,3]:
         # optimize trivial cases
         if b==-1: return 1/a
         else: return reduce(lambda x,y:x*y,[a]*b,1)
      else: return exp(log(a)*b)

   def __rpow__(a,b):
      if type(b) in IntTypes:
         # Set b to same relative precision as a (this is plenty).
         sa = bitsize(a.mant)-bitsize(b)
         # was: if sa>=0: b = Real(long(b)<<sa,-sa)
         if sa>=0: b = Real( longShift(sa,long(b)) ,-sa)
         else: b = Real(b>>-sa,-sa)
         return exp(a*log(b))
      else: return r(b)**a

   def __lshift__(a,b):
      "fast multiply by power of 2"
      return Real(a.mant,a.expon+b)
   def __rshift__(a,b):
      "fast divide by power of 2"
      return Real(a.mant,a.expon-b)

   def __long__(a):
      "convert to long. Mostly truncated, sometimes rounded"
      if a.expon>=0: raise PrecError, "Not enough accuracy for int/long"
      # Python says the long() operator may round or truncate.
      # That's exactly what we're going to do:
      # It will be truncated most of the time, but sometimes rounded!
      return a.mant>>-a.expon

   def __int__(a): return int(long(a))

   def __hex__(a):
      "Binary floating point representation in hex"
      bs = bitsize(a.mant)
      ae = a.expon+bs
      head = a.mant<0 and "-0x." or "0x."
      mant = hex(abs(a.mant)<<-bs%4)[2:-1]
      err = "+-" + `1<<-bs%4`
      if ae>0: exp = "<<" + hex(ae)[2:]
      elif ae<0: exp = ">>" + hex(-ae)[2:]
      else: exp = ""
      # why does hex() do longs in upper case, and ints in lower case?
      return head + mant + err + string.upper(exp)

############################################
# the standard functions
############################################

def floor(x):
   "largest integral number smaller than x"
   # Number is very inaccurate: floor doesn't do much.
   if x.expon>=0: return x
   # Check if it almost is an integer.
   r = x.mant>>-x.expon
   # was: if r<<-x.expon == x.mant: return Real(r,0)
   if longShift(-x.expon,r) == x.mant: return Real(r,0)
   # Now we can make a long, since we know the exact value.
   return r

# Why this works, is left as an exercise for the reader!
# Hint: what if x is an exact integer?
ceil = lambda x: floor(x)+1

# square root
# -----------

def sqrt(x):
   "square root"
   if type(x) != InstanceType: x = r(x)   # force a Real
   if x.mant<-1: raise ValueError, "Root of negative number"
   elif x.mant<1: raise PrecError, "Root of possible negative number"
   # Output exponent, generating resulting error of at most 1.
   outExp = x.expon-bitsize(x.mant)>>1
   rm = _sqrtAux(x.mant, x.expon, outExp)
   # Due to rounding, we now are 2 too low, 1 too high.
   return Real(rm+1>>1, outExp+1)

# Auxiliary function that does square root (xm,re)
# with output exponent re.
def _sqrtAux(xm, xe, re):
   "square root with given output exponent"
   xl = bitsize(xm)
   #nnumber of bits used for estimate must be even.
   estBits = min(xl-1,(xe+xl>>1)-re,64)&-2
   rm = long(math.sqrt(xm>>xl-estBits-(xe+xl&1))*(1L<<estBits/2))
   # Shift mantissa for division.
   sa = xe-2*re
   # was: if sa>0: m = xm<<sa
   if sa>0: m = longShift(sa,xm)
   else: m = xm>>-sa
   # Start up Newton's iteration.
   # was: rm,rold = rm<<(xe+xl>>1)-re-estBits,-1
   rm,rold = longShift((xe+xl>>1)-re-estBits,rm),-1
   while not rold-1 <= rm <= rold:
      rold,rm = rm,m/rm+rm>>1
   return rm

def root(power,base):
   "root to any base. root(x,y) is x-th root of y"
   return exp(log(base)/power)

# exponential functions
# ---------------------

def exp(x):
   "power of e"
   if type(x) != InstanceType: x = r(x)   # Force a Real.
   if x.expon>0: raise PrecError, "Not enough accuracy for exp"
   # prec is our output exponent: it is kept low for x near 0.
   prec = min(-x.expon,2*_relBits(x))
   # Make sure x is small enough to guarantee fast convergence.
   # First term: make x<1; second term: heuristic speedup.
   d = bitsize(x.mant>>-x.expon)+int(math.sqrt(bitsize(x.mant)/2))
   # Now represent x/2**d as xm,-xs.
   xs = prec+_underEst
   if xs+x.expon>0: xm = x.mant<<xs+x.expon
   else: xm = x.mant>>-xs-x.expon
   xs = xs+d
   # Compute exp(x/2**d) with the Taylor series.
   # was: rs, term, n = 1L<<xs, xm, 2
   rs, term, n = longShift(xs), xm, 2
   while not -1 <= term <= 0:
      rs,term,n = rs+term, (term*xm>>xs)/n, n+1
   re = -xs
   # Now we have exp(x/2**d) as rs,re.
   # Square for the halving of the input; keep bitsize about xs.
   while d:
      # Use low estimate of bitsize.
      # was: b = 2*bitsize(rs)-xs
      b = 2*int((len(marshal.dumps(rs))-_offset)*_slope)-xs
      rs, re = rs*rs>>b, re*2+b
      d = d-1
   # Determine shift amount of result: make result precision prec-1.
   b = bitsize(rs)-prec+1
   # Now we are 2 too low, 1 too high.
   return Real((rs>>b)+1>>1, re+b+1)

def log(x):
   "natural logarithm"
   if type(x) != InstanceType: x = r(x)   # Force a Real.
   xl = bitsize(x.mant)
   if x.mant<0: raise ValueError, "Log of negative number"
   if x.mant<1: raise PrecError, "Log of possible negative number"
   # First estimate the result; represent est as em,-ee.
   em = long((
      math.log(x.mant<<20>>xl)+_log2*(xl+x.expon-20)
      )*0x100000)
   # Now compute exponent of estimate;
   # make sure est is smaller than 1/4 to speed up convergence.
   # (This 1/4 improves speed a bit.)
   d = bitsize(em>>20)+2
   ee = 20+d
   # Compute exp(est/2**d) using xl+_underEst+d bits.
   re = max(ee,xl+_underEst+d)
   # was: rm, term, n = 1L<<re, em<<re-ee, 2
   rm, term, n = longShift(re), longShift(re-ee,em), 2
   while not -1 <= term <= 0:
      rm, term, n = rm+term, (term*em>>ee)/n, n+1
   # Now we have exp(est/2**d) as rm,-re;
   # square for the halving of the input.
   bs = re
   while d:
      # Use low estimate of bitsize.
      # was: b = 2*bitsize(rm)-bs
      b = 2*int((len(marshal.dumps(rm))-_offset)*_slope)-bs
      rm, re = rm*rm>>b, re*2-b
      d = d-1
   # Now we have r=exp(est) as rm,-re;
   # compute the log of x/r in precision of result.
   # First compute x/r-1 = (x-r)/r as xrm,-bs.
   # was: xrm = ((x.mant<<x.expon+re)-rm <<bs)/rm
   xrm = longShift(bs,longShift(x.expon+re,x.mant)-rm)/rm
   # was: rm, term, n = em<<bs-20, xrm, 1
   rm, term, n = longShift(bs-20,em), xrm, 1
   while not -1 <= term <= 0:
      rm, term, n = rm+term/n, -term*xrm>>bs, n+1
   # Now we can return: log(x)=log(x/exp(est))+log(exp(est))=
   # =log(x/rs)+est.
   # Output exponent based on input precision: -xl+2.
   return Real((rm>>bs-xl+1)+1 >> 1, 2-xl)

def log10(x):
   "base 10 logarithm"
   if type(x) != InstanceType: x = r(x)   # force a Real
   b = bitsize(x.mant)+bitsize(x.expon)+5
   # was: return log(x)/log(Real(10L<<b, -b))
   return log(x)/log(Real(longShift(b,10L), -b))

# trigonometric functions
# -----------------------

def atan(x):
   "arc tangent"
   if type(x) != InstanceType: x = r(x)   # force a Real
   xl = bitsize(x.mant)
   # Heuristic optimization: determine upper limit for x.
   limit = int(math.sqrt(xl))/5
   # The loop below reduces x using tangent of half angle
   # xs=-xe is the output precision: absolute precision of x,
   # plus twice the bits before the point.
   # Limit it to twice the relative precision of x.
   if xl+x.expon<=0:
      # If abs(x)<1, prepare for loop immediately.
      # Precision: abs. of x, maximum twice rel. precision.
      xs = -x.expon+limit+10
      xm = longShift(xs+x.expon,x.mant)
      halfs = 0
   else:
      # Do one reduction before entering loop to get x below 1.
      # Determine precision.
      xs = 2*xl+min(x.expon,0)+limit+10
      # Tangent of half angle: x/(sqrt(x**2+1)+1).
      x1m, x1e = x.mant*x.mant<<limit+10, 2*x.expon-limit-10
      # was: if x1e<=0: x1m = x1m+(1L<<-x1e)
      if x1e<=0: x1m = x1m+longShift(-x1e)
      x1m = _sqrtAux(x1m,x1e,-xs)
      # was: x1m = x1m+(1L<<xs)
      x1m = x1m+longShift(xs)
      # was: xm = (x.mant<<2*xs+x.expon)/x1m
      xm = longShift(2*xs+x.expon,x.mant)/x1m
      halfs = 1
   # We now have x in xm,-xs. Keep xs steady from here.
   # was: one, one1, xlimit = 1L<<xs, 1L<<2*xs, 1L<<xs-limit
   one,one1,xlimit = longShift(xs),longShift(2*xs),longShift(xs-limit)
   # Reduce x below limit using half angle formula.
   while not -xlimit < xm < xlimit:
      # Tangent of half angle: x/(sqrt(x**2+1)+1)
      x1m = _sqrtAux(xm*xm+one1, -2*xs, -xs)
      # was: xm = (xm<<xs)/(x1m+one)
      xm = longShift(xs,xm)/(x1m+one)
      halfs = halfs+1
   # Now compute tangent of the reduced number.
   rs, term, n = 0L,xm, 1
   xm = xm*xm>>xs
   while term:
      rs,term,n = rs + term/n, -term*xm>>xs, n+2
   # Now we are 2 too low, 1 too high.
   # Double result for each halving.
   return Real((rs>>limit+10-halfs)+1>>1, -xs+limit+11)

asin = lambda x: atan(x/sqrt(1-x*x))

acos = lambda x: atan(sqrt(1-x*x)/x)

def _sincos(x):
   "Auxiliary function for sin(), cos()"
   # Take x mod pi for trigonometry functions.
   # returns x%pi.mant, x%pi.expon, rounding bits, oddflag.
   if type(x) != InstanceType: x = r(x)   # force a Real
   # Precision is more than a half circle: return nothing informative.
   if x.expon > 2: return 0, 0, 1, 0
   # First, compute pi/2 to correct precision (see pi()).
   xl = max(bitsize(x.mant), -x.expon)
   # was: p = 2*atan(Real(1L<<xl+10,-xl-10))
   p = 2*atan(Real(longShift(xl+10),-xl-10))
   if p.expon > x.expon: raise "Internal error in sin/cos"
   # Take x mod pi, nearest to zero.
   xm, pm = x.mant<<x.expon-p.expon, p.mant
   q = xm/pm+1>>1
   # Precision: same as x.
   return xm-(q<<1)*pm, -p.expon, x.expon-p.expon, q&1

def sin(x):
   "sine (in radians)"
   rm, re, rb, f = _sincos(x)
   # Now use the series to approximate sin.
   rs, term, n = 0L, rm, 1
   rm = rm*rm>>re
   while term:
      rs, term, n = rs+term, -term*rm/((n+1)*(n+2))>>re, n+2
   if f: rs = -rs
   return Real((rs>>rb)+1>>1,-re+rb+1)

def cos(x):
   "cosine (in radians)"
   rm, re, rb, f = _sincos(x)
   # Now use the series to approximate cos.
   rm = rm*rm>>re
   # was: rs, term, n = 1L<<re, -rm>>1, 2
   rs, term, n = longShift(re), -rm>>1, 2
   while term:
      rs, term, n = rs+term, -term*rm/((n+1)*(n+2))>>re, n+2
   if f: rs = -rs
   return Real((rs>>rb)+1>>1,-re+rb+1)

tan = lambda x: sin(x)/cos(x)

# hyperbolic functions
# --------------------
def sinh(x):
   "Hyperbolic sine (implemented using exp())"
   ex = exp(x)
   return (ex-1/ex)>>1

def cosh(x):
   "Hyperbolic cosine (implemented using exp())"
   ex = exp(x)
   return (ex+1/ex)>>1

def tanh(x):
   "Hyperbolic tangent (implemented using exp())"
   ex = exp(x<<1)
   return (ex-1)/(ex+1)

asinh = lambda x: log(x+sqrt(x*x+1))
acosh = lambda x: log(x+sqrt(x*x-1))
atanh = lambda x: log((1+x)/(1-x))>>1

# factorial function
# ------------------
# This is the most complicated function of all!

# auxiliary function: greatest common divisor
def gcd(a,b):
   "greatest common divisor of int/long"
   a,b = abs(a),abs(b)
   while a: a,b = b%a,a
   return b

# auxiliary function: bernoulli numbers for even arguments
# cached output values
_bernCache = [(1L,1L)]
# least common multiple of denominators of bernCache (and 2 for bern(1))
_bernLCM = 2L

# compute the Bernoulli number B(n)
def bern(m):
   "Compute Bernoulli number as (numerator,denominator)"
   if m&1:
      if m==1: return -1,2
      else: return 0,1
   global _bernCache, _bernLCM
   # fill the cache until enough
   while 2*len(_bernCache) <= m:
      n = 2*len(_bernCache)
      # since sum(k,choose(n,k)*B(k))=B(n) for n>1:
      # compute sum(k=0,n-2,choose(n+1,k)*B(k))/(n+1)
      n1choosek = 1L       # keeps choose(n+1,k)
      # start with the only nonzero odd bernoulli number
      sum = -_bernLCM*(n+1)/2
      # no do the even terms
      for k in range(0,n,2):
         tn,td = _bernCache[k/2]
         sum = sum + n1choosek*(_bernLCM/td*tn)
         # compute next binomial coefficient
         # use only one long multiplication and division
         n1choosek = n1choosek*((n+1-k)*(n-k))/((k+1)*(k+2))
      # reduce result
      bd = _bernLCM*(n+1)
      g = gcd(bd,sum)
      sum,bd = -sum/g,bd/g
      _bernCache.append((sum,bd)) # KMB was sum,bd
      g = gcd(bd,_bernLCM)
      _bernLCM = _bernLCM*bd/g
   # return the last value added
   # if cached: just return the cached value
   return _bernCache[m/2]

# the factorial function!
def fact(x):
   "factorial: fact(x) = gamma(x+1)"
   if type(x) != InstanceType: x = r(x)
   # compute value of pi to correct precision
   # we need this in two places
   myPi = atan(Real(longShift(-x.expon),x.expon))<<2
   if int(x)<-100:
      # use mirror formula: x!*(-x)! = pi*x/sin(pi*x)
      pix = myPi*x
      return pix/(sin(pix)*fact(-x))
   # determine minimun value for stirling
   # decrease factor below (e.g., 0.5) for faster repeated computations
   # may not be lower than .12 (infinite loop results!)
   # increase (e.g., 1.1) to save memory, and faster first computation
   minx = int(-x.expon*0.9)+1
   # compute number of division steps and the start value
   steps = max(minx-int(x),0)
   minx = x+steps
   # compute the stirling formula
   res = log(2*myPi*minx)/2+(log(minx)-1)*minx
   # now sum the terms. We use fixed precision for speed and accuracy
   # as precision, we use the precision of x,
   # +underEst bits extra for last loop
   # keep this number for later
   auxBits = res.expon-minx.expon+_underEst
   rm,rs = res.mant<<auxBits,_underEst-minx.expon
   # compute sum(k>0,B(k+1)/(k*(k+1)*x**k))) with k shifted
   k = 0
   # keep x**(k+1) in xpow
   xpow = minx.mant<<_underEst
   x2 = xpow*xpow>>rs
   t = 1L   # last term: to check convergence
   while t:
      k = k+2
      n,d = bern(k)
      # was: t = (n<<rs*2)/(d*((k-1)*k)*xpow)
      # we inline longShift here for speed
      ns,sh = n,2*rs
      while sh>30000:ns,sh = ns<<30000,sh-30000
      t = (ns<<sh)/(d*((k-1)*k)*xpow)
      rm = rm+t
      xpow = x2*xpow>>rs
   res = exp(Real(rm,-rs))
   # now convert from minx! to x! by division
   rm,re = res.mant,res.expon
   xm,xe = minx.mant, minx.expon
   # auxiliary values to speed up loop
   xs,one = max(bitsize(minx.mant)-_underEst,0),longShift(-xe)
   for i in range(steps):
      sa = rs-int((len(marshal.dumps(rm))-_offset)*_slope)+xs
      rm,re = (rm<<sa)/xm, re-sa-xe
      xm = xm-one
   return Real( (rm>>auxBits+9)+1 >>1,re+auxBits+10)

############################################
# constants
############################################

# compute pi using atan
def pi(*args):
   "Compute Real pi with given number of digits"
   if len(args) == 0: p = default
   else: (p,) = args
   return 4*atan(r(1,p+1))

def e(*args):
   "Compute Real e with given number of digits"
   if len(args) == 0: p = default
   else: (p,) = args
   return exp(r(1,p+2))

############################################
# user interface
############################################

# Function to repeatedly evaluate an expression with more and more
# digits (function is evaluated with growing values of default).
# You can also pass a function that takes as (only) argument the
# number of digits to use for the calculation.

# Computing any amount of digits using rep takes less than twice the
# time to compute only that precision, so you don't waste too much time.
def rep(s, p=default):
   "Evaluate expression with increasing precisions"
   global default
   save_default = default
   if type(s) == StringType:
      import __main__
      sc = compile(s, "<rep()>", "eval")
   try:
      while p<100000000:   # you won't get this far
         try:
            if type(s) == StringType:
               default = p
               print eval(sc,__main__.__dict__)
            else: print s(p)
         except PrecError, err: print "(" + err + ")"
         except KeyboardInterrupt: return
         p = -(-p*4/3)     # round up
   finally: default = save_default

# The function to print "all" digits of a number!
# Input is the same as rep.
def pa(s, p=default):
   "Compute 'all' digits of expression"
   global default
   if type(s) == StringType:
      import __main__
      sc = compile(s, "<pr()>", "eval")
   noExpon = 1
   save_default = default
   try:
      # First step: compute enough digits to figure out the exponent.
      while noExpon:
         try:
            if type(s) == StringType:
               default = p
               x = eval(sc,__main__.__dict__)
            else: x = s(p)
            expon = int(ceil(log10(abs(x))))
         except PrecError: p = p*2
         except KeyboardInterrupt: return
         else: noExpon = 0
      if x.mant<0: sys.stdout.write('-')
      if expon < -100:
         sys.stdout.write("10**"+`expon`+"*\n.")
      elif expon < 0: sys.stdout.write("."+"0"*-expon)
      elif expon > 100:
         sys.stdout.write(
            "(Decimal point after "+`expon`+" digits.)\n")
      sys.stdout.flush()
      # Now print the mantissa.
      try:
         printedDigits, allDigits, theDigits = 0, 0, 0
         while p < 10000000:
            # Compute power of exponent.
            rm, re = _pow10(abs(expon), bitsize(x.mant))
            rr = Real((rm>>9)+1>>1, re+10)
            # Compute the mantissa.
            if expon>0: mantissa = abs(x)/rr
            elif expon<0: mantissa = abs(x)*rr
            else: mantissa = abs(x)
            # We have a mantissa: let's compute some digits.
            # Estimate the number of digits we can print.
            newDigits = int(math.floor(bitsize(mantissa.mant)*_log102))
            while newDigits > printedDigits:
               # Use previously computed digits for comparison.
               digitsPower = 10L**(newDigits-printedDigits)
               oldDigits = allDigits
               try:
                  # Compute the new digits.
                  allDigits = long(floor(mantissa*10L**newDigits))
                  # Subtract previously computed digits.
                  theDigits = allDigits - digitsPower*oldDigits
                  if not 0<=theDigits<digitsPower:
                     raise "Internal error in pa()"
                  # All OK: print them!
                  outData = string.zfill(
                     `theDigits`[:-1],newDigits-printedDigits)
                  if printedDigits<=expon<newDigits:
                     # Put a period at the proper location.
                     sys.stdout.write(outData[:expon-printedDigits])
                     sys.stdout.write('.')
                     sys.stdout.write(outData[expon-printedDigits:])
                  else: sys.stdout.write(outData)
                  sys.stdout.flush()
               except PrecError:
                  # There is still doubt about the digits: try fewer.
                  newDigits = newDigits - 1
               else: printedDigits = newDigits
            # We can't print any more: compute more digits of x.
            p = -(-p*4/3)     # round up
            if type(s) == StringType:
               default = p
               x = eval(sc,__main__.__dict__)
            else: x = s(p)
      except KeyboardInterrupt:
         sys.stdout.write('\n')
         return
   finally: default = save_default

# Now make methods for the functions to help NumPy users.
# The names are in correspondence with the "umath" module.
# These do not always match the function names in the math library!
Real.floor,Real.ceil=floor,ceil
Real.sqrt=sqrt
Real.exp,Real.log,Real.log10=exp,log,log10
Real.arctan,Real.arcsin,Real.arccos=atan,asin,acos
Real.sin,Real.cos,Real.tan=sin,cos,tan
Real.sinh,Real.cosh,Real.tanh=sinh,cosh,tanh
# These aren't in the umath library (yet)
Real.arcsinh,Real.arccosh,Real.arctanh=asinh,acosh,atanh