#!/usr/bin/python
from __future__ import generators

import gmpy

class cache:
    def __init__(self, f):
        self.d = {}
        self.f = f
    def __call__(self, *i):
        try:
            return self.d[i]
        except:
            r = apply(self.f, i)
            self.d[i] = r
            return r

def n_partitions(j, k = 1):
    total = 1
    j -= k
    while (j>=k):
        total += n_partitions(j, k)
        j, k = j-1, k+1
    return total

n_partitions = cache(n_partitions)

def Time(f,*l):
    import time
    start = time.time()
    x = apply(f, l)
    return (x, time.time() - start)

fact = gmpy.fac

from operator import mul,add
prod = lambda l:reduce(mul,[1L]+l)
sum = lambda l:reduce(add,[0L]+l)
multinomial = lambda l: fact(sum(l))/prod(map(fact, l))
neg_pow = lambda j: [1,-1][gmpy.getbit(j,0)]
bit_length = lambda x:gmpy.numdigits(x, 2)

binomial = gmpy.bincoef
Sqrt = gmpy.fsqrt


def B(n):
    if n == 0: return 1
    if n == 1: return gmpy.mpq(-1,2)
    if (n%2) == 1:  return 0
    return -gmpy.mpq(sum(map(lambda x:binomial(n+1, x)*B(x), range(n))), n+1)

B = cache(B)

import time, math

def Exp(x):
    prec = x.getprec()
    y = gmpy.mpf(0,prec+1)
    t = gmpy.mpf(1,prec+1)
    x = gmpy.mpf(x,prec+1)
    oy = y
    i = 1
    sign = 0
    if x < 0:
        sign = 1
        x = -x
    while(1):
        y += t
        t *= x
        t /= i
        if (oy - y) == 0:
            break
        oy = y
        i += 1
    if(sign):
        y = 1/y
    return gmpy.mpf(y, prec)

def Sin(x):
    prec = x.getprec()
    y = gmpy.mpf(0,prec+10)
    t = gmpy.mpf(x,prec+10)
    x = gmpy.mpf(x,prec+10)
    oy = y
    i = 2
    while(1):
        y += t
        t *= -x*x
        t /= i*(i+1)
        if (oy - y) == 0:
            #print "pre = %d\n" % i
            break
        oy = y
        i += 2
    return gmpy.mpf(y,prec)

def Cos(x):
    y = gmpy.mpf(0)
    t = gmpy.mpf(1)
    oy = y
    i = 1
    while(1):
        y += t
        t *= -x*x
        t /= i*(i+1)
        if (oy - y) == 0:
            #print "pre = %d\n" % i
            break
        oy = y
        i += 2
    return y

Cosh = lambda x:(Exp(x) + Exp(-x))/2
Sinh = lambda x:(Exp(x) - Exp(-x))/2
Tan = lambda x: Sin(x)/Cos(x)

def Ln_normalise(x, f, eps=0.01):
    
    # Optimization for convergence: the taylor series converges faster
    # when x is close to zero.  So a trick is to first take the square
    # root a couple of times, until x is sufficiently close to 1.
    # This function does this normalisation, passing the actual
    # evaluation on to function f(1+x)
    
    xx = gmpy.mpf(x)
    shifts=0
    precision = 2*xx.getprec()
    y = gmpy.mpf(xx,precision)
    minusone = gmpy.mpf(-1)
    x = gmpy.mpf(precision)

    delta = gmpy.mpf(eps)
    while 1:
        y = Sqrt(y)
        shifts+=1
        x = y-1
        if x < delta:
            break
    z = f(x,precision=precision)
    return z*(gmpy.mpf(2)**shifts)

def Ln_taylor(x,precision):
    # Ln(y) = Ln(1+x) = Sum(i=1 to inf) (-1)^(i+1) * x^(i) / i
    # thus y=1+x => x = y-1
    
    # i <- 0
    i = gmpy.mpf(0)
    # sum <- 1
    aResult = gmpy.mpf(0)
    # term <- 1
    term = gmpy.mpf(-1,precision);
    
    one = gmpy.mpf(1,precision);
    eps = gmpy.mpq(2)**-precision/2
    # While (term>epsilon)
    while abs(term) > eps:
        term *= -x
        i += 1
        aResult += term / i
    return aResult

Ln = lambda x:Ln_normalise(x, Ln_taylor)

def Pow(x, y):
    return Exp(y*Ln(x))
    

