Category Archives: Math

Drawing the Mandelbrot Fractal in C#

I suppose you could say I had a deprived coding childhood; I already wrote about how I never calculated pi and now I am about to tell you that until this afternoon, I had never drawn the Mandelbrot set. I don’t know why I never took the time to do so, but as with the pi calculator, when I decided to plot the Mandelbrot set, I found a very interesting coding challenge waiting for me.

The first step in the process was to fix the notoriously slow Bitmap.SetPixel() method. I haven’t decompiled System.Drawing to look at .NET’s implementation but I assume it involves some sort of Marshal.Copy()ing of unmanaged bitmap data to a managed byte array, modifications to pixels in the managed byte array, and recopying the array to unmanaged memory. That’s well and good for setting three pixels. We’re going to render Mandelbrot sets of an arbitrary size; 2800×1600 is about 4.5 million pixels. Bitmap.SetPixel() would make this process way, way longer than necessary.

I created a class with a readable, writable bitmap as well as a method to set any of its pixels using unsafe code.

public class FastBitmap
{
    public FastBitmap(int width, int height)
    {
        this.Bitmap = new Bitmap(width, height, PixelFormat.Format24bppRgb);
    }
 
    public unsafe void SetPixel(int x, int y, Color color)
    {
        BitmapData data = this.Bitmap.LockBits(new Rectangle(0, 0, this.Bitmap.Width, this.Bitmap.Height), ImageLockMode.ReadWrite, PixelFormat.Format24bppRgb);
        IntPtr scan0 = data.Scan0;
 
        byte* imagePointer = (byte*)scan0.ToPointer(); // Pointer to first pixel of image
        int offset = (y * data.Stride) + (3 * x); // 3x because we have 24bits/px = 3bytes/px
        byte* px = (imagePointer + offset); // pointer to the pixel we want
        px[0] = color.B; // Red component
        px[1] = color.G; // Green component
        px[2] = color.R; // Blue component
 
        this.Bitmap.UnlockBits(data); // Set the data again
    }
 
    public Bitmap Bitmap
    {
        get;
        set;
    }
}

This class isn’t very useful outside of this program; it would need some modifications to use pixel formats other than 24bpp RGB. However, it gets the job done! The SetPixel method is what makes this class useful. It locks the provided Bitmap’s pixels, then gets a pointer to where the first color byte of the first pixel is in memory (Scan0.ToPointer()). It then figures out what the offset is to the pixel we want to modify. Finally, it creates a pointer to the pixel we want. It sets the blue, green and red to the specified color’s blue, green and red, then finally unlocks the bitmap data. When the data gets unlocked, it is saved in memory.

The SetPixel method presented here does not provide tremendous performance. It does, however, let you work with a bitmap fast enough that you won’t get bored. Most Mandelbrot rendering time is in calculations, anyway.

Now that we have an image class, we can get down to the math. The Mandelbrot fractal is based on the following formula:

Where C is the starting point on the complex plane and the starting value is zero.

Where C is the starting point on the complex plane and the starting value is zero.

The formula is given a starting value, C, that lies somewhere in the complex plane. If the formula (using C) does not tend to infinity, then C is part of the Mandelbrot set. If C tends to infinity, it is part of the Mandelbrot set. As a quick programmatic note, if the absolute value of Z sub N (the value of the formula after N iterations) ever becomes greater than two, the formula will tend to infinity.

Knowing that, we can start coming up with some pseudo-code:

for x on the screen:
  for y on the screen:
    scale (x, y) to the complex plane to figure out where we are
    c = where we are on the screen
    z = 0 + 0i
    b = a bitmap
    for i = 0 to however many iterations:
      z = z^2 + c
      if |z| > 2:
        plot (x, y) on b as some color that indicates how much "infinite tendency" there is at C
        break

If you read the Wikipedia article about rendering the Mandelbrot set, you will notice that it uses two real numbers to accomplish the same thing as one complex one. C# has a complex number class, which we will use in our code. However, before we talk about the actual rendering, we need a color palette. The members of the Mandelbrot set will be rendered black (or actually, based on the flow of the pseudo-code, not rendered at all), while the rest of the set will be rendered some color based on how many iterations it takes for the absolute value of Z to eclipse two. The more iterations it takes, the whiter the pixel. The fewer iterations, the bluer the pixel. It all comes together in the following (very simple) method:

