Once I could print out iCal Tasks using Python, I though I’d move on to marking tasks as complete using Python. The following program will take a list of Task UIDs as parameters on the command line, and mark the corresponding tasks as complete:

#!/usr/bin/python
from CalendarStore import CalCalendarStore
import sys

if len(sys.argv) < 2:
    print "Syntax:nt%s  ..." % (sys.argv[0])
    sys.exit(0)

store = CalCalendarStore.defaultCalendarStore()

for taskUID in sys.argv[1:]:
    task = store.taskWithUID_(taskUID)
    task.setIsCompleted_(True)
    store.saveTask_error_(task, None)
    print "%s marked complete" % (task.title())

This program attempts to retrieve each Task from the CalendarStore, mark each as complete, and print out a simple message. This does assume that all of the UIDs are valid, and skips the requisite error checking you’d need to make this a robust and useful tool, but the idea is there. The output of the program looks like:

./todoComplete.py A771147B-3F98-4895-8CB1-8466A7433DD0
Task #3 marked complete

Digging around i PyObjC and the CalendarStore (which is the local shared calendar database for which ever user is logged in), I thought it would be fun to pull the Tasks from my ToDo list into Python. This small chunk of code will extract all available of your tasks from all of your calendars:

#!/usr/bin/python
from CalendarStore import CalCalendarStore

store = CalCalendarStore.defaultCalendarStore()
calendars = CalCalendarStore.calendars(store)
predicate = CalCalendarStore.taskPredicateWithCalendars_(calendars)
tasks = store.tasksWithPredicate_(predicate)

for task in tasks:
    completedMark = " "
    completeDate = task.completedDate()
    if completeDate is not None:
        completedMark = "#"
    print "%s(%s:%s) %s" % (completedMark, task.calendar().title(), task.priority(), task.title())
    print "   uid: %s" % (task.uid())

The end result is a list of tasks:

#(Home:5) Look, a task!
   uid: 1AD57B3A-C31A-4BBD-BCB7-19E78E4060B9
 (Home:1) Task #1
   uid: 240C7449-B74B-4666-98DF-1A6C6F6564CF
 (Work:5) Task #2
   uid: A36AFE1B-5BC6-4A3A-B18F-7B5DB4E4A5F3
 (Home:9) Task #3
   uid: A771147B-3F98-4895-8CB1-8466A7433DD0

The values in parentheses are the Calendar from which the task is taken, and the task priority. Tasks without a priority have a value of 0, priority of High is 1, Medium is 5, and Low is 9. The text following the parenthesized values is the title of the task, below which is the Task’s UID, which is used by the CalendarStore to reference specific tasks. Any tasks that have a completedDate of None are not yet marked as complete in iCal. Those tasks that are complete are printed with a hash (#) sign in front of them.

Among various other tasks for a project this week I’ve been writing some Python classes to query Apple’s AddressBook database. OS X 10.5 ships with PyObjC, a Python to Objective-C bridge, preinstalled, which makes this kind of task fairly straight forward. Or, at least, it should. A few minutes of digging discovered that what few examples of working with Python and AddressBook were easy to find were incomplete, broken, or counter to the task I was trying to accomplish. Add to that a lack of any good documentation on the PyObjC methods, and it fell to experimenting and digging through Objective-C documentation.



What follows is a method for retrieving a list of People from the local AddressBook as a list of Python Dictionaries.

from AddressBook import *
import pprint

def addressBookToList():
        """
        Read the current user's AddressBook database, converting each person
        in the address book into a Dictionary of values. Some values (addresses,
        phone numbers, email, etc) can have multiple values, in which case a
        list of all of those values is stored. The result of this method is
        a List of Dictionaries, with each person represented by a single record
        in the list.
        """
        # get the shared addressbook and the list of
        # people from the book.
        ab = ABAddressBook.sharedAddressBook()
        people = ab.people()

        peopleList = []

        # convert the ABPerson to a hash
        for person in people:
                thisPerson = {}
                props = person.allProperties()
                for prop in props:

                        # skip some properties
                        if prop == "com.apple.ABPersonMeProperty":
                            continue
                        elif prop == "com.apple.ABImageData":
                            continue

                        # How we convert the value depends on the ObjC
                        # class used to represent it
                        val = person.valueForProperty_(prop)
                        if type(val) == objc.pyobjc_unicode:
                                # Unicode String
                                thisPerson[prop.lower()] = val
                        elif issubclass(val.__class__, NSDate):
                                # NSDate
                                thisPerson[prop.lower()] = val.description()
                        elif type(val) == ABMultiValueCoreDataWrapper:
                                # List -- convert each item in the list
                                # into the proper format
                                thisPerson[prop.lower()] = []
                                for valIndex in range(0, val.count()):
                                        indexedValue = val.valueAtIndex_(valIndex)
                                        if type(indexedValue) == objc.pyobjc_unicode:
                                                # Unicode string
                                                thisPerson[prop.lower()].append(indexedValue)
                                        elif issubclass(indexedValue.__class__, NSDate):
                                                # Date
                                                thisPerson[prop.lower()].append(indexedValue.description())
                                        elif type(indexedValue) == NSCFDictionary:
                                                # NSDictionary -- convert to a Python Dictionary
                                                propDict = {}
                                                for propKey in indexedValue.keys():
                                                        propValue = indexedValue[propKey]
                                                        propDict[propKey.lower()] = propValue
                                                thisPerson[prop.lower()].append(propDict)
                peopleList.append(thisPerson)
        return peopleList

So, given an entry in the AddressBook that looked like:




This method will have the following Dictionary in the list of returned People:

{   u'address': [   {   u'city': u'Anytown',
                        u'country': u'USA',
                        u'countrycode': u'us',
                        u'state': u'NY',
                        u'street': u'123 Fake Street',
                        u'zip': u'10111'}],
    u'aiminstant': [u'john_doe_aim'],
    u'creation': u'2008-04-09 13:33:17 -0500',
    u'email': [u'john@doe.com'],
    u'first': u'John',
    u'last': u'Doe',
    u'modification': u'2008-04-09 13:33:17 -0500',
    u'phone': [u'555-555-1212'],
    u'uid': u'BBFAB17F-591A-4D3C-BE75-4FE5B25B984D:ABPerson'}

The Python Decimal class is a handy arbitrary precision numeric type for dealing with extraordinarily large or small numbers, or numbers requiring exact precision. The strength of this class is in it’s use of overridden methods that allow it to be used in standard mathematical constructs:

>>> Decimal(2) ** 16 - 1
Decimal("65535")

Because of these overridden methods, the class can in fact be passed to many of the math library’s methods:

>>> from math import *
>>> floor(Decimal("1.21"))
1.0
>>> fmod(Decimal("10.25"), Decimal("10"))
0.25
>>> modf(Decimal("1.234"))
(0.23399999999999999, 1.0)

Of course the short coming of this is what is returned: a built in float type. When using the math library methods we lose the precision of the Decimal type. We could attempt to use the math library log method to calculate the logarithm of a Decimal type:

>>> log(Decimal(100), 10)
2.0
>>> log(Decimal(100), 3.14)
4.0247145803329714

For a log that returns an integer this works well, but for anything requiring decimal precision we lose the accuracy and precision of the Decimal type. We can instead use an adaptation of the common logarithm, reconstructed to work as a class method for the Decimal type. The Python code:

def dec_log(self, base=10):
	cur_prec = getcontext().prec
	getcontext().prec += 2
	baseDec = Decimal(10)
	retValue = self

	if isinstance(base, Decimal):
		baseDec = base
	elif isinstance(base, float):
		baseDec = Decimal("%f" % (base))
	else:
		baseDec = Decimal(base)

	integer_part = Decimal(0)
	while retValue < 1:
		integer_part = integer_part - 1
		retValue = retValue * baseDec
	while retValue >= baseDec:
		integer_part = integer_part + 1
		retValue = retValue / baseDec

	retValue = retValue ** 10
	decimal_frac = Decimal(0)
	partial_part = Decimal(1)
	while cur_prec > 0:
		partial_part = partial_part / Decimal(10)
		digit = Decimal(0)
		while retValue >= baseDec:
			digit += 1
			retValue = retValue / baseDec
		decimal_frac = decimal_frac + digit * partial_part
		retValue = retValue ** 10
		cur_prec -= 1
	getcontext().prec -= 2

	return integer_part + decimal_frac

This algorithm calculates the integer and fractional portions of the logarithm in separate passes, composing the two when returning the value. Taking a cue from a comment by Jake Voytko on adding nth roots to Python, this algorithm adjusts the working precision of the Decimal types by adding 2 extra digits of precision during the calculation, and trimming two digits from the precision when the calculation is complete, to ensure that the trailing digits of the result are as accurate as possible (a technique I should have used in previous examples, had I consulted the Decimal numeric recipes).

The logarithm can be added to the Decimal type and used with no parameters, a float, an integer or a Decimal parameter. Without a parameter the logarithm defaults to log₁₀:

>>> Decimal.log = dec_log
>>> Decimal(100).log()
Decimal("2.000000000000000000000000000000000000000000000000000000000000000000000
000000000000000000000000000000")
>>> Decimal(100).log(3.14)
Decimal("4.024714580332970585367941961553706972817670933108784336660631220853887
020695465320530995672957013298")
>>> Decimal(1234).log(Decimal("1.1234"))
Decimal("61.17246795014182187657829372589717334008147684454521224410597688872638
875290237619255864354194950937")

The simplest means of working with large numbers in Python is using the Decimal type, which supports exact representation of numbers with arbitrary precision. We’ve already used the Decimal type when working with π, γ, and ζ.

The Decimal type is especially useful as it overrides all of the basic mathematical operators, so instead of using method calls to perform operations on a Decimal number, you can use normal operators:

>>> Decimal(2) ** 2
Decimal("4")
>>> 1 / Decimal(4)
Decimal("0.25")

One of the built in methods for the Decimal type is sqrt, which generates the square root of the current Decimal value. But what if you want the cube root, or the fifth root, or even the 713^th root of a number? We can add generic roots to the Decimal class with an n-th root algorithm.

This root finding algorithm is fast converging, so we can get results fairly quickly, and the method is fairly simple. For a given n-th root:
ⁿ√A
Where we want the n^th root of A, we can use the following pattern:

1. Make an initial guess at the value, which we’ll call x₀
2. At every iteration we build on the previous values for x, such that x_{k + 1} = 1 / n * ( (n – 1) * x_k + A / x_k^{n – 1})
3. Repeat step 2 until the last iteration and this iteration have the same result, at which point we’ve reached the proper value in the current precision

The algorithm looks like this in Python:

from decimal import *

def nth_root_gen(num, n, digits, iter):
	getcontext().prec = digits
	a = Decimal(num)
	oneOverN = 1 / Decimal(n)
	nMinusOne = Decimal(n) - 1
	curVal = Decimal(num) / (Decimal(n) ** 2)
	if curVal <= Decimal("1.0"):
		curVal = Decimal("1.1")
	lastVal = 0
	curN = Decimal(1)
	while (lastVal != curVal) and (curN <= iter):
		lastVal = curVal
		curVal = oneOverN * ( (nMinusOne * curVal) + (a / (curVal ** nMinusOne)))
		yield curVal
		curN = curN + 1

def nth_root(num, n, digits):
	getcontext().prec = digits
	a = Decimal(num)
	oneOverN = 1 / Decimal(n)
	nMinusOne = Decimal(n) - 1
	curVal = Decimal(num) / (Decimal(n) ** 2)
	if curVal <= Decimal("1.0"):
		curVal = Decimal("1.1")
	lastVal = 0
	while lastVal != curVal:
		lastVal = curVal
		curVal = oneOverN * ( (nMinusOne * curVal) + (a / (curVal ** nMinusOne)))
	return curVal

def dec_ext_nth_root(self, root):
	a = self
	oneOverN = 1 / Decimal(root)
	nMinusOne = Decimal(root) - 1
	curVal = self / (Decimal(root) ** 2)
	if curVal <= Decimal("1.0"):
		curVal = Decimal("1.1")
	lastVal = 0
	while lastVal != curVal:
		lastVal = curVal
		curVal = oneOverN * ( (nMinusOne * curVal) + (a / (curVal ** nMinusOne)))
	return curVal

Two of these methods can be used directly in code or in the Python interpreter, nth_root_gen and nth_root. The former is a generator that yields it’s current value at every iteration, so you can monitor the value as it converges on the actual root value:

>>> for r in nth_root_gen(100, 7, 20, 100):
...     print r
...
1.9470038113262390670
1.9311017464441574735
1.9306979823757723882
1.9306977288833500144
1.9306977288832501670
1.9306977288832501670

or using the nth_root method:

>>> print nth_root(100, 7, 20)
1.9306977288832501670

With both of these we’re calculating the 7^th root of 100, to an accuracy of 20 decimal places. In the generator version we also specify a maximum number of iterations to use before giving up the calculation, in this case 100.

The third method, dec_ext_nth_root, has a considerably different method signature. That’s because it’s designed to extend the built in Decimal class, by adding a generic nth_root method. It can be used like this:

Decimal.nth_root = dec_ext_nth_root

That only needs to be executed once for any script that wants to extend the Decimal class. After that, you can use the nth_root method attached to any Decimal object:

>>> Decimal.nth_root = dec_ext_nth_root
>>> a = Decimal(125)
>>> a.sqrt()
Decimal("11.180339887498948482")
>>> a.nth_root(2)
Decimal("11.180339887498948482")
>>> a.nth_root(3)
Decimal("5.0000000000000000000")
>>> a.nth_root(15)
Decimal("1.3797296614612148324")
>>> a.nth_root(120)
Decimal("1.0410563801900299289")

So we’ve fairly easily added a fast converging method to the Python Decimal object to handle arbitrary roots. There are, of course, some limitations: the algorithm we use depends on the Decimal classes exponentiation operator, which only supports whole number powers, which means that we an only calculate whole number roots. This means that we can find the 3^rd root, but not the 3.14^th root. This doesn’t effect the number we’re finding the root of, however. Any number that can be represented by the Decimal class we can take the root of, but the root type must be a whole number. We’re also limited to positive roots; we can’t calculate the -3^rd root.

Building on the code developed for calculating the Riemann zeta function, we can write a fairly simple set of functions for finding the value of γ, the Euler-Mascheroni Constant. Not to be confused with e, γ is used in number theory and numeric analysis problems.

It’s interesting to calculate γ because it’s a very slowly convergent number, and it’s built on top of a function (in our implementation) that is also slowly convergent. The function we’ll use for calculating gamma is:

γ = Σ where n = [2 -> ∞) of (-1ⁿ * ζ(n)) / n

When choosing how accurate a result of γ to calculate, we'll need to take into account how accurate we need ζ(n) calculated. If we attempt too many digits of precision for ζ(n), each iteration of γ could take minutes, which would drastically slow our calculation of an already slowly converging number.

Our Python code to calculate γ:

#!/usr/bin/python

from decimal import *
import time

def zeta(n, digits, iter):
	getcontext().prec = digits
	curSum = Decimal(0)
	curN = Decimal(1)
	while curN <= iter:
		curSum = curSum + (1 / (curN ** n))
		curN = curN + 1
	return curSum

def gamma_gen(digits, iter, zeta_iter):
	getcontext().prec = digits
	curSum = Decimal(0)
	curN = Decimal(2)
	negone = Decimal(-1)
	while curN < iter:
		curSum += ((negone ** curN) * zeta(curN, digits, zeta_iter)) / curN
		yield curSum
		curN += 1

def gamma(digits, iter, zeta_iter):
	getcontext().prec = digits
	curSum = Decimal(0)
	curN = Decimal(2)
	negone = Decimal(-1)
	while curN < iter:
		curSum += ((negone ** curN) * zeta(curN, digits, zeta_iter)) / curN
		curN += 1
	return curSum

def find_gamma(digits, iter, zeta_iter):
	g_const = "0.57721566490153286060"
	looking_at = 1
	last_g = Decimal(0)
	curiter = 1
	start = time.time()
	for g in gamma_gen(digits, iter,zeta_iter):
		while (g.to_eng_string()[:looking_at + 2] == g_const[:looking_at + 2]) and (g.to_eng_string()[:looking_at + 2] == last_g.to_eng_string()[:looking_at + 2]):
			curtime = time.time()
			print “Found %d digits at iteration %d (%f seconds)” % (looking_at, curiter, curtime – start)
			looking_at += 1
		curiter += 1
		last_g = g

We've left out the zeta_gen method because we're only using the non-generator version of ζ calculation. Two gamma methods are listed, one a generator, the other a direct calculation. When calling either method, we specify how many digits after the decimal point our numbers should have, how many iterations we should use to calculate γ, and how many iterations each calculation of ζ should use. The higher the ζ iterations, the slower our γ values are generated.

The methods can be called in the following manner, using 100 digits of precision, with 10,000 iterations to create γ, with 300 iterations to generate ζ values (which will give us approximately 5 digits of accuracy for our ζ values):

y = gamma(100, 10000, 300)

for y in gamma_gen(100, 10000, 300):
    print y

A third method is included, which will track the generation of γ and report as the accuracy of the calculated value increases. How accurately γ can be calculated is heavily dependent on how accurately each ζ value is calculated. The lower the number of iterations used when finding ζ values, the fewer accurate digits of γ that can be found.

The find_gamma method can be called in a similar manner to the other two gamma methods:

find_gamma(100, 10000, 300)

Knowing that the value of γ is approximately 0.57721 56649 01532 86060, we can test various values for our iterations using the following Python code combined with the gamma function defined above:

realy = Decimal("0.57721566490153286060")
yiters = [10, 50, 100, 200, 250, 1000]
ziters = [10, 50, 100, 200, 250, 1000]
for yiter in yiters:
	for ziter in ziters:
		st = time.time()
		y = gamma(20, yiter, ziter)
		et = time.time()
		yaccuracy = (Decimal("1.0") - ((realy - y) / realy)) * 100
		print "gamma(20, %d, %d)" % (yiter, ziter)
		print "    y = %s" % (y.to_eng_string())
		print "    %s%% accurate" % (yaccuracy.to_eng_string())
		print "    %f seconds to calculate" % (et - st)

This test generates these values:

gamma(20, 10, 10)
    y = 0.47851665659521345415
    82.900843773330986011% accurate
    0.078000 seconds to calculate
gamma(20, 10, 50)
    y = 0.51482338102399920668
    89.190819364166572838% accurate
    0.422000 seconds to calculate
gamma(20, 10, 100)
    y = 0.51970067621726072567
    90.035788669373065474% accurate
    0.813000 seconds to calculate
gamma(20, 10, 200)
    y = 0.52216971548126864071
    90.463538540719524532% accurate
    1.718000 seconds to calculate
gamma(20, 10, 250)
    y = 0.52266598571969530717
    90.549515112147350409% accurate
    2.203000 seconds to calculate
gamma(20, 10, 1000)
    y = 0.52415975665402424336
    90.808304161918506284% accurate
    8.750000 seconds to calculate
gamma(20, 50, 10)
    y = 0.52097300115391053324
    90.256213202873356618% accurate
    0.641000 seconds to calculate
gamma(20, 50, 50)
    y = 0.55727972558912639506
    96.546188794822930761% accurate
    3.719000 seconds to calculate
gamma(20, 50, 100)
    y = 0.56215702078238791914
    97.391158100029424279% accurate
    7.578000 seconds to calculate
gamma(20, 50, 200)
    y = 0.56462606004639583418
    97.818907971375883337% accurate
    15.562000 seconds to calculate
gamma(20, 50, 250)
    y = 0.56512233028482250064
    97.904884542803709214% accurate
    19.750000 seconds to calculate
gamma(20, 50, 1000)
    y = 0.56661610121915143684
    98.163673592574865091% accurate
    84.251000 seconds to calculate
gamma(20, 100, 10)
    y = 0.52604798241963353033
    91.135430724904526274% accurate
    1.594000 seconds to calculate
gamma(20, 100, 50)
    y = 0.56235470685484939215
    97.425406316854100416% accurate
    9.750000 seconds to calculate
gamma(20, 100, 100)
    y = 0.56723200204811091623
    98.270375622060593935% accurate
    18.812000 seconds to calculate
gamma(20, 100, 200)
    y = 0.56970104131211883127
    98.698125493407052992% accurate
    38.875000 seconds to calculate
gamma(20, 100, 250)
    y = 0.57019731155054549773
    98.784102064834878869% accurate
    48.891000 seconds to calculate
gamma(20, 100, 1000)
    y = 0.57169108248487443393
    99.042891114606034746% accurate
    205.783000 seconds to calculate
gamma(20, 200, 10)
    y = 0.52856673124800451833
    91.571792553165140269% accurate
    3.953000 seconds to calculate
gamma(20, 200, 50)
    y = 0.56487345568322038015
    97.861768145114714412% accurate
    23.266000 seconds to calculate
gamma(20, 200, 100)
    y = 0.56975075087648190423
    98.706737450321207931% accurate
    43.766000 seconds to calculate
gamma(20, 200, 200)
    y = 0.57221979014048981927
    99.134487321667666988% accurate
    92.032000 seconds to calculate
gamma(20, 200, 250)
    y = 0.57271606037891648573
    99.220463893095492865% accurate
    116.563000 seconds to calculate
gamma(20, 200, 1000)
    y = 0.57420983131324542193
    99.479252942866648742% accurate
    477.159000 seconds to calculate
gamma(20, 250, 10)
    y = 0.52906898120188240022
    91.658805083215508989% accurate
    4.860000 seconds to calculate
gamma(20, 250, 50)
    y = 0.56537570563709826204
    97.948780675165083132% accurate
    28.719000 seconds to calculate
gamma(20, 250, 100)
    y = 0.57025300083035978612
    98.793749980371576650% accurate
    56.578000 seconds to calculate
gamma(20, 250, 200)
    y = 0.57272204009436770116
    99.221499851718035707% accurate
    117.501000 seconds to calculate
gamma(20, 250, 250)
    y = 0.57321831033279436762
    99.307476423145861584% accurate
    152.220000 seconds to calculate
gamma(20, 250, 1000)
    y = 0.57471208126712330382
    99.566265472917017462% accurate
    718.733000 seconds to calculate
gamma(20, 1000, 10)
    y = 0.53057273117000842393
    91.919322955401557198% accurate
    35.735000 seconds to calculate
gamma(20, 1000, 50)
    y = 0.56687945560522428575
    98.209298547351131341% accurate
    204.599000 seconds to calculate
gamma(20, 1000, 100)
    y = 0.57175675079848580983
    99.054267852557624859% accurate
    422.542000 seconds to calculate
gamma(20, 1000, 200)
    y = 0.57422579006249372487
    99.482017723904083916% accurate
    696.751000 seconds to calculate
gamma(20, 1000, 250)
    y = 0.57472206030092039133
    99.567994295331909793% accurate
    951.662000 seconds to calculate
gamma(20, 1000, 1000)
    y = 0.57621583123524932753
    99.826783345103065670% accurate
    3445.229000 seconds to calculate

Some interesting selections from these numbers:

• The simplest of the calculations (gamma(20, 10, 10)) is already over 80% accurate.
• 90% accuracy is only slightly more complex, still taking less than a second (gamma(20, 10, 100))
• The jump to 95% accurate is almost immediate (gamma(20, 50, 50)), and takes only a little under 4 seconds.
• 99% accuracy happens first at gamma(20, 100, 1000), but again at gamma(20, 200, 200). The second occurrence takes less than half as long to calculate.
• The difference in 99.134% accuracy (gamma(20, 200, 200)) and 99.826% accuracy (gamma(20, 1000, 1000)) is 3352 seconds of calculation.
• Even the most accurate calculation is still only precise to two digits (0.57)

So to reach better than 90% accuracy takes only a second of calculations, but getting close to 100% accuracy is an exercise in diminishing returns.

Today I’m working with the ζ function, also known as the Riemann zeta function. ζ shows up in a lot of mathematical constructs, and a few of the constants and formulas I’d like to work with depend on it for their calculations. The general form of the ζ function should be familiar, it’s a summation:

ζ(s) = Σ where n = [1 -> ∞] of 1/n^s

There is one important difference from some of our earlier work with summations (Working with π, Working with e); the Riemann zeta function converges very slowly to an approximate value. Whereas our π calculation could generate 2000 decimals of accuracy in a few minutes, the zeta function can take well over a minute to reach 8 or 9 decimal places of accuracy. Take for instance Apery’s constant, which is defined as ζ(3). The approximate value of ζ(3) to 20 decimal places is:

1.20205690315959428539

Reaching any level of accuracy using our implementation of the Riemann zeta function takes quite a while:

>>> find_apery()
Found 5 digits at iteration 269 (0.329000 seconds)
Found 6 digits at iteration 744 (0.938000 seconds)
Found 7 digits at iteration 12580 (17.485000 seconds)
Found 8 digits at iteration 12580 (17.500000 seconds)
Found 9 digits at iteration 55973 (79.001000 seconds)
Found 10 digits at iteration 91597 (132.845000 seconds)
Found 11 digits at iteration 228286 (350.190000 seconds)

Here we’ve reached 5 digits of accuracy in 269 iterations (or n=269), and 10 digits of accuracy after 91,597 iterations. To work efficiently with the Riemann zeta function, we can specify the number of digits in our resulting number, and the number of iterations to spend reaching for accuracy.

Here’s the Python code used for calculating Apery’s constant, with the accompanying methods for creating a Riemann function generator, or just retrieving a Riemann function value after a set number of iterations:

#!/usr/bin/python

from decimal import *
import time

def zeta_gen(n, digits, iter):
	getcontext().prec = digits
	curSum = Decimal(0)
	curN = Decimal(1)
	while curN <= iter:
		curSum = curSum + (1 / (curN ** n))
		yield (curN, curSum)
		curN = curN + 1

def zeta(n, digits, iter):
	getcontext().prec = digits
	curSum = Decimal(0)
	curN = Decimal(1)
	while curN <= iter:
		curSum = curSum + (1 / (curN ** n))
		curN = curN + 1
	return curSum

def apery_gen(digits, iter):
	agen = zeta_gen(3, digits, iter)
	for (i, a) in agen:
		yield (i, a)

def apery(digits, iter):
	return zeta(3, digits, iter)

def find_apery():
	a_const = "1.20205690315959428539"
	looking_at = 5
	start = time.time()
	for (i, a) in apery_gen(40, Decimal("1000000000")):
		while a.to_eng_string()[:looking_at + 2] == a_const[:looking_at + 2]:
			curtime = time.time()
			print "Found %d digits at iteration %d (%f seconds)" % (looking_at, i, curtime - start)
			looking_at += 1

The zeta and zeta_gen methods can be used to find the Riemann zeta function value for a number. This usage will retrieve the Riemann value ζ(3) with 100 digits after 10,000 iterations (accurate to about 6 digits):

z = zeta(3, 100, 10000)

Or you can use the generator version to watch the value converge:

for z in zeta(3, 100, 10000):
    print z

Using the code from the Calculating π, I’ve generated the first 4,000 digits of π. This is, of course, an incredibly small number of digits compared to the 1.2 trillion digits of π calculated by Yasumasa Kanada. Still, it’s a fun bit of calculation, and somewhat rewarding to see the digits flash by in the background as real work gets done.

The digits:

3.14159265 3589793238 4626433832 7950288419 7169399375 1058209749
4459230781 6406286208 9986280348 2534211706 7982148086 5132823066
4709384460 9550582231 7253594081 2848111745 0284102701 9385211055
5964462294 8954930381 9644288109 7566593344 6128475648 2337867831
6527120190 9145648566 9234603486 1045432664 8213393607 2602491412
7372458700 6606315588 1748815209 2096282925 4091715364 3678925903
6001133053 0548820466 5213841469 5194151160 9433057270 3657595919
5309218611 7381932611 7931051185 4807446237 9962749567 3518857527
2489122793 8183011949 1298336733 6244065664 3086021394 9463952247
3719070217 9860943702 7705392171 7629317675 2384674818 4676694051
3200056812 7145263560 8277857713 4275778960 9173637178 7214684409
0122495343 0146549585 3710507922 7968925892 3542019956 1121290219
6086403441 8159813629 7747713099 6051870721 1349999998 3729780499
5105973173 2816096318 5950244594 5534690830 2642522308 2533446850
3526193118 8171010003 1378387528 8658753320 8381420617 1776691473
0359825349 0428755468 7311595628 6388235378 7593751957 7818577805
3217122680 6613001927 8766111959 0921642019 8938095257 2010654858
6327886593 6153381827 9682303019 5203530185 2968995773 6225994138
9124972177 5283479131 5155748572 4245415069 5950829533 1168617278
5588907509 8381754637 4649393192 5506040092 7701671139 0098488240
1285836160 3563707660 1047101819 4295559619 8946767837 4494482553
7977472684 7104047534 6462080466 8425906949 1293313677 0289891521
0475216205 6966024058 0381501935 1125338243 0035587640 2474964732
6391419927 2604269922 7967823547 8163600934 1721641219 9245863150
3028618297 4555706749 8385054945 8858692699 5690927210 7975093029
5532116534 4987202755 9602364806 6549911988 1834797753 5663698074
2654252786 2551818417 5746728909 7777279380 0081647060 0161452491
9217321721 4772350141 4419735685 4816136115 7352552133 4757418494
6843852332 3907394143 3345477624 1686251898 3569485562 0992192221
8427255025 4256887671 7904946016 5346680498 8627232791 7860857843
8382796797 6681454100 9538837863 6095068006 4225125205 1173929848
9608412848 8626945604 2419652850 2221066118 6306744278 6220391949
4504712371 3786960956 3643719172 8746776465 7573962413 8908658326
4599581339 0478027590 0994657640 7895126946 8398352595 7098258226
2052248940 7726719478 2684826014 7699090264 0136394437 4553050682
0349625245 1749399651 4314298091 9065925093 7221696461 5157098583
8741059788 5959772975 4989301617 5392846813 8268683868 9427741559
9185592524 5953959431 0499725246 8084598727 3644695848 6538367362
2262609912 4608051243 8843904512 4413654976 2780797715 6914359977
0012961608 9441694868 5558484063 5342207222 5828488648 1584560285
0601684273 9452267467 6788952521 3852254995 4666727823 9864565961
1635488623 0577456498 0355936345 6817432411 2515076069 4794510965
9609402522 8879710893 1456691368 6722874894 0560101503 3086179286
8092087476 0917824938 5890097149 0967598526 1365549781 8931297848
2168299894 8722658804 8575640142 7047755513 2379641451 5237462343
6454285844 4795265867 8210511413 5473573952 3113427166 1021359695
3623144295 2484937187 1101457654 0359027993 4403742007 3105785390
6219838744 7808478489 6833214457 1386875194 3506430218 4531910484
8100537061 4680674919 2781911979 3995206141 9663428754 4406437451
2371819217 9998391015 9195618146 7514269123 9748940907 1864942319
6156794520 8095146550 2252316038 8193014209 3762137855 9566389377
8708303906 9792077346 7221825625 9966150142 1503068038 4477345492
0260541466 5925201497 4428507325 1866600213 2434088190 7104863317
3464965145 3905796268 5610055081 0665879699 8163574736 3840525714
5910289706 4140110971 2062804390 3975951567 7157700420 3378699360
0723055876 3176359421 8731251471 2053292819 1826186125 8673215791
9841484882 9164470609 5752706957 2209175671 1672291098 1690915280
1735067127 4858322287 1835209353 9657251210 8357915136 9882091444
2100675103 3467110314 1267111369 9086585163 9831501970 1651511685
1714376576 1835155650 8849099898 5998238734 5528331635 5076479185
3589322618 5489632132 9330898570 6420467525 9070915481 4165498594
6163718027 0981994309 9244889575 7128289059 2323326097 2997120844
3357326548 9382391193 2597463667 3058360414 2813883032 0382490375
8985243744 1702913276 5618093773 4440307074 6921120191 3020330380
1976211011 0044929321 5160842444 8596376698 3895228684 7831235526
5821314495 7685726243 3441893039 6864262434 1077322697 8028073189
1544110104 4682325271 6201052652 2721116590 6

For reference, this calculation took a little over 3.25 hours (11802 seconds) to complete. If 4,000 digits aren’t enough for you there are a few more resources you could check.

Continuing the kata from earlier in the week, where I calculated digits of e, now I’m working through digits of π. Another irrational number calculation? But why?

I find calculation like these interesting to work through, because they typically provide an extreme case for playing with optimization. In the case of calculating digits of π, I’m using a Liebniz calculation simplified for arctangent, which takes the form:

π / 2 = summation [n = 0 -> ∞) of n! / (2n + 1)!!

This equation is a prime candidate for a naive and optimized solution because of the factorial and double factorial. In a naive implementation every step of the summation, which will be considerable in number, would require a recalculation of two factorials.

In the code here, this is bypassed by keeping a running factorial calculation going throughout the calculation of the summation, resulting in just two multiplications at every step, instead of an entire factorial. This is a fairly simple optimization, and the algorithm itself is slow in comparison with the slew of modern algorithms for calculating π, but it provides an interesting kata.

The code below is used for the calculations. As with the Euler constant calculation code, this can be used in two ways. You can watch as the π calculation converges with:

for pi in pi_gen(1000):
    print pi

Or you can wait for the calculation (which is much quicker), by using the pi, or calculate_pi methods:

print pi(1000)
calculate_pi(1000)

The calculate_pi method only adds some timing information at the end of the calculation, so the two are interchangeable. The code:

from decimal import *
import time

def factorial(n):
 	if n <= 1:
 		return Decimal(1)
 	t = n
 	r = Decimal(1)
 	while t > 1:
 		r = r * t
 		t = t - 1
	return r

def doubleFactorial(n):
	if n <= 1:
		return Decimal(1)
	t = n
	r = Decimal(1)
	while t > 1:
		r = r * t
		t = t - 2
	return r

def pi_gen(digits):
	getcontext().prec = digits
	lastSum = Decimal(-1)
	curSum = Decimal(1)
	curN = Decimal(1)
	lastFacN = Decimal(1)
	lastDFacN = Decimal(1)
	while lastSum != curSum:
		lastSum = curSum
		lastFacN = lastFacN * curN
		lastDFacN = lastDFacN * ( (curN * 2) + 1)
		curSum = curSum + (lastFacN / lastDFacN)
		curPi = curSum * 2
		yield curPi
		curN = curN + 1

def pi_with_fac(digits):
	getcontext().prec = digits
	lastSum = Decimal(-1)
	curSum = Decimal(1)
	curN = Decimal(1)
	while lastSum != curSum:
		lastSum = curSum
		curSum = curSum + (factorial(curN) / doubleFactorial( (2 * curN) + 1))
		curN = curN + 1
	return curSum * 2

def pi(digits):
	getcontext().prec = digits
	lastSum = Decimal(-1)
	curSum = Decimal(1)
	curN = Decimal(1)
	lastFacN = Decimal(1)
	lastDFacN = Decimal(1)
	while lastSum != curSum:
		lastSum = curSum
		lastFacN = lastFacN * curN
		lastDFacN = lastDFacN * ( (curN * 2) + 1)
		curSum = curSum + (lastFacN / lastDFacN)
		curN = curN + 1
	return curSum * 2

def wrapNumber(s, chunkSize=10, lineSize=60):
	s = s.replace("n", "")
	st = ""
	n = 0
	for c in s:
		st = st + c
		n = n + 1
		if (n % chunkSize) == 0:
			st = st + " "
		if (n % lineSize) == 0:
			st = st + "n"
	return st

def calculate_pi_with_fac(n):
    start = time.time()
    p = pi_with_fac(n)
    end = time.time()
    print wrapNumber(p.to_eng_string())
    print
    print "(pi) to %d decimal places" % (n)
    print "Calculated in %f seconds" % (end - start)

def calculate_pi(n):
    start = time.time()
    p = pi(n)
    end = time.time()
    print wrapNumber(p.to_eng_string())
    print
    print "(pi) to %d decimal places" % (n)
    print "Calculated in %f seconds" % (end - start)

Along with the pi method is the naive implementation using an inline factorial, pi_with_fac, which can be run for comparison. Even with a relatively short representation of π, the results are interesting:

>>> calculate_pi(100)
3.14159265 3589793238 4626433832 7950288419 7169399375 1058209749
4459230781 6406286208 9986280348 2534211706 0

(pi) to 100 decimal places
Calculated in 1.094000 seconds

>>> calculate_pi_with_fac(100)
3.14159265 3589793238 4626433832 7950288419 7169399375 1058209749
4459230781 6406286208 9986280348 2534211706 0

(pi) to 100 decimal places
Calculated in 57.626000 seconds

So our optimized method to a little more than a second, compared to the naive implementation, which took closer to a minute.

Using the code from Calculating e I’ve generated the first 10,000 digits of Eulers constant:

2.71828182 8459045235 3602874713 5266249775 7247093699 9595749669
6762772407 6630353547 5945713821 7852516642 7427466391 9320030599
2181741359 6629043572 9003342952 6059563073 8132328627 9434907632
3382988075 3195251019 0115738341 8793070215 4089149934 8841675092
4476146066 8082264800 1684774118 5374234544 2437107539 0777449920
6955170276 1838606261 3313845830 0075204493 3826560297 6067371132
0070932870 9127443747 0472306969 7720931014 1692836819 0255151086
5746377211 1252389784 4250569536 9677078544 9969967946 8644549059
8793163688 9230098793 1277361782 1542499922 9576351482 2082698951
9366803318 2528869398 4964651058 2093923982 9488793320 3625094431
1730123819 7068416140 3970198376 7932068328 2376464804 2953118023
2878250981 9455815301 7567173613 3206981125 0996181881 5930416903
5159888851 9345807273 8667385894 2287922849 9892086805 8257492796
1048419844 4363463244 9684875602 3362482704 1978623209 0021609902
3530436994 1849146314 0934317381 4364054625 3152096183 6908887070
1676839642 4378140592 7145635490 6130310720 8510383750 5101157477
0417189861 0687396965 5212671546 8895703503 5402123407 8498193343
2106817012 1005627880 2351930332 2474501585 3904730419 9577770935
0366041699 7329725088 6876966403 5557071622 6844716256 0798826517
8713419512 4665201030 5921236677 1943252786 7539855894 4896970964
0975459185 6956380236 3701621120 4774272283 6489613422 5164450781
8244235294 8636372141 7402388934 4124796357 4370263755 2944483379
9801612549 2278509257 7825620926 2264832627 7933386566 4816277251
6401910590 0491644998 2893150566 0472580277 8631864155 1956532442
5869829469 5930801915 2987211725 5634754639 6447910145 9040905862
9849679128 7406870504 8958586717 4798546677 5757320568 1288459205
4133405392 2000113786 3009455606 8816674001 6984205580 4033637953
7645203040 2432256613 5278369511 7788386387 4439662532 2498506549
9588623428 1899707733 2761717839 2803494650 1434558897 0719425863
9877275471 0962953741 5211151368 3506275260 2326484728 7039207643
1005958411 6612054529 7030236472 5492966693 8115137322 7536450988
8903136020 5724817658 5118063036 4428123149 6550704751 0254465011
7272115551 9486685080 0368532281 8315219600 3735625279 4495158284
1882947876 1085263981 3955990067 3764829224 4375287184 6245780361
9298197139 9147564488 2626039033 8144182326 2515097482 7987779964
3730899703 8886778227 1383605772 9788241256 1190717663 9465070633
0452795466 1855096666 1856647097 1134447401 6070462621 5680717481
8778443714 3698821855 9670959102 5968620023 5371858874 8569652200
0503117343 9207321139 0803293634 4797273559 5527734907 1783793421
6370120500 5451326383 5440001863 2399149070 5479778056 6978533580
4896690629 5119432473 0995876552 3681285904 1383241160 7226029983
3053537087 6138939639 1779574540 1613722361 8789365260 5381558415
8718692553 8606164779 8340254351 2843961294 6035291332 5942794904
3372990857 3158029095 8631382683 2914771163 9633709240 0316894586
3606064584 5925126994 6557248391 8656420975 2685082307 5442545993
7691704197 7780085362 7309417101 6343490769 6423722294 3523661255
7250881477 9223151974 7780605696 7253801718 0776360346 2459278778
4658506560 5078084421 1529697521 8908740196 6090665180 3516501792
5046195013 6658543663 2712549639 9085491442 0001457476 0819302212
0660243300 9641270489 4390397177 1951806990 8699860663 6583232278
7093765022 6014929101 1517177635 9446020232 4930028040 1867723910
2880978666 0565118326 0043688508 8171572386 6984224220 1024950551
8816948032 2100251542 6494639812 8736776589 2768816359 8312477886
5201411741 1091360116 4995076629 0779436460 0585194199 8560162647
9076153210 3872755712 6992518275 6879893027 6176114616 2549356495
9037980458 3818232336 8612016243 7365698467 0378585330 5275833337
9399075216 6069238053 3698879565 1372855938 8349989470 7416181550
1253970646 4817194670 8348197214 4888987906 7650379590 3669672494
9925452790 3372963616 2658976039 4985767413 9735944102 3744329709
3554779826 2961459144 2936451428 6171585873 3974679189 7571211956
1873857836 4475844842 3555581050 0256114923 9151889309 9463428413
9360803830 9166281881 1503715284 9670597416 2562823609 2168075150
1777253874 0256425347 0879089137 2917228286 1151591568 3725241630
7722544063 3787593105 9826760944 2032619242 8531701878 1772960235
4130606721 3604600038 9661093647 0951414171 8577701418 0606443636
8154644400 5331608778 3143174440 8119494229 7559931401 1888683314
8328027065 5383300469 3290115744 1475631399 9722170380 4617092894
5790962716 6226074071 8749975359 2127560844 1473782330 3270330168
2371936480 0217328573 4935947564 3341299430 2485023573 2214597843
2826414216 8487872167 3367010615 0942434569 8440187331 2810107945
1272237378 8612605816 5668053714 3961278887 3252737389 0392890506
8653241380 6279602593 0387727697 7837928684 0932536588 0733988457
2187460210 0531148335 1323850047 8271693762 1800490479 5597959290
5916554705 0577751430 8175112698 9851884087 1856402603 5305583737
8324229241 8562564425 5022672155 9802740126 1797192804 7139600689
1638286652 7700975276 7069777036 4392602243 7284184088 3251848770
4726384403 7953016690 5465937461 6193238403 6389313136 4327137688
8410268112 1989127522 3056256756 2547017250 8634976536 7288605966
7527408686 2740791285 6576996313 7897530346 6061666980 4218267724
5605306607 7389962421 8340859882 0718646826 2321508028 8286359746
8396543588 5668550377 3131296587 9758105012 1491620765 6769950659
7153447634 7032085321 5603674828 6083786568 0307306265 7633469774
2956346437 1670939719 3060876963 4953288468 3361303882 9431040800
2968738691 1706666614 6800015121 1434422560 2387447432 5250769387
0777751932 9994213727 7211258843 6087158348 3562696166 1980572526
6122067975 4062106208 0649882918 4543953015 2998209250 3005498257
0433905535 7016865312 0526495614 8572492573 8620691740 3695213533
7325316663 4546658859 7286659451 1364413703 3139367211 8569553952
1084584072 4432383558 6063106806 9649248512 3263269951 4603596037
2972531983 6842336390 4632136710 1161928217 1115028280 1604488058
8023820319 8149309636 9596735832 7420249882 4568494127 3860566491
3525267060 4623445054 9227581151 7093149218 7959271800 1940968866
9868370373 0220047531 4338181092 7080300172 0593553052 0700706072
2339994639 9057131158 7099635777 3590271962 8506114651 4837526209
5653467132 9002599439 7663114545 9026858989 7911583709 3419370441
1551219201 1716488056 6945938131 1838437656 2062784631 0490346293
9500294583 4116482411 4969758326 0118007316 9943739350 6966295712
4102732391 3874175492 3071862454 5432220395 5273529524 0245903805
7445028922 4688628533 6542213815 7221311632 8811205214 6489805180
0920247193 9171055539 0113943316 6815158288 4368760696 1102505171
0073927623 8555338627 2553538830 9606716446 6237092264 6809671254
0618695021 4317621166 8140097595 2814939072 2260111268 1153108387
3176173232 3526360583 8173151034 5957365382 2353499293 5822836851
0078108846 3434998351 8404451704 2701893819 9424341009 0575376257
7675711180 9008816418 3319201962 6234162881 6652137471 7325477727
7834887743 6651882875 2156685719 5063719365 6539038944 9366421764
0031215278 7022236646 3635755503 5655769488 8654950027 0853923617
1055021311 4741374410 6134445544 1921013361 7299628569 4899193369
1847294785 8072915608 8510396781 9594298331 8648075608 3679551496
6364489655 9294818785 1784038773 3262470519 4505041984 7742014183
9477312028 1588684570 7290544057 5106012852 5805659470 3046836344
5926525521 3700806875 2009593453 6073162261 1872817392 8074623094
6853678231 0609792159 9360019946 2379934342 1068781349 7346959246
4697525062 4695861690 9178573976 5951993929 9399556754 2714654910
4568607020 9901260681 8704984178 0791739240 7194599632 3060254707
9017745275 1318680998 2284730860 7665368668 5551646770 2911336827
5631072233 4672611370 5490795365 8345386371 9623585631 2618387156
7741187385 2772292259 4743373785 6955384562 4680101390 5727871016
5129666367 6445187246 5653730402 4436841408 1448873295 7847348490
0030194778 8802046032 4660842875 3518483649 5919508288 8323206522
1281041904 4804724794 9291342284 9519700226 0131043006 2410717971
5027934332 6340799596 0531446053 2304885289 7291765987 6016667811
9379323724 5385720960 7582277178 4833616135 8261289622 6118129455
9274627671 3779448758 6753657544 8614076119 3112595851 2655759734
5730153336 4263076798 5443385761 7153334623 2527057200 5303988289
4990342595 6623297578 2488735029 2591668258 9445689465 5992658454
7626945287 8051650172 0674785417 8879822768 0653665064 1910973434
5288783386 2172615626 9582654478 2056729877 5642632532 1594294418
0399432170 0009054265 0763095588 4658951717 0914760743 7136893319
4690909819 0450129030 7099566226 6203031826 4936573369 8419555776
9637876249 1885286568 6607600566 0256054457 1133728684 0205574416
0308370523 1224258722 3438854123 1794813885 5007568938 1124935386
3186352870 8379984569 2619981794 5233640874 2959118074 7453419551
4203517261 8420084550 9170845682 3682008977 3945584267 9214273477
5608796442 7920270831 2150156406 3413416171 6644806981 5483764491
5739001212 1704154787 2591998943 8253649505 1477137939 9147205219
5290793961 3762110723 8494290616 3576045962 3125350606 8537651423
1153496656 8371511660 4220796394 4666211632 5515772907 0978473156
2782775987 8813649195 1257483328 7937715714 5909106484 1642678309
9497236744 2017586226 9402159407 9244805412 5536043131 7992696739
1575424192 9660731239 3763542139 2306178767 5395871143 6104089409
9660894714 1834069836 2993675362 6215452472 9846421375 2891079884
3813060955 5262272083 7518629837 0667872244 3019579379 3786072107
2542772890 7173285487 4374355781 9665117166 1833088112 9120245204
0486822000 7234403502 5448202834 2541878846 5360259150 6445271657
7000445210 9773558589 7622655484 9416217149 8953238342 1600114062
9507184904 2778925855 2743035221 3968356790 1807640604 2138307308
7744601708 4268827226 1177180842 6643336517 8000217190 3449234264
2662922614 5600433738 3868335555 3434530042 6481847398 9215627086
0956506293 4040526494 3244261445 6659212912 2564889356 9655009154
3064261342 5266847259 4914314239 3988454324 8632746184 2846655985
3323122104 6625989014 1712103446 0842716166 1900125719 5870793217
5696985440 1339762209 6749454185 4071184464 3394699016 2698351607
8489245140 5894094639 5267807354 5797003070 5116368251 9487701189
7640028276 4841416058 7206184185 2971891540 1968825328 9309149665
3457535714 2731848201 6384644832 4990378860 6900807270 9327673127
5819665639 4114896171 6832980455 1397295066 8760474091 5420428429
9935410258 2911350224 1690769431 6685742425 2250902693 9034814856
4513030699 2519959043 6384028429 2674125734 2244776558 4177886171
7372654620 8549829449 8946787350 9295816526 3207225899 2368768457
0178230380 9656788311 2289305809 1405726108 6588484587 3101658151
1675333276 7488701482 9167419701 5125597825 7270740643 1808601428
1490241467 8047232759 7684269633 9357735429 3018673943 9716388611
7642090040 6866339885 6841681003 8723892144 8317607011 6684503887
2123643670 4331409115 5733280182 9779887365 9091665961 2402021778
5588548761 7616198937 0794380056 6633648843 6508914480 5571039765
2146960276 6258359905 1987042300 1794655366 2

Of course, many, many more digits have been calculated using smarter and more efficient means. This does raise the question, “Why calculate 10,000,000 decimal places of a constant?” How many decimal places are actually useful?

If we consider that the diameter of the known universe is approximately 10²⁶ meters, or 10³⁹ nanometers, than if we attempted to calculate the area of the universe, the fist 40 or so digits of pi would be significant in the calculation. If the same level of accuracy is true for e, then of the 100,000,000,000 digits known, only 4 * 10^-8 percent of the digits are really useful in any theoretical sense.