In [1]:

    

%matplotlib inline


    

In [2]:

    

import sympy
from sympy import log, exp, integrate, EulerGamma, pi, symbols, Function, simplify
from sympy.functions.special.gamma_functions import digamma, gamma
from sympy.plotting import plot


    

In [3]:

    

x = symbols('x')
F = Function('F')
f = digamma(-1/x) + EulerGamma.evalf() # - digamma(1).evalf()


    

In [4]:

    

plot(f, (x, 1.7377, 5))


    

    













    
Out[4]:





<sympy.plotting.plot.Plot at 0x109b5b250>

In [5]:

    

plot(f, (x, -10, -1e-5))


    

    













    
Out[5]:





<sympy.plotting.plot.Plot at 0x10dded290>

In [6]:

    

plot(f, (x, -100, -10))


    

    













    
Out[6]:





<sympy.plotting.plot.Plot at 0x109b5b3d0>

In [7]:

    

[[i, v, v2,  v - v2] for i in xrange(2,25) for v, v2 in [(f.evalf(subs={x: i}), f.evalf(subs={x: i/(1.0-i)}))]]


    

    
Out[7]:





[[2, 0.613705638880109, -1.38629436111989, 2.00000000000000],
 [3, 2.25898124911494, -0.741018750885056, 3.00000000000000],
 [4, 3.49135478511506, -0.508645214884939, 4.00000000000000],
 [5, 4.61220709819539, -0.387792901804606, 5.00000000000000],
 [6, 5.68648625222927, -0.313513747770728, 6.00000000000000],
 [7, 6.73681978717086, -0.263180212829140, 7.00000000000000],
 [8, 7.77319859335384, -0.226801406646162, 8.00000000000000],
 [9, 8.80072726463218, -0.199272735367815, 9.00000000000000],
 [10, 9.82228871495448, -0.177711285045519, 10.0000000000000],
 [11, 10.8396352530268, -0.160364746973197, 11.0000000000000],
 [12, 11.8538937363563, -0.146106263643719, 12.0000000000000],
 [13, 12.8658219386562, -0.134178061343823, 13.0000000000000],
 [14, 13.8759483056153, -0.124051694384721, 14.0000000000000],
 [15, 14.8846526483315, -0.115347351668485, 15.0000000000000],
 [16, 15.8922149904058, -0.107785009594165, 16.0000000000000],
 [17, 16.8988463322028, -0.101153667797221, 17.0000000000000],
 [18, 17.9047086908058, -0.0952913091941544, 18.0000000000000],
 [19, 18.9099285390088, -0.0900714609912065, 19.0000000000000],
 [20, 19.9146060486962, -0.0853939513037924, 20.0000000000000],
 [21, 20.9188215915782, -0.0811784084218070, 21.0000000000000],
 [22, 21.9226404021718, -0.0773595978281745, 22.0000000000000],
 [23, 22.9261159820042, -0.0738840179957958, 23.0000000000000],
 [24, 23.9292926244895, -0.0707073755105217, 24.0000000000000]]

In [8]:

    

f.evalf(subs={x:1e+15}) / 1e+15


    

    
Out[8]:





0.999999996372506

In [9]:

    

plot(f, (x, -3, -0.5))


    

    













    
Out[9]:





<sympy.plotting.plot.Plot at 0x10df9e250>

In [10]:

    

digamma(1)


    

    
Out[10]:





-EulerGamma

In [11]:

    

digamma(1).evalf()


    

    
Out[11]:





-0.577215664901533

In [12]:

    

f.evalf(subs={x:-1})


    

    
Out[12]:





4.94291515243064e-18

In [13]:

    

gamma(1)


    

    
Out[13]:





1

In [14]:

    

[[y, y.diff(), simplify(y.diff(x)/y), simplify(log(y).diff(x))] for y in [-1/x, 1/(1-x), x/(1-x), 1-1/x]]


    

    
Out[14]:





[[-1/x, x**(-2), -1/x, -1/x],
 [1/(-x + 1), (-x + 1)**(-2), -1/(x - 1), -1/(x - 1)],
 [x/(-x + 1), x/(-x + 1)**2 + 1/(-x + 1), -1/(x*(x - 1)), -1/(x*(x - 1))],
 [1 - 1/x, x**(-2), 1/(x*(x - 1)), 1/(x*(x - 1))]]

In [15]:

    