public static List<Color> GenerateColorPalette()
{
    List<Color> retVal = new List<Color>();
    for (int i = 0; i <= 255; i++)
    {
        retVal.Add(Color.FromArgb(255, i, i, 255));
    }
    return retVal;
}

Okay, we can finally translate the pseudo-code into a real Mandelbrot-rendering method:

public static Bitmap DrawMandelbrot(int width, int height, double rMin, double rMax, double iMin, double iMax)
{
    List<Color> Palette = GenerateColorPalette();
    FastBitmap img = new FastBitmap(width, height); // Bitmap to contain the set
 
    double rScale = (Math.Abs(rMin) + Math.Abs(rMax)) / width; // Amount to move each pixel in the real numbers
    double iScale = (Math.Abs(iMin) + Math.Abs(iMax)) / height; // Amount to move each pixel in the imaginary numbers
 
    for (int x = 0; x < width; x++)
    {
        for (int y = 0; y < height; y++)
        {
            Complex c = new Complex(x * rScale + rMin, y * iScale + iMin); // Scaled complex number
            Complex z = c;
            for (int i = 0; i < Palette.Count; i++) // 255 iterations with the method we already wrote
            {
                if (z.Magnitude >= 2.0)
                {
                    img.SetPixel(x, y, Palette[i]); // Set the pixel if the magnitude is greater than two
                    break; // We're done with this loop
                }
                else
                {
                    z = c + Complex.Pow(z, 2); // Z = Zlast^2 + C
                }
            }
        }
    }
 
    return img.Bitmap;
}

Some notable points about this method:

  1. The rScale and iScale variables represent the change in real or imaginary value in the complex plane of each pixel on the screen.
  2. Some good parameters to pass to this method are rMin = -2.5, rMax = 1.0, iMin = -1.0 and iMax (no copyright infringement intended!) = 1.0.
  3. This method could be significantly improved on multi-core machines by dividing up the x-values to render and running a thread on each core. However, any rendering would have to be done after all the threads finished to avoid any cross-thread bit-locking issues. Luckily, the part that takes a while is the iterating, not the rendering.

That’s all it takes to draw a pretty good-looking Mandelbrot fractal. Here’s a fairly large one that you can give to your math teacher. I promise it will make him or her proud. I’m totally giving one to my math team coach as a New Year’s present.

Very blue...

Very blue…

Calculating Pi

One programming task that I had never taken on (up until this morning, anyway), was a pi calculator. Finding pi is not an overwhelmingly daunting task; one can get an extremely precise approximation from the System.Math.PI constant. However, that’s no fun. The real fun is in rolling your own pi calculator. As it turns out, this was no beginner project. It was math-heavy, theory-heavy and CPU-heavy. Oh, and it was a blast!

The trouble was that at first, I had no idea where to begin. In C#, there is no “arbitrary precision” decimal type; the decimal class can hold 28 or so significant figures (remember the calculator tutorial?). No self-respecting pi calculator would ever stop at 28 significant figures. 28,000 would be a bit better! The first challenge in coding the pi calculator would be coming up with a way to store an unlimited number of significant figures. I’d read about some pi calculator authors implementing “BigFraction” or “BigRational” classes based off of the fairly new BigInteger (arbitrarily-sized .NET integer) class, so I thought I’d take a crack at writing one of those. I also decided that before I wrote a single line of code, I needed to choose a pi calculation formula.

The requirements for a formula are simple: that it converges relatively quickly and requires relatively few calculations per iteration. Unfortunately, the use of a “BigRational” class adds another requirement: no fractional powers of numbers that have irrational fractional powers. That’s a real drag, because it means we can’t use the Chudnovsky Algorithm, which is the Holy Grail of fast pi calculation algorithms:

The trouble is with the 3k + 3/2 part; 640320 is not a perfect square (it is almost 800.2^2, but not quite).

The trouble is with the 3k + 3/2 part; 640320 is not a perfect square (it is almost 800.2^2, but not quite).

With the really good algorithm out of the way, I started trying to find another good algorithm. My first thought took me back to my trig identity days a few years ago (okay, the past three years of high school):

So then if you take the deriv... actually, no.

So then if you take the derivat… actually, no.

This is actually a sound, correct thought. In fact it has a name: the Gregory-Leibniz formula. This will eventually converge to pi. However, after following a link off the handy, dandy Wikipedia page I used to decide which algorithm I would use in my program, I learned that it takes five billion iterations of arctan(1) (which I will get to in a moment) to produce ten correct digits of pi. YIKES! It was a good thought.