def Atan(x,precision=-1):
    #  1. Reduce x to positive by Atan(x) = -Atan(-x).
    #  2. According to the integer k=4t+0.25 chopped, t=x, the argument
    #      is further reduced to one of the following intervals and the
    #      arctangent of t is evaluated by the corresponding formula:
    #
    #      [0,7/16]      Atan(x) = t-t^3*(a1+t^2*(a2+...(a10+t^2*a11)...)
    #      [7/16,11/16]  Atan(x) = Atan(1/2) + Atan( (t-0.5)/(1+t/2) )
    #      [11/16.19/16] Atan(x) = Atan( 1 ) + Atan( (t-1)/(1+t) )
    #      [19/16,39/16] Atan(x) = Atan(3/2) + Atan( (t-1.5)/(1+1.5t) )
    #      [39/16,INF]   Atan(x) = Atan(INF) + Atan( -1/t )
    if precision < 0: precision = x.getprec()
    if x < 0: return -Atan(-x,precision)
    if x == 1: return Pi(precision)/4
    if x > 1: return Pi(precision)/2 + Atan(-1/x,precision)
    if x > 0.75: return Pi(precision)/4 + Atan((x-1)/(1+x),precision)
    y = gmpy.mpf(0,precision)
    t = gmpy.mpf(x,precision)
    oy = y
    i = 1
    while(1):
        y += t/i
        t *= -x*x
        if (oy - y) == 0:
            break
        oy = y
        i += 2
    return gmpy.mpf(y,precision)

def Pi(precision=-1):
    return 16*Atan(1/gmpy.mpf(5),precision) - 4*Atan(1/gmpy.mpf(239),precision)

import math

G = lambda x: fact(x-1)

def stirling(x):
    x2 = gmpy.mpf(x + 0.5)
    gi = lambda k:neg_pow(k)*kow(k)*G(2*k-1)/(((2*x2)**(2*k)/x2))
    l = (x2)*Ln(x2) + sum(map(gi, range(1, 3))) - Ln(x)/2-x
    return Exp(gmpy.mpf(l))*Sqrt(Pi())

def Plot(f, l):
    file = open("tmp", "w")
    for i in l:
        file.write("%s %s\n" % (str(i), str(f(i))))
    file.close()
    import os
    os.system("""echo "plot 'tmp'" | gnuplot -persist """)
    

if 1:
    print Sqrt(Pi())**2, Pi()

if 1:
    z = Pi()
    print "Ln (z)",Ln(z)
if 1:
    print "sin(pi) = 0 ~= ", Sin(Pi())
    print Pow(Pow(2,100), gmpy.mpf("0.01"))

if 1:
    prec = 170
    x = 0.7
    print x, Atan(gmpy.mpf(x),prec), math.atan(x)
    print x, Atan(gmpy.mpf(x),prec*2)
    print "\n\n"
    for i in range(20):
        x = i*0.1
        print x, Atan(gmpy.mpf(x),prec), math.atan(x)
    for i in range(2,20):
        x = i
        print x, Atan(gmpy.mpf(x),prec), math.atan(x)
    for i in range(2,100):
        x = 10**i
        print "10^%d" %i, Atan(gmpy.mpf(x),prec), math.atan(x)

def compute_eulerconst_besselintegral1 (prec):
    # shamelessly stolen from the ginac code :)
    # We compute f(x) classically and g(x) using the partial sums of f(x).
    actualprec = prec+1; # 1 Schutz-Digit
    intDsize = 10 # no idea what this number does?
    sx = int(0.25*0.693148*intDsize*actualprec)+1  # 1/4 * ln(2)
    N = (3.591121477*sx)   # what is this number?
    x = int(sx)**2;
    eps = gmpy.mpf(2,prec)**(-int(sx*2.88539+10))
    fterm = gmpy.mpf(1,actualprec)
    fsum = gmpy.mpf(fterm,actualprec)
    gterm = gmpy.mpf(0,actualprec)
    gsum = gmpy.mpf(gterm,actualprec)
    # After n loops
    #   fterm = x^n/n!^2, fsum = 1 + x/1!^2 + ... + x^n/n!^2,
    #   gterm = H_n*x^n/n!^2, gsum = H_1*x/1!^2 + ... + H_n*x^n/n!^2.
    for n in range(1, N):
        fterm = (fterm*x)/(n**2)
        gterm = ((gterm*x)/n + fterm)/n;
        #if (prec < 10 or n <= sx):
        #	fsum = fsum + fterm;
        #	gsum = gsum + gterm;
        #else:
        #	# For n > sx, the terms are decreasing.
        #	# So we can reduce the precision accordingly.
        fsum = fsum + fterm;
        gsum = gsum + gterm;
    result = gsum/fsum - Ln(sx)#,actualprec)
    return gmpy.mpf(result,prec) # verkürzen und fertig


for i in range(50, 1000, 50):
    print i,compute_eulerconst_besselintegral1(i)