Algorithms Geometry Math Programming Python

Two Mathematical Bugs

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.

computer science Math Programming Python Statistics Utility Functions

Solution for the Random Selection Challenge

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.

Math Programming Python

Fun with Matrices

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
array([[ 1.,  2.],
       [ 0.,  1.]])
In [84]: my_sqrt(m*m, 10)
array([[ 1.,  1.],
       [ 0.,  1.]])
In [85]: m
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.

Challenges Math Programming Python

Small Python Challenge No. 3 – Random Selection

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).

Compilation Math Programming Python

Interesting links – 4

“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.

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.

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.

Fractals Graphics Programming Python

A second look at the dragon fractal

At first when I drew the (twin)dragon fractal, I had a small bug. I used the base 1+i instead of 1-i. This also generated a very similar looking beast. Thinking about that for a while, made me curious. Just like the Mandelbrot set, maybe other interesting fractals could be generated using exactly the same method, with a different complex ‘seed’.

So I patched up my code, and made it output a series of images for complex numbers on a path. I thought the reasonable choice would be a spiral, so using t*(cost+isint) as a base I wrote a spiral that would go around the origin several times and finally land on 1-i.
It might seem obvious to you to try and make this interactive instead – why not move the mouse around and watch the different fractals that emerge? Well, I wanted a video.

I also wanted the fractal to convey a little more information. So each point in the set was given a color according to its generation. I decided after some trial and error that white was better for the older generation, and red for younger generations.
After some PIL and ffmpeg work, it was ready.

While working on the video, I witnessed some very interesting fractals with different seeds.
I was very curious as to why certain shapes emerged. For example, a square pattern was very common:
Squares fractal
This turned out to be not that surprising. Its generator number is of the shape 0+ia. So it made sense. I still didn’t figure out how come hexagon shapes were so common:
Hexagons1 Hexagons2 Hexagons1

Maybe it had something to do with pi/6, I’m not too sure about that. If it did, I would expect to see many more regular polygons.
Here is another curious one:

Another interesting phenomenon was what I started to call ‘the dragon’s heart’. You see, in the final dragon, the starting generations were spread about pretty evenly. However, with other bases, even ones generating something which is pretty similar to the dragon, the starting generations are clustered together – sometimes at the side, sometimes at the middle.

I’ve got a feeling (not proven or anything) that to be space filling, a fractal’s starting generations should be spread about. What do you think?

Fractals Math Programming Python

Fractals in 10 minutes no. 3 – The Dragon

When first I looked through the pages of the book “Hacker’s Delight”, I found myself looking at the chapter about bases. There I learned a very curious fact – with the digits of 0,1 and the base of -2, you can represent any integer. Right afterwards I learned something even more interesting – with the digits of 0,1 and the base of 1-i, you can represent and number of the form a+bi where a and b are integers. Having nothing to do with this curious fact, I let the subject go.
Some time later, I was reading through Knuth’s “Art of Computer Programming”, and found that with the base of (1-i)^-1, and digits of 0,1 you can generate the dragon fractal!

The dragon fractal

Generating the fractal is quite simple actually:

def create_dragon_set(n):
    """calculate the dragon set, according to Knuth"""
    s = set([0.0+0.0j])
    for i in range(n):
        new_power = (1.0-1.0j)**(-i)
        s |= set(x+new_power for x in s)
    return s

(By the way, can you do it better?)

The annoying part is converting the complex numbers to drawable integer points. After doing so, I used PIL to draw the jpeg.
Here’s a link to the code.

Cryptography Math Programming Python

Using zpint in Cryptography Homework

Originally, I wrote zpint to help me with my algebric structures homework. For those not familiar with it, it allows you to do computations modulo p the same way you do for ints. Not surprisingly, I found myself using it in my cryptography homework as well. At first I used it for trying to breaking the Hill Cipher. Then for some computations such as using the Chinese Remainder Theorem to solve equations, or to compute RSA values. Last it was for checking if algorithms I wrote were correct.

Usually I prefer to write a script than do all the computations by hand, even if writing the script is more work. Mostly because it is work I enjoy more. This time I had a script ready. I improved it a little bit during the semester, and it is still available. It is certainly not the fastest way to compute, but it is a fast way of doing computations.

Algorithms Fractals Programming Python

Optimizing the recursive spring

Remember the recursive spring? Well, when I worked on it, it was quite slow. I mean very slow. Especially at the higher levels.

The original reason was the recursion. Consider this – to ‘springify’ a track, you had to sample it twice, so that you could get the heading. Here is an illustration:

(lower track is ‘inner’ track)

Now it is clear – the number of recursive calls made on all levels grow exponentially. There is some hope. The recursive calls use points calculated previously. It seems that with good caching we’ll be able to cancel out the exponential growth, and make it linear.

Actually, good caching did much more than that – I cached the last track as well. Which meant that the second time around drawing the track didn’t cost anything – just some dictionary lookups. This yielded a pretty elegant result.

To those curious, here is the code. Before I started working on it I hacked together a small 3d explorer with PyOpenGL, so that I’ll be able to play with the 3d models. This module is useful by itself as well, and someday soon I’ll dedicate a post to it. (Only after making it better. Right now I don’t consider it high-quality code.)

To run it, use, and play with it a little. If you look at the code, there are many parameters to play with, yielding interesting results.

Fractals Graphics Programming Python

Fractals in 10 minutes No. 2 – Recursive Spring

Imagine a straight wire. Bend this wire until its ends meet. You get a ring. Next stage.
Take another straight wire, bend it as before, but this time, don’t let the ends meet, instead continue bending the wire more and more, until you get a spring. Now, Think of this spring as the first wire, and bend it until the spring’s ends meet. We get a spring-ring. (See pictures below)

At this point it should be obvious – instead of making the ends of the spring meet, we can continue bending the spring into a meta-spring, and again make this meta-spring’s ends meet. Rinse and repeat ad-infinitum.

At first I had a hard time picturing it, and I thought it would be tough to write about it without drawing it. So I set with pen and paper, and tried to draw it. That was even harder than picturing it. So I thought I could write a script to draw it, and that turned out to be pretty easy.

Here are the first four stages. The second and third are shown from the side. Click for bigger versions.

To implement this, I wrote a function that draws any track. A track is defined to be a function takes a float argument between 0 and 1, and returns the appropriate 3d point in between, where 0 is one end of the track and 1 is the other. (Actually, I also decided that tracks should return a “right hand” vector as well). Then I wrote a ring:

def ring(t):
    pi = math.pi
    cos = math.cos
    sin = math.sin
    x = numpy.array((2*cos(2*pi*t),2*sin(2*pi*t),0))
    return (x,-x/3)

Lastly I wrote a wrapper function, that takes any track, and springifies it. It was called springify_track. At this point I could do this:
(The constants are the radius of the spring loops, the sampling delta, and the number loops.)

t1 = ring
t2 = springify_track(t1, 0.25, 0.1, 50)
t3 = springify_track(t2, 0.06, 0.0011, 5000)
t4 = springify_track(t3, 0.011, 0.000012, 500000)