Anyway, I eventually decided on the original Machin formula, from 1706:

Although it was created in 1706, Wikipedia says it is still one of the fastest-converging pi-calculation algorithms!

Despite its age, it converges quickly and is relatively easy to compute.

Awesome. We have a formula. Now we need to figure out how an arctangent can be computed. It is actually rather simple, especially with programmatic help:

See. Not that bad.

See. Not that bad.

Alright! We have the formula for computing pi and we have the formula for computing arctangent, so now we’re done. Not so fast. Remember when I mentioned the BigRational class? As it turns out, writing a BigRational class is the bulk of this project. I won’t bore you with too many tiny details of the way I implemented it; I will tell you that it was designed for simplicity over speed. It definitely needs some refining, but it works.

The basic concept behind the BigRational class is simple: it is a fraction. It has a numerator and a denominator. The numerator and denominator are arbitrarily-precise, so the entire ratio is arbitrarily-precise. The class has a bunch of methods to help the numerator and denominator do their fraction-y duty:

As it turns out, a lot of help.

As it turns out, a lot of help.

Most of the methods listed in the class diagram are pretty self-explanatory; I won’t go into detail here (but you can download the source code to the project at the bottom of the page). All the operators help the BigRationals play well together, and the * operator has an overload that allows “scalar” multiplication. That is, it allows multiplication by a BigInteger.

The two methods that I will delve into are the Reduce method and the ToDecimalString method. The Reduce method reduces a BigRational into its simplest form. It implements Euclid’s Algorithm, which is a means of reducing a fraction. It works in steps:

  1. The larger of the numerator and the denominator is set as the “big number.” The smaller is set as the “small number.”
  2. The program finds the remainder when the big number is divided by the small number.
  3. If the remainder is not zero, the small number becomes the big number and the remainder from step two becomes the small number.
  4. If the remainder is zero, both the numerator and denominator may be divided by the small number. At this point, the small number is the greatest common factor of the numerator and denominator.
  5. Steps 2, 3 and 4 are repeated until step four terminates with the small number equal to one. At this point, the fraction is in simplest form.

In code, the above steps look like this:

public void Reduce()
{
    // Set up big and small number
    BigInteger bigNumber, smallNumber = 0;
    if (Numerator > Denominator)
    {
        bigNumber = Numerator;
        smallNumber = Denominator;
    }
    else
    {
        bigNumber = Denominator;
        smallNumber = Numerator;
    }
 
    // Now divide a bunch of times
    while (smallNumber != 1)
    {
        // Find the remainder
        BigInteger rem = (bigNumber % smallNumber);
        if (rem != 0)
        {
            bigNumber = smallNumber; // Set up the numbers for the next iteration
            smallNumber = rem;
        }
        else
        {
            // We can divide both the numerator and denominator by the previous remainder
            Numerator /= smallNumber;
            Denominator /= smallNumber;
 
            // Re-assign bigNumber and smallNumber
            if (Numerator > Denominator)
            {
                bigNumber = Numerator;
                smallNumber = Denominator;
            }
            else
            {
                bigNumber = Denominator;
                smallNumber = Numerator;
            }
        }
    }
 
    // One last thing -- if the numerator is positive and the denominator is negative, the numerator needs to become negative and the denominator needs to become positive.
    // If they're both negative, they should both be positive
    if (Denominator < 0 && Numerator > 0 || Numerator < 0 && Denominator < 0)
    {
        Numerator = -Numerator;
        Denominator = -Denominator;
    }
}

The only change between the pseudo-code and actual implementation was the addition of the small if-statement at the bottom. I noticed that the denominator would sometimes be negative instead of the numerator being negative. That would probably throw off some calculations, so I made sure that only the numerator of the BigRational could ever be negative.

The ToDecimalString method converts a BigRational into its decimal equivalent — a handy helper when it comes to calculating pi (since it’s no fun unless it starts with 3.14159…). I considered two different implementation styles:

  1. Multiplying the entire fraction by 10^however many digits I wanted.
  2. Creating my own long division algorithm that would return a string of digits of a specified length.

Contrary to what seems sane, I chose the second method. I suppose it was partially because I had always wanted to write an actual long division algorithm, but by any means, it took some head-scratching at first. Long division isn’t something I have to do very much, so I started by doing out a few long division problems. As it turns out, there is a nice pattern. The non-decimal part of the number may be calculated in one fell swoop by using integer division. It may then be committed to a string. After that, a decimal point should be committed to the string. Then, a “new numerator” should be created. Its value is the remainder from the previous iteration multiplied by ten. The next decimal digit is equal to the “new numerator” minus the remainder when it is divided by the denominator of the fraction, all divided by the denominator of the fraction. This process can be repeated until the desired number of digits is reached.

Okay, it makes more sense in code. Wow.

public string ToDecimalString(BigInteger decimalDigits)
{
    StringBuilder rv = new StringBuilder();
    this.Reduce(); // Go as fast as possible
    BigInteger remainder = Numerator % Denominator;
 
    // Get the non-decimal part
    rv.Append(((Numerator - remainder) / Denominator).ToString() + ".");
    remainder = Numerator % Denominator;
 
    BigInteger newNum = remainder * 10;
 
    // Now get the decimal part
    for (BigInteger i = 0; i < decimalDigits; i++)
    {
        rv.Append(((newNum - (newNum % Denominator)) / Denominator).ToString()); // This literally just does long division
        newNum = (newNum % Denominator) * 10;
    }
 
    return rv.ToString();
}

The only real advantage this method presents is the lack of any string parsing. It is all basic use of a StringBuilder. Plus, it seems to work very well!

With the BigRational class done, the next step was the actual pi calculation. The first step was to implement an arctangent method, in accordance with the formula I mentioned earlier.

public static BigRational ArcTangent(BigRational input, BigInteger iterations)
{
    BigRational retVal = input;
    for (BigInteger i = 1; i < iterations; i++)
    {
        // arctan(x) = x - x^3/3 + x^5/5 ...
        // = summation(n->infinity) (-1)^(n) * x^(2n+1)/(2n+1)
        BigRational powRat = input.Pow((2 * i) + 1);
        retVal += new BigRational(powRat.Numerator * (BigInteger)Math.Pow(-1d, (double)(i)), ((2 * i) + 1) * powRat.Denominator);
        if (i % 100 == 0)
        {
            Console.WriteLine("ArcTangent {0}: {1}/{2} iterations complete.", input, i, iterations); // Status update.
        }
    }
 
    return retVal;
}

Wait. There’s a parameter called “iterations.” How do we figure out how many iterations are necessary? This page provided the answer. As it turns out…

I find it's a good idea to add twenty or so digits to the end to make sure everything is correct.

I find it’s a good idea to add twenty or so digits to the end to make sure everything is correct.

In other words…

public static BigInteger GetConversionIterations(BigInteger digits, BigRational q)
{
    return (BigInteger)((double)digits / (2 * Math.Log10((double)q.GetReciprocal())));
}

That was quick. That also finished laying the foundation for actually calculating pi. Phew. That was a lot of keystrokes.

public static BigRational GetPi(BigInteger numDigits)
{
    // pi = 16 arctan(1/5) − 4 arctan(1/239)
    BigRational oneFifth = new BigRational(1,5);
    BigRational oneTwoThirtyNine = new BigRational(1, 239);
    BigRational arcTanOneFifth = PiCalc.ArcTangent(oneFifth, PiCalc.GetConversionIterations(numDigits + 1, oneFifth)); // Start computing
    BigRational arcTanOneTwoThirtyNine = PiCalc.ArcTangent(oneTwoThirtyNine, PiCalc.GetConversionIterations(numDigits + 1, oneTwoThirtyNine));
 
    return (arcTanOneFifth * 16) - (arcTanOneTwoThirtyNine * 4);
}

This code block implements Machin’s formula as discussed above. An improved version of this program could take advantage of multi-core computing by calculating one arctangent on one thread and the other arctangent on another. More formulas should also be tested to determine which is the fastest; this program really slows down as the iteration count gets higher because the fractions get really, really, really huge. Of course, this program was never going to break any records anyway because it was written in C#. JIT-compiled, super high-level languages aren’t exactly speed demons. Then again, my Phenom II x4 (overclocked to 3.7 gHz) is also beginning to show its age. I’d be curious to see how performance compares on one of the latest-generation i7s.

Although it is sluggish, it does seem to compute pi properly. Here are the first 10,001 digits:

