Archive for the ‘Programming’ Category

Mandelbrot

Saturday, November 8th, 2014

mandelbrot1

Math can be beautiful. One of the mathematical things that has always intrigued me has been the visualization of the Mandelbrot set. I never tried to generate it myself until just recently. The only understanding I had previously of the Mandelbrot set was vague. I looked for some explanation on how to do the computation, but I couldn’t find an adequately succinct explanation. My brother helped me take the Mandelbrot equation and massage it into an understandable format that I could then plug into a program to generate my own image of the Mandelbrot set.

Generating the Mandelbrot set is straight forward once you understand what the Mandelbrot set really is. The fractal depicts a grid of coordinates along the x and y axis that either fall into the set or not. The formula gets computed using the x and y position for the pixel in question, and then the resulting number get put back into the equation. The resulting number either gets larger or smaller. If it passes a certain threshold, then we decide it is not a part of the set, and move on. The resulting image is a bit map of numbers that are either part of the Mandelbrot set, or not. The fun gradients and colors, with which we are all familiar, represent how many times the calculation had to be performed before the number passed the threshold.

The Mandelbrot equation is as follows, where z_(n+1) is the value for which we are solving, (z_n)² is our accumulated value so far, and C is our complex coordinate.

z_(n+1) = (z_n)² + C

We are just focusing on the right side of the equation for now, because the left contains the values that we want to solve for. Let’s expand the right side of the equation with complex coordinates so we can simplify the answer.

= (x′ + y′ × i)² + x + y × i

The next step is to expand the accumulated value. We can use our old friend FOIL.

= x′² + 2(x′ × y′ × i) + (y′ × i)² + x + y × i

The i² simplifies down to -1, reducing the number of imaginary numbers to worry about.

= x′² + 2(x′ × y′ × i) – y′² + x + y × i

Next we will move all the components with an imaginary number (i) to the right.

= (x′² + x – y′² )+( 2(x′ × y′ × i) + y × i)

Now we can the simplify by isolating the i.

= (x′² + x – y′²)+(2 × x′ × y′ + y) × i

Now the right side of the equation resembles a complex coordinate, which happens to be what we are solving for on the left.

x′′ + y′′ × i = …

So, each element of the left hand side of the equation can be solved with its counterpart element on the right.

x′′ = x′² + x – y²

y′′ = 2 × x′ × y′ + y

i = i

Since the i is equal to i, and remains independent from the rest of the variables, we can just ignore it.

Now, all we need to do is plug x and y to generate our new x′′ and y′′, and then keep plugging in the new values until the numbers go out of bounds, or we reach our threshold.

With the equation in hand, I started out by prototyping the program in Python, and outputting a Portable PixelMap (PPM) file because it is a dead-simple file format. When I first ran the program, it generated this image:

man1

I was really disappointed that I didn’t get a meaningful image the first time I ran it. I verified the equation I was using, and the basic logic of the code, but it still came out like this. I kicked up the threshold a little bit to see what would happen, and suddenly I realized what the problem was. I was only rendering the upper left quadrant of the image. I finagled the bounds of the image, and ended up with the recognizable shape here:

man4

This certainly looked better, but something was still off. In all the Mandelbrot images I had seen, the different threshold levels looked smooth and curved, not jagged like mine were looking. It turned out I had misunderstood how to check if the calculation was out-of-bounds. This code is the culprit:

if (abs(xT) + abs(yT)) > 10:

This is the conditional statement that decides whether the number has fallen out-of-bounds. The key is to check if the number has left the circular area that contains the entire Mandelbrot set, which has a radius of 2. This would mean that the enclosing circle is something like x′′² + y′′² = 2². The corrected out-of-bounds check looks like this:

if ((xT**2) + (yT**2)) > 4:

With that change, the images started coming out beautifully. With that, I started adding in command line options to set the bounds of the image, the depth, the gradient colors, etc. in order to easily render custom images of the set.

wide

Now that I did it all in Python, I decided I should do it in C. Rendering 1024 pixel square image in Python takes about 15 seconds. To compare, rendering the same image in C takes .1 seconds; substantially faster!

The following code from my function, mandel(), calculates whether any given point is in the Mandelbrot set. The function takes an x and a y coordinate, and the threshold level which I call “depth”. It returns a -1 if the number is in the set, or it will return an integer representing the number of iterations it took for the point to become out-of-bounds.

int mandel(double x, double y, int depth)
{
    double xP=0, yP=0; // Prime Values
    double xT=0, yT=0; // Temporary Values
    int i=0;
    for (i=0; i<depth; i++)
    {
        xT = (pow(xP, 2)) + x - (pow(yP, 2));
        yT = 2 * xP * yP + y;
        if (pow(fabs(xT),2) + pow(fabs(yT),2) > 4)
            return i;
        xP = xT;
        yP = yT;
    }
    return -1;
}

mandelbrot_angled_green_scaled

The current code for both the Python and C versions are available from my GitHub page.

Project Euler 12

Friday, September 5th, 2014

I was recently going through some old Project Euler solutions, and noticed that my solution to problem 12 was never optimized. The problems asks:

The sequence of triangle numbers is generated by adding the natural numbers. So the 7th triangle number would be 1 + 2 + 3 + 4 + 5 + 6 + 7 = 28. The first ten terms would be:

1, 3, 6, 10, 15, 21, 28, 36, 45, 55, …

Let us list the factors of the first seven triangle numbers:

1: 1
3: 1,3
6: 1,2,3,6
10: 1,2,5,10
15: 1,3,5,15
21: 1,3,7,21
28: 1,2,4,7,14,28

We can see that 28 is the first triangle number to have over five divisors.

What is the value of the first triangle number to have over five hundred divisors?

The way I originally solved this problem was a brute-force approach following this basic logic:

  1. Calculate next triangle number
  2. Try dividing number by every integer smaller than the number
  3. Count the number of integers that divide evenly
  4. Go to step 1 and repeat until divisors are greater than 500

This process yielded the results, but it is far from elegant or optimized. it took 14 minutes to calculate the answer! I knew that dividing by each number was a waste; there had to be a better way.

Although I didn’t know it, I felt fairly certain that there was a relationship between prime factors and the number of divisors. It turns out I was right. There is a trick to calculate the number of divisors if you know the prime factors.

Using the number 28 as an example, we can find that the prime factors are 2, 2, and 7, which can be written as 22 × 71. Interestingly, if we take the exponents, add 1 to each, and then multiply them together, we can calculate the number of divisors. For this example, this becomes (2+1) × (1+1), which is 2 × 3, which equals 6. (This is also mentioned in the Wikipedia article for divisor briefly in the section titled ‘Further notions and facts’).

Equipped with this knowledge, I was able to devise a more efficient method to solve this problem:

  1. Pre-calculate primes up to x amount (to save effort recomputing each time)
  2. Calculate next triangle number that is divisible by 3 or 5 (I found a pattern that the most divisible triangle numbers were divisible by 3 or 5)
  3. Go through primes and divide number by primes until number equals 1
  4. Convert the prime factors to divisors
  5. Go to step 1 and repeat until divisors are greater than 500

While this method is more complex, it yielded great results. Once I was done refining this method, I got the compute time down to .086 seconds. That is a very substantial improvement!

The interesting part of the code is the findDivisors function, seen below.

int findDivisors(int input, int *primes, int primesToFind)
{
    int i, count=1;
    int primeCount[primesToFind];
    for (i=0; i<primesToFind; i++)
    primeCount[i] = 0;

    // Prime factorization
    while (input != 1)
        for (i=0; i<primesToFind; i++)
            if (input % primes[i] == 0)
            {
                input = input / primes[i];
                primeCount[i]++;
            }

    // Convert prime factorization to divisors
    for (i=0; i<primesToFind; i++)
        if (primeCount[i] > 0)
            count *= (primeCount[i]+1);

    return count;
}

In this function, the inputs are:

  • The triangle number (input)
  • The array of pre-calculated prime numbers (primes)
  • The number of primes in the array (primesToFind)

As the division of primes occurs, I keep track in a separate array (primeCount) how many times a given prime index is used to divide. This keeps track of the prime factorization I talked about earlier so we can finally compute the number of divisors. Once computed, we pass back the count, and the program moves on to the next triangle number.

So there you have it. My entire solution is available on my GitHub page.

Emacs on Mac OS X

Monday, May 19th, 2014

carbon-emacs-iconAt work I use a Linux work station a lot, and one of the things I really enjoy is that when I use Emacs from the command line, it opens up the graphical user interface. This, of course, is not the behavior on OS X. Additionally, the version that ships with Mac OS X is woefully out of date (version 22.1.1 at the time of writing).

I decided to do something about it, and grabbed an up-to-date build of Emacs from emacsformacosx.com. By default, it gets installed in /Applications, which will work just fine. It provides a graphical version of Emacs with support for normal Mac OS X short cuts, like command-c, etc.

But, Michael, how does this solve the problem of the Emacs from the command line being out of date? Well, Inside the application bundle is the Emacs executable that can be run from the command line, sans GUI. All you have to do is create an alias with this code (in my case, I made a bash script that resides in my ~/.bin named ’emacs’) with the following code:

#!/bin/sh
if
    [ -z "$SSH_CONNECTION" ]
then
    /Applications/Emacs.app/Contents/MacOS/Emacs "$@"
else
    /Applications/Emacs.app/Contents/MacOS/Emacs -nw "$@"
fi

This makes it so that if the terminal detects that you are running over an SSH (non-graphical) connection, it will run the non-windowed version of Emacs, otherwise it will launch the graphical version. It works just like it does in most Linux distros, and is a real pleasure to use!