Estimate e to "N" decimal places

e-estim.py

# Print e, (believed) correct to N decimal places.
#
# [Not absolutely guaranteed accurate, because
# I'm not an expert or professional. Casually
# believed to work based on tests with N < 3000
# with some e-decimals sources on the Internet.]
# 
# [e.g. to get identical formatting for comparison 
#  with NASA's "e.2mil", pipe through
#	 sed 's/2\./ 2\./g' | fold -60
#  I don't think my program is efficient enough to
#  churn out 2 million digits yet, but you can
#  still compare for smaller digit counts.]
#
# Can find 20000 digits in 3.451s on a 1.50 GHz/dual-core 
# "Intel(R) Celeron(R) B800" processor running 64-bit Linux.
# (This timing includes the time spent brute-forcing the
#  necessary term count). (This is probably quite lame
#  performance).
#
# The algorithm is quite simple-minded and based
# on the MacLaurin 1/n! series. It involves
# dividing and multiplying Pythonic big integers,
# which is sort of awkward.
# 
# I have plans for improvements... this should
# be much faster and less memory-intensive and
# have no more giant integer divisions, if
# all this is possible.
#
# tl;dr I love programming and calculus but am
#       not very good at either.
#	
# -- blockeduser, September 26-28, 2012.

# The error estimation formula used is based on the Lagrange 
# Remainder Term  (for finding the error in a Taylor/MacLaurin series 
# estimate). This is (essentially) what I was taught recently in 
# my Calculus class.
# 
# I found this nice web page with the formulas on it:
# http://www.millersville.edu/~bikenaga/calculus/tayerr/tayerr.html
#
#		              (n+1)
#  error for n-th term  =   f      (c)	      n + 1
#			    ---------- (x - a)
#			     (n + 1)!  
#
#	where c is "a point between x and a".
#
#				  (n+1)
# We have a = 0, f(x) = e^x, so f      (x) = e^x
# we use x = 1 to get e^1 = 1 so the error is
#
#		 e^c	     n + 1	     e^c
#	        ----- (1 - 0)		=   -----
#		(n+1)!			    (n+1)!
#
# I don't bother finding c, but since it is between 0 and 1,
# it is certain that the error is less than
#	 e^1					3
#       ------  , which is in turn less than   ---
#	(n+1)!				      (n+1)!
#
#
# My algorithm skips the first two terms of the series (1/0! + 1/1! = 2)
# to focus only on the decimal places. As a consequence, I must
# use
#
#		 	 3		3	
#	error   <=     ------ 	 =   -------
#	       	       (n+1-2)!	      (n-1)!



import sys	# (for argv)
from math import factorial, sqrt

# Brute-force approximate needed term count
# to get N digits by looking at the number
# of digits in the inverse of the error formula.
# Smarter algorithm probably exists, and I'll 
# be looking for it. This gets quite slow as n increases :(
def bf_n(dig):
	n = 1
	k = 1
	while n <= 100000:
		acc = len(str(k))
		if acc >= dig:
			return n
		for i in range(20):
			# smartass way of computing
			# successive values of n! / 3
			if n != 3:
				k *= n
			n += 1
	return None		# way too much !!!

def do_estimate_e(n, dig_limit=None):
	# under-estimate of decimal digits of accuracy,
	# based on max error estimate formula. there's
	# the logarithm way of doing this also, but
	# I don't know if it's significantly faster, and
	# things get funny:
	#	>>> log(1000) / log(10)
	#	2.9999999999999996
	acc = len(str(factorial(n - 1) / 3))

	if acc < 1:
		print "too small n !"
		return

	# the string we make
	estimate = '2.'

	# Funky algorithm based on MacLaurin theory.
	# Meant to reduce repeated computations and obtain
	# estimate-fraction fast.
	#
	# If e.g. the term count n = 5, it looks something like this:
	#
	#		1		= 1
	#		5		= 5
	#		5 * 4 		= 20
	#		5 * 4 * 3 	= 60	 -- sum so far = 86
	#		-------------------------------------------
	#		5 * 4 * 3 * 2 	= 120
	#
	#		=> e  ~  2 + 86 / 120
	#		      =
	#
	#
	#		 120
	#		----- 	= 8	 [ = (n-1)! / 3 >= error ]
	#		3 * 5    \_/
	#			  1 digit long
	#		   	      -> approx. 1 decimal digit
	#			         of accuracy

	sum = 1
	mul = n
	k = n - 1
	while k >= 2:
		sum += mul
		mul *= k
		k -= 1

	# Now we churn out the digits one by one
	# by integer-dividing quite big integers
	# and doing some tricks. [This part of the
	# program *really* needs to be redone smarter.]
	#
	# e.g. if n = 6, we started out with
	#
	#		 517
	#   e  ~   2  +  ---
	#      =	 720
	#
	# and 2 decimals of guaranteed accuracy.
	# [ (n-1)! / 3 = 5! / 3 = 40 is 2 digits long ]
	# 
	# so we do
	# 517 -> 5170
	# 5170 \ 720 = 7			*7*
	# 5170 - 7 * 720 = 130
	# 130 -> 1300
	# 1300 \ 720 = 1			*1*
	#	and now we stop because we
	#	have our guaranteed 2 decimals,
	#	beyond that could be garbage,
	#	as in the previous version of
	#	this program where my error-estimation
	#	was off by a factorial.
	#
	# [ \ is my stupid symbol for integer division ]
	#
	# For larger n, we end up doing lots of integer divisions
	# of big-ish integers, which is not very smart.

	b = sum * 10
	dc = 0
	for i in range(acc):
		# stop if we were only asked for some
		# number of digits and have them.
		if dig_limit != None and dc >= dig_limit:
			break
		dc += 1

		estimate += str(b/mul)	# add this digit to printout string
		
		# trick to go to next digit
		b -= mul * (b / mul)
		b *= 10

	info = "%d terms; accuracy: %d decimal places" % (n, acc)
	if dig_limit != None:
		info += " (%d printed)" % dig_limit
	print info
	print ""
	print estimate

if len(sys.argv) != 2:
	print "use: %s digits" % sys.argv[0]
	sys.exit(1)
else:
	digits = int(sys.argv[1])

# this can be long, so warn the user
print "figuring out needed term count..."
tc = bf_n(digits)

# do the main thing!
do_estimate_e(tc, digits)