Although it looks like a lot of pi, it's not a lot of pi.

Although it looks like a lot of pi, it’s not a lot of pi.

So that’s it. It’s done. It works. Now it just needs to run for a few days so I can brag to my friends! The funny thing about calculating pi is how although we believe pi is completely irrational and will never end, we’ve calculated enough digits that it might as well end. If my calculations are correct, it only takes about 80 digits of pi to calculate the volume of the observable universe to the nearest cubic meter. Thus, ten thousand would calculate it to the nearest 1/10^9920 cubic meter, whatever that would be. That’s some mind-blowing precision. Imagine what a million digits would do.

Download this Project (MIT License)

Will it Rain this Weekend?

While browsing MIT+K12 open projects, I came across an entire course proposal about probability. Now, I have very little formal probability training outside of one unit in precalculus last year, but I was intrigued by the idea of writing a program to simulate 10,000 trials. Obviously, 10,000 trials isn’t enough to get a really good data set so I added a few zeros to my C# for loop and ran 1,000,000 trials.

So, in keeping with the MIT course proposal, I checked the forecast for the weekend and found out that this weekend, the chance of rain on Saturday is 10% and the chance of rain on Sunday is 30%. My first reaction was mild aggravation because I went camping last weekend in the rain and I simply could’ve waited a week. My second reaction was, “Alright. Let’s get started.”

Probability Tree Image

Mr. Whalen, I hope you’re proud of me.

The next step was to write a program to test my calculations. I came up with a program structure that I hope is “random” enough. I am going to use System.Random instead of RNGCryptoServiceProvider because it is way faster and I believe it fine for something this simple. When I originally wrote this program, I used lists, populated them with true or false values, randomized them, and picked from them over and over again. I am not sure what kind of awful state of mind I was in when I wrote that code (whatever it was, there wasn’t enough caffeine), but I’ve come up with a way more elegant solution – Random.NextDouble(). If the double is less than the probability of rain, it rained. Otherwise, it didn’t. We also need three counters – one for when it rains sometime during the weekend at least once, one for when it rains exactly once and one for when it rains twice.

There is one other case: when it doesn’t rain, but 1,000,000 – the number of times it rained at least once takes care of that. Some operators make determining which case we’re dealing with easy – XOR is true if exactly one of two cases is true, OR is true if at least one case is true, and AND is true if both cases are true.

With that taken care of, all that was left was to write the code in a console application:

static void Main(string[] args)
{
    double probRainSat = 0.1;
    double probRainSun = 0.3;
 
    int rainyWeekends = 0;
    int oneDayRainyWeekends = 0;
    int bothDayRainyWeekends = 0;
 
    Console.WriteLine("Simulating one million weekends...");
    Console.WriteLine();
 
    Random rand = new Random(Environment.TickCount);
 
    // one million trials
    for (int i = 0; i < 1000000; i++)
    {
        if (i % 100000 == 0)
        {
            Console.WriteLine("{0}...", i);
        }
 
        bool rainedSat = rand.NextDouble() < probRainSat;
        bool rainedSun = rand.NextDouble() < probRainSun;
        if (rainedSun || rainedSat)
        {
            rainyWeekends++;
        }
        if (rainedSat ^ rainedSun)
        {
            oneDayRainyWeekends++;
        }
        else if(rainedSat && rainedSun) // We can use else if because it can't rain on exactly one and both days.
        {
            bothDayRainyWeekends++;
        }
    }
 
    Console.WriteLine("Rain on only 1 day: {0}/1,000,000 ({1})", oneDayRainyWeekends, oneDayRainyWeekends / 1000000);
    Console.WriteLine("Rain during the weekend: {0}/1,000,000 ({1})", rainyWeekends, rainyWeekends / 1000000);
    Console.WriteLine("Rain on both days: {0}/1,000,000 ({1})", bothDayRainyWeekends, bothDayRainyWeekends / 1000000);
    Console.WriteLine("Rain on neither day: {0}/1,000,000 ({1})", 1000000 - rainyWeekends, (1000000 - rainyWeekends) / 1000000);
    Console.ReadLine(); // Keep the program from closing.
}

Now you can run your new creation:

Probability Simulator Results

You just simulated about 19,200 years worth of weekends.

Not bad, huh? If you round, you’re spot on with the calculations. The math predicted almost exactly what the program would do. That, my friends, is pretty cool.