[[y, simplify(exp(1j*pi + integrate(log(y).diff(x), x)))] for y in [-1/x, 1/(1-x), x/(1-x), 1-1/x]]


    

    
Out[15]:





[[-1/x, -1/x],
 [1/(-x + 1), -1/(x - 1)],
 [x/(-x + 1), -x/(x - 1)],
 [1 - 1/x, -(x - 1)/x]]

In [16]:

    

[[y, exp(integrate(simplify(log(1/(1-y)).diff(x)), x))] for y in [-1/x, 1/(1-x), x/(1-x), 1-1/x]]


    

    
Out[16]:





[[-1/x, x/(x + 1)],
 [1/(-x + 1), (x - 1)/x],
 [x/(-x + 1), (x - 1)/(x - 1/2)],
 [1 - 1/x, x]]

In [17]:

    

log(-1).evalf()


    

    
Out[17]:





3.14159265358979*I

In [18]:

    

[[y, simplify(exp(log(-1) + integrate(log(y).diff(x), x)))] for y in [-1/x, 1/(1-x), x/(1-x), 1-1/x]]


    

    
Out[18]:





[[-1/x, -1/x],
 [1/(-x + 1), -1/(x - 1)],
 [x/(-x + 1), -x/(x - 1)],
 [1 - 1/x, (-x + 1)/x]]

In [19]:

    

exp(integrate(simplify(log(integrate(F(-1/x), x)).diff(x)), x))


    

    
Out[19]:





exp(Integral(F(-1/x)/Integral(F(-1/x), x), x))

In [20]:

    

h = simplify(f - x)


    

In [21]:

    

plot(h, (x, -1e7, -1e5))


    

    













    
Out[21]:





<sympy.plotting.plot.Plot at 0x10e6a3f90>

In [22]:

    

plot(h, (x, 1e5, 1e7))


    

    













    
Out[22]:





<sympy.plotting.plot.Plot at 0x10e179110>

In [23]:

    

[[i, v, v - i, v - v2 - i] for i in xrange(int(1e6),int(1e6+25)) for v, v2 in [(f.evalf(subs={x: i}), f.evalf(subs={x: i/(1.0-i)}))]]


    

    
Out[23]:





[[1000000, 999999.999998355, -1.64494849741459e-6, 0],
 [1000001, 1000000.99999836, -1.64494849741459e-6, 0],
 [1000002, 1000001.99999836, -1.64494849741459e-6, 0],
 [1000003, 1000002.99999836, -1.64494849741459e-6, 0],
 [1000004, 1000003.99999836, -1.64494849741459e-6, 0],
 [1000005, 1000004.99999836, -1.64494849741459e-6, 0],
 [1000006, 1000005.99999836, -1.64494849741459e-6, 0],
 [1000007, 1000006.99999836, -1.64494849741459e-6, 0],
 [1000008, 1000007.99999836, -1.64494849741459e-6, 0],
 [1000009, 1000008.99999836, -1.64494849741459e-6, 0],
 [1000010, 1000009.99999836, -1.64494849741459e-6, 0],
 [1000011, 1000010.99999836, -1.64494849741459e-6, 0],
 [1000012, 1000011.99999836, -1.64494849741459e-6, 0],
 [1000013, 1000012.99999836, -1.64494849741459e-6, 0],
 [1000014, 1000013.99999836, -1.64494849741459e-6, 0],
 [1000015, 1000014.99999836, -1.64494849741459e-6, 0],
 [1000016, 1000015.99999836, -1.64494849741459e-6, 0],
 [1000017, 1000016.99999836, -1.64494849741459e-6, 0],
 [1000018, 1000017.99999836, -1.64494849741459e-6, 0],
 [1000019, 1000018.99999836, -1.64494849741459e-6, 0],
 [1000020, 1000019.99999836, -1.64494849741459e-6, 0],
 [1000021, 1000020.99999836, -1.64494849741459e-6, 0],
 [1000022, 1000021.99999836, -1.64494849741459e-6, 0],
 [1000023, 1000022.99999836, -1.64494849741459e-6, 0],
 [1000024, 1000023.99999836, -1.64494849741459e-6, 0]]

In [24]:

    

f.subs(x, 2)


    

    
Out[24]:





-2*log(2) - EulerGamma + 2.57721566490153

In [25]:

    

f.subs(x, 2).evalf()


    

    
Out[25]:





0.613705638880109

In [26]:

    

[[i, f.subs(x, i)] for i in xrange(10)]


    

    
Out[26]:





[[0, polygamma(0, zoo) + 0.577215664901533],
 [1, zoo],
 [2, -2*log(2) - EulerGamma + 2.57721566490153],
 [3, -3*log(3)/2 - EulerGamma + sqrt(3)*pi/6 + 3.57721566490153],
 [4, -3*log(2) - EulerGamma + pi/2 + 4.57721566490153],
 [5, 0.577215664901533 + polygamma(0, -1/5)],
 [6, 0.577215664901533 + polygamma(0, -1/6)],
 [7, 0.577215664901533 + polygamma(0, -1/7)],
 [8, 0.577215664901533 + polygamma(0, -1/8)],
 [9, 0.577215664901533 + polygamma(0, -1/9)]]

In [27]:

    

[[i, f.subs(x, i)] for i in xrange(0,-10,-1)]


    

    
Out[27]:





[[0, polygamma(0, zoo) + 0.577215664901533],
 [-1, -EulerGamma + 0.577215664901533],
 [-2, -2*log(2) - EulerGamma + 0.577215664901533],
 [-3, -3*log(3)/2 - sqrt(3)*pi/6 - EulerGamma + 0.577215664901533],
 [-4, -3*log(2) - pi/2 - EulerGamma + 0.577215664901533],
 [-5, polygamma(0, 1/5) + 0.577215664901533],
 [-6, polygamma(0, 1/6) + 0.577215664901533],
 [-7, polygamma(0, 1/7) + 0.577215664901533],
 [-8, polygamma(0, 1/8) + 0.577215664901533],
 [-9, polygamma(0, 1/9) + 0.577215664901533]]

In [28]:

    

plot(f-x, (x, 1e5, 1e7))


    

    













    
Out[28]:





<sympy.plotting.plot.Plot at 0x10e165510>

In [29]:

    

plot(f.subs(x, x/(1-x)), (x, 1e5, 1e7))


    

    













    
Out[29]:





<sympy.plotting.plot.Plot at 0x10eb52b90>

In [30]:

    

plot(f-x-f.subs(x, x/(1-x)), (x, 1e5, 1e7))


    

    













    
Out[30]:





<sympy.plotting.plot.Plot at 0x10e417110>

In [31]:

    

EulerGamma.evalf()


    

    
Out[31]:





0.577215664901533

In [32]:

    

plot(f-x-f.subs(x, x/(1-x)), (x, 1.1, 9))


    

    













    
Out[32]:





<sympy.plotting.plot.Plot at 0x10ef69490>

In [33]:

    

plot(f-f.subs(x, x/(1-x)), (x, 1.1, 9))


    

    













    
Out[33]:





<sympy.plotting.plot.Plot at 0x10ef3d890>

In [34]:

    

plot(f-f.subs(x, x/(1-x)), (x, -9, 9))


    

    













    
Out[34]:





<sympy.plotting.plot.Plot at 0x10f0f5890>

In [35]:

    

plot(f.subs(x, x/(1-x)), (x, 1e3, 1e6))


    

    













    
Out[35]:





<sympy.plotting.plot.Plot at 0x10ef69510>

In [36]:

    

plot(f.subs(x, x/(1-x)), (x, -1e2, -1))


    

    













    
Out[36]:





<sympy.plotting.plot.Plot at 0x10f1cb710>

In [37]:

    

f.subs(x, x/(1-x)).subs(x, -1).evalf()


    

    
Out[37]:





1.00000000000000

In [38]:

    

f.subs(x, x/(1-x)).subs(x, 1).evalf()


    

    
Out[38]:





zoo

In [39]:

    

f.subs(x, -1).evalf()


    

    
Out[39]:





4.94291515243064e-18

In [40]:

    

f.subs(x, 1).evalf()


    

    
Out[40]:





zoo

In [41]:

    

simplify((-1/x).subs(x, x/(1-x)))


    

    
Out[41]:





(x - 1)/x

In [42]:

    

h = digamma(-1/x) - digamma((x-1)/x) - x


    

In [43]:

    

[simplify(h.subs(x, i)) for i in xrange(2, 25)]


    

    
Out[43]:





[0,
 0,
 0,
 -5 - polygamma(0, 4/5) + polygamma(0, -1/5),
 -6 - polygamma(0, 5/6) + polygamma(0, -1/6),
 -7 - polygamma(0, 6/7) + polygamma(0, -1/7),
 -8 - polygamma(0, 7/8) + polygamma(0, -1/8),
 -9 - polygamma(0, 8/9) + polygamma(0, -1/9),
 -10 - polygamma(0, 9/10) + polygamma(0, -1/10),
 -11 - polygamma(0, 10/11) + polygamma(0, -1/11),
 -12 - polygamma(0, 11/12) + polygamma(0, -1/12),
 -13 - polygamma(0, 12/13) + polygamma(0, -1/13),
 -14 - polygamma(0, 13/14) + polygamma(0, -1/14),
 -15 - polygamma(0, 14/15) + polygamma(0, -1/15),
 -16 - polygamma(0, 15/16) + polygamma(0, -1/16),
 -17 - polygamma(0, 16/17) + polygamma(0, -1/17),
 -18 - polygamma(0, 17/18) + polygamma(0, -1/18),
 -19 - polygamma(0, 18/19) + polygamma(0, -1/19),
 -20 - polygamma(0, 19/20) + polygamma(0, -1/20),
 -21 - polygamma(0, 20/21) + polygamma(0, -1/21),
 -22 - polygamma(0, 21/22) + polygamma(0, -1/22),
 -23 - polygamma(0, 22/23) + polygamma(0, -1/23),
 -24 - polygamma(0, 23/24) + polygamma(0, -1/24)]

In [44]:

    

[simplify(h.subs(x, i)).evalf() for i in xrange(2, 25)]


    

    
Out[44]:





[0,
 0,
 0,
 0.e-125,
 0.e-124,
 0.e-124,
 -0.e-126,
 -0.e-126,
 0.e-125,
 0.e-125,
 -0.e-126,
 0.e-125,
 0.e-125,
 0.e-125,
 0.e-127,
 -0.e-125,
 -0.e-127,
 -0.e-127,
 0.e-126,
 -0.e-127,
 -0.e-126,
 -0.e-125,
 -0.e-125]

In [45]:

    

def newton(h, x, x0, eps=1e-15):
    c = x0
    while True:
        if abs(h.evalf(subs={x: c})) < eps:
            break
        c += -(h / h.diff(x)).evalf(subs={x: c})
    return c


    

In [46]:

    

[[r, fr] for r in [newton(f, x, 2)] for fr in [f.evalf(subs={x: r})]]


    

    
Out[46]:





[[1.76256875808632, -2.54808341503461e-16]]

In [47]:

    

[[r, fr] for r in [newton(f - EulerGamma.evalf(), x, 2)] for fr in [f.evalf(subs={x: r}) - EulerGamma.evalf()]]


    

    
Out[47]:





[[1.98380025433306, -3.33066907387547e-16]]

In [48]:

    

[f.evalf(subs={x: -1}), 0, (digamma(1) + EulerGamma).evalf(), 0]


    

    
Out[48]:





[4.94291515243064e-18, 0, 0, 0]

In [49]:

    

[f.evalf(subs={x: -1.0/2}), 1, (digamma(2) + EulerGamma).evalf(), 1]


    

    
Out[49]:





[1.00000000000000, 1, 1.00000000000000, 1]

In [50]:

    

(x/(1-x)).integrate(x)


    

    
Out[50]:





-x - log(x - 1)

In [51]:

    

exp((x/(1-x)).integrate(x))


    

    
Out[51]:





exp(-x)/(x - 1)

In [52]:

    

simplify((exp(-x)/(x - 1)).diff() / exp(-x)/(x - 1))


    

    
Out[52]:





(-x - (x - 1)**2 + 1)/(x - 1)**4

In [53]:

    

simplify((-x - log(x - 1) + 123456789).diff(x))


    

    
Out[53]:





-x/(x - 1)

In [54]:

    

simplify((-x - (x - 1)**2 + 1)/(x - 1)**4 - -x/(x - 1))


    

    
Out[54]:





(x*(x - 1)**3 - x - (x - 1)**2 + 1)/(x - 1)**4

In [55]:

    

simplify(exp(pi*1j + (-1/x).integrate(x)))


    

    
Out[55]:





-1/x

In [56]:

    

simplify(exp(pi*1j + (x/(1-x)).integrate(x)))


    

    
Out[56]:





exp(-x + 1.0*I*pi)/(x - 1)

In [57]:

    

simplify(exp(pi*1j + (1/(1-x)).integrate(x)))


    

    
Out[57]:





-1/(x - 1)

In [58]:

    

simplify(exp(pi*1j + (x**2).integrate(x)))


    

    
Out[58]:





-exp(x**3/3)

In [ ]: