I've always been fascinated by how to compute high precission approximations to the more interesting functions such as sin, exponential and the gamma functions. Here are some simple methods of computing approximations to these functions. Thanks to Victor Kowalenko for invaluable contributions!
The most important part is to import a library to access the gmp multiprecission library:
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)
Notice how I use dynamic programming to accellerate the computation? apparently the number of partitions grows as e^(sqrt(n))
def Time(f,*l):
import time
start = time.time()
x = apply(f, l)
return (x, time.time() - start)
Although I'm really only interested in how to compute functions, it is useful to know which of two methods are faster so that I can build new functions on old functions.
fact = gmpy.fac from operator import mul,add prod = lambda l:reduce(mul,l,1) sum = lambda l:reduce(add,l,0) 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
These functions are interesting to implement in their own right, but I'm leaving that to the experts. I've defined sum and product to correctly handle the empty case. Multinomial probably isn't the most efficient method possible, but it is well fast enough. neg_pow is a convienience function to compute -1^k, I could possibly have just used pow (that would be log(n) time, but neg_pow is easier to type :) Bit length is merely an approximation to lg(x). Binomial and Sqrt are already defined in gmp.
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)
This is a relatively fast way to compute Bernoulli numbers, basically derived from the definition. The approach is to rearrange terms of the Bernoulli polynomials in terms of the previous values and solve (they all end up being linear functions). I have cached the intermediate values, which speeds thing up immensely. Also notice that all the odd Bs are zero.
Ok, onto some interesting functions!
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)
Probably the easiest function to evaluate - the series converges everywhere, and usually quite fast. One subtlety is the fact that for negative values it converges much more slowly. To get around this I first tranform the input (negate it), then evaluate the positive version, finally reciprocating the answer. Victor says it may converge slowly for large x, I haven't tried.
Ok, I found that for very large values, Exponential does indeed take a long time to evaluate:
Time(Exp(1000000)) =
(mpf('3.03321539680208754509e434294',64), 34.884858965873718)
That's nearly 35 seconds! The problem is that we have to wait until the terms get to around 1000000 before the series makes its convergence. Before that, each term in the series is growing.
I realised that I could take out a factor of 2^n quite cheaply and then compute the value by squaring repeatedly - with this modification Exp's evaluation time becomes only really dependant on the requested precission and the size of the output. Exp2 source
Time(Exp2(1000000)) =
(Number of squares taken out = 20)
(mpf('3.03321539680208754509e434294',64), 0.0014319419860839844)
Only 1.4ms for a number with almost half a million digits :)
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
Again, the power series is a very good starting point for Sin and Cosine. These simply implement these. I should use a reduction method for large magnitude x, as the series becomes less and less stable with increasing x. But to do so would require me to compute an accurate version of Pi.
Tan is pretty obvious: Tan = lambda x: Sin(x)/Cos(x)
Cosh = lambda x:(Exp(x) + Exp(-x))/2 Sinh = lambda x:(Exp(x) - Exp(-x))/2 Tanh = lambda x: Sinh(x)/Cosh(x)
This is the first function to present some difficulty - there is a power series, but it converges very slowly, and only for |x-1| < 1. The 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)
def Ln_normalise(x, f, eps=0.01):
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)
Raising a number to an integer power is easy - you just multiply enough times, but how do you compute to the power 0.7? The exponential function provides the key - with a bit of rearranging we get:
def Pow(x, y): return Exp(y*Ln(x))
It's not a perfect solution because it would be nice if, e.g., 256^(1/4) gave 16 exactly, but we're doing floating point, not number theory, so I'm leaving that as an exercise to the reader :)
Arctangent was tricky because the standard power series converges really slowly much past 1/2. Looking in some libraries for doing simply floating point I found this remark:
# 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 )
Really helpful, except it doesn't say how to get those constants...
def Atan(x,precision=-1):
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)
Except we don't know how to compute pi?! Well Tan(pi/4) = 1, so maybe we can use that? Well it works, but it is so slow that it takes a few seconds to compute something as accurate as my calculator!
It turns out that someone(should relook his name up) discovered, via a reduction formula, another expression for Pi using Atan of much smaller numbers, which converge fast enough to be useful.
def Pi(precision=-1):
return 16*Atan(1/gmpy.mpf(5),precision) - 4*Atan(1/gmpy.mpf(239),precision)
You can see that I can use my Atan formula to compute Pi, and use Pi to compute Atans!
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())
This is still very much a work in progress, apparently there aren't many good ways to work out the gamma function, but this is an implementation of Stirling's approximation. The terms inside the sum are based on the Bernoulli numbers but I've used another function, kow. I'll fix this up sometime...
This is a simple function to plot a function using gnuplot.
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 """)
def compute_eulerconst_besselintegral1 (prec):
# shamelessly stolen from the cln 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)