Archive for the ‘Math’ Category

Fractal Memory Usage and a Big Number

Thursday, June 12th, 2008

In a previous post I said I’d talk about 4**(4**(4**4)).
First, about the number. I first saw it mentioned in a math lesson back in high-school, when the teacher wanted to demonstrate estimation abilities.
He wrote it on the board, like so:

4444

and then asked students to estimate how big that number is. Turns out this number is much bigger than most students thought. (For comparison, the number of atoms in the observable universe is 1080.)

Given that this number is so big, it should come as no surprise that computing it will consume all available memory.
What is interesting is the way the memory is consumed:

(It died with a MemoryError just before the next “jump”.)

When I first saw this I was fascinated. To those interested it is generated by the long pow function in Python, implemented in longobject.c. According to the comment there, the algorithm is a left to right 5-ary exponentiation, taken from “Handbook of Applied Cryptography”.
In short, the simpler algorithm “left-to-right binary exponentiation” (used for smaller exponents) repeatedly squares the accumulator, and according to the binary digits of the exponents, multiplies it by the base. 5-ary exponentiation basically does the same thing, but in each iteration it ‘eats’ five digits of the exponent instead of just one.

It might be very interesting to study algorithms and code according to their memory usage. Once at work I used a similar memory usage graph to optimize the memory usage of some algorithm that was eating too much memory. I don’t remember the details as it happened a few years ago, but I do remember that I recognized the “type of the memory eating”, which helped me to narrow the problem down to a specific area of code, where indeed the problem was located. (I consider this a case of programming-fu :)

Notes:
1) The reason I’m talking about Python is because it has arbitrarily sized longs built-in.
2) This graph was generated using Windows’ task manager. On my Ubuntu machine it looks quite different, almost a straight line. Any guesses why?

PyKoan – The Logic Game

Sunday, June 1st, 2008

As you can probably tell, I'm back from my undeclared hiatus. I've got lots of stuff to talk about, and I'll be starting with PyKoan, one small project I've been working on lately in my spare time.

A few weeks ago I stumbled upon an article in wikipedia, regarding a logic game. This game fascinated me, especially because of the "Godel Escher Bach" connection. Quoting from Wikipedia:

Zendo is a game of inductive logic designed by Kory Heath in which one player (the "Master") creates a rule for structures ("koans") to follow, and the other players (the "Students") try to discover it by building and studying various koans which follow or break the rule. The first student to correctly state the rule wins.

As it happens I'm also taking a mathematical logic course this semester. Having read about the game, I wanted to write a similar computer game. Hence - PyKoan.

PyKoan is a game where the goal is to discover some logical rule, for example, "For each x holds x%2 == 0". This rule is applied to a koan - a list of integers. An example koan that "has Buddha nature" (follows the rule) is [0,2,8]. One which doesn't is [1].

To implement this idea, I wrote an implementation of an expression tree very similar to the one used in vial, but with a different implemented language. I also did some experimentation with the design. Since I've been talking a lot about expression trees without giving a full explanation, in a future post I'll write about the implementation used in PyKoan.

So far I didn't code a lot of the game, just the expression tree framework, and a simple rule builder. When using Python's interactive prompt, one can get a taste of how the game might feel:

In [19]: r = rulegen.create_rule(rulegen.easy_grammer, rulegen.easy_grammer_start)

In [20]: r.eval([])
Out[20]: False

In [21]: r.eval([0])
Out[21]: False

In [22]: r.eval([1])
Out[22]: False

In [23]: r.eval([2])
Out[23]: True

In [24]: r.eval([2,2])
Out[24]: True

In [25]: r.eval([2,2,3])
Out[25]: True

In [26]: r.eval([3])
Out[26]: False

In [27]: r.eval([4])
Out[27]: False

In [28]: print r
x exists such that x == 2

Here is how I would generate such an expression manually:

In [3]: from exptree import exps

In [4]: exps.Exists('x',exps.Eq('x',2))
Out[4]: exps.Exists(exps.Var(x), exps.Eq(exps.Var(x), exps.Imm(2)))

In [5]: print _
x exists such that x == 2

The game has many interesting possibilities for research, for example, computer players. Other possibilities include "just" guessing koans (not the rule itself), creating interesting and playable rules, and so on. There's a lot to do.

This time, instead of just publishing the (unfinished) code, I decided to do something different. I've opened a space in assembla, with public read access. I'm opening this project for participation: if you want to join then leave a comment, or send me an email.

(Since Assembla seems to be going through some connectivity issues right now, here's a link to a copy).

Two Mathematical Bugs

Saturday, March 8th, 2008

A few days ago, I discovered I had at least one, maybe two mathematical bugs. The first bug is in the line clipping algorithm. Here's a quick reminder of what's going on there:

There are six points of interest: the four intersection points, and the two endpoints of the line. We sort the six points according to their parameter, and take the two innermost points as the endpoints of the clipped line. To check if the line is outside the box we make sure that the result points' parameters are between 0 and 1 (which correspond to the original endpoints).

Turns out this is not always true - here's a counter example:

Counter Example for the Line Clipping Algorithm

In the counter example you can see a line, whose endpoints lie between the intersections. The algorithm described above will say this line lies completely inside the box. Clearly this is not true. How do we fix this? Simple actually. We check for each intersection point whether or not it is on the edges of the clipping rectangle. If at most one point is on - our line is outside. I hope this is the end of it. It will require more rigorous testing.

The second bug was in my solution to the random selection challenge. I wasn't too sure of my solution, but I really liked it, as its code is rather elegant. So when meeting someone I used to work with, I told her about it. After thinking about my solution a bit, she said it is probably wrong, and gave me the reason: While the probability of a point x passing the first stage is indeed (proportional to) p(x), the probability of it passing the second stage is not 1 (or constant). I still need to think about this one, so no easy fix yet. I guess this reopens the challenge, so if anyone feels like writing a good proof or a counter example for my solution, I'd be happy to see it. I'll be even happier to see a completely different solution.

Solution for the Random Selection Challenge

Saturday, February 23rd, 2008

A few days ago, I wrote up two small Python Challenges. Several people have presented solutions for the first challenge, and I also posted my solution in the comments there.

However, the second challenge remained unsolved, and I will present a solution for it in this post.

(more...)

Fun with Matrices

Wednesday, February 20th, 2008

I'll let the code speak for itself:

In [81]: m = Matrix(array([[1.0,1.0],[0.0,1.0]]))

In [82]: def my_sqrt(x, num_iters):
   ....:     r = 0.5*x
   ....:     for i in xrange(num_iters):
   ....:             r = 0.5*(r+x/r)
   ....:     return r
   ....:

In [83]: m*m
Out[83]:
array([[ 1.,  2.],
       [ 0.,  1.]])

In [84]: my_sqrt(m*m, 10)
Out[84]:
array([[ 1.,  1.],
       [ 0.,  1.]])

In [85]: m
Out[85]:
array([[ 1.,  1.],
       [ 0.,  1.]])

It's always fun to see the math work out. At first when I learned that e^A for a matrix A may also be defined using the Taylor series in the usual way, I was really amazed. It still amazes me that this stuff works out so well. This is another kind of beauty.

Small Python Challenge No. 3 – Random Selection

Sunday, February 10th, 2008

This time I'll give two related problems, both not too hard.

Lets warm up with the first:

You have a mapping between items and probabilities. You need to choose each item with its probability.

For example, consider the items ['good', 'bad', 'ugly'], with probabilities of [0.5, 0.3, 0.2] accordingly. Your solution should choose good with probability 50%, bad with 30% and ugly with 20%.

I came to this challenge because just today I had to solve it, and it seems like a common problem. Hence, it makes sense to ask 'what is the best way?'.

The second problem is slightly harder:

Assume a bell shaped function p(x) that you can 'solve'. This means that given a value y, you can get all x such that p(x)=y. For example, sin(x)^2 in [0,pi] is such a function. Given a function such as Python's random.random() that yields a uniform distribution of values in [0,1), write a function that yields a distribution proportional to p(x) in the appropriate interval.

For example, consider the function p(x) = e^(-x^2) in [-1,1]. Since p(0) = 1, and p(0.5)~0.779, the value 0 should be p(0)/p(0.5)~1.28 times more common than 0.5.

As usual, the preferred solutions are the elegant ones. Go!

note: please post your solutions in the comments, using [ python]...[ /python] tags (but without the spaces in the tags).

Interesting links – 4

Monday, February 4th, 2008

http://ivory.idyll.org/articles/advanced-swc/

"Intermediate and Advanced Software Carpentry in Python" is an excellent reading by Titus Brown. If you feel you're good with Python but want to improve it, or if you are an experienced programmer that wants to get better, this is a good place to go. I liked it.

http://c2.com/cgi/wiki?AlternateHardAndSoftLayers

During the time I was working on my Compilation course, I was thinking about the challenge of writing yacc in yacc. Well, I went searching for "yacc in yacc" and stumbled across this page. It is a part of a very strange 'pattern wiki'. It has fascinating discussions of the subjects in it, but I don't like the old-style wiki navigation. Still worth a taste.

http://www.scottaaronson.com/writings/bignumbers.html

Do you know the kind of duel where two people need to name the biggest number? If you thought something along the lines of 'hey, I'll write something with a lot of nines', well, you are in for a surprise. This article has a nice solution for the problem. Excellent read, especially if you are into computation theory.  I really liked that one.