Complex Numbers

From Winamp Developer Wiki
Revision as of 02:14, 3 October 2009 by Culix (Talk | contribs)

(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)
Jump to: navigation, search

Originally I'd planned to just write a short guide to generating a fractal using complex numbers, but complex numbers are a huge field with many applications in computer graphics. For now I'll keep everything in this article, but at some point it'll probably be split into different articles.

Introduction to Complex Numbers

The best way to describe a complex number is to say that it's a 2 dimensional number (actually complex numbers can have any number of dimensions, but we'll stick with 2 for now). It has two components, a real value and an 'imaginary' value. It's not really important to know the differences between the two components, except to know that all our standard operations (like addition, multiplication...) might work a little differently when we're working with 2 dimensional numbers.

A complex number can be represented in two ways, the first way is called the Cartesian form and looks like a + bi, where a is the real part and b is the imaginary part (the i isn't really a variable, it's just there to say b is the imaginary component). This is very handy for us working with a 2D display like Milkdrop, because it means we can think of a point on the display, P=(x,y), as a complex number, x is the real part and y is the imaginary part. In this form point P becomes a two dimensional number, and we can use it in an equation as if it were a single number (as long as we remember that the operations work differently), like this:

  F(P) = P^2 + c

which is the equation for the Julia fractal. c is a second complex number, which for us is just another float2 with arbitrary values. We'll often assign some oscillating values to c in the form of q variables that we define in the per-frame code, with equations similar to 0.5*sin(time), because this will cause the fractal to change over time. We'll come back to the Julia fractal a bit later.

So the important thing to remember is that for us a complex number is nothing more than a float2, with one variable of any equation we use being the current uv coordinates (uv.x,uv.y).

Operations with Complex Numbers

I'll give an example of each operation and how it might look in Milkdrop. The main thing to note with all these rules is that the form a + bi doesn't mean a plus bi, the + in there is just there to signify that this is a two dimensional number. So when we map this to a float2 we simply ignore the + sign, a becomes x and b becomes y.

For these examples we're going to use these two variables:

  float2 uv1 = float2(uv.x, uv.y);
  float2 c = float2(q1, q2);

Note uv1 is just our normal uv coordinates with a different name, and c is just two arbitrary values that we're importing from the per-frame code. All of the following operations will create a third float2 that would normally get assigned to a variable.

Addition: (a + bi) + (c + di) = (a + c) + (b + d)i

  uv1 + c => float2(uv1.x + c.x, uv1.y + c.y)

Subtraction: (a + bi) - (c + di) = (a - c) + (b - d)i

  uv1 - c => float2(uv1.x - c.x, uv1.y - c.y)

So addition and subtraction is pretty straight forward, it works just like adding two float2's together because that's all we're doing!

Multiplication: (a + bi) * (c + di) = (ac - bd) + (bc + ad)i

  uv1 * c => float2(uv1.x*c.x - uv1.y*c.y, uv1.y*c.x + uv1.x*c.y)

We'll see this one pop up quite a bit, because rotations on complex numbers are all done with multiplication.

Division: (a + bi) / (c + di) = ( ac + bd / c^2 + d^2 ) + ( bc - ad / c^2 + d^2 )i

  uv1 / c => float2( (uv1.x*c.x + uv1.y*c.y)/(c.x*c.x + c.y*c.y), (uv1.y*c.x + uv1.x*c.y)/(c.x*c.x + c.y*c.y) )

Division is a little messy, but on the upside we won't use it very often.

Exponents are tricky with complex numbers in Cartesian form. To avoid a complete shit storm of math we're limited to integer components, and realistically we're limited to the exponents 2 or 3. We achieve this simply by treating it as a multiplication that gets repeated for each exponent, uv1^3 = uv1*uv1*uv1. Each time we multiply we need to use the multiplication rule from above and feed the answer into another iteration of the multiplication rule. This quickly gets out of control, and there's no easy way to calculate something like uv1^1.5. For that we'll need to express complex numbers in polar form, but before we get to that lets take a look at how the Julia fractal can be computed in Milkdrop.

The Julia Fractal

The equation for the Julia fractal is z' = z^2 + c, where z is our uv coordinates and c is some arbitrary numbers that control how the resulting fractal looks. Using the operator rules from above, it's pretty easy to work this out:

  float1 zoom = 1.6;
  float2 uv1 = (uv-0.5)*zoom;
  float2 c = float2(0.2, 0.45);
  uv1 = float2(uv1.x*uv1.x - uv1.y*uv1.y, 2*uv1.x*uv1.y)+c;
  ret = tex2D(sampler_fc_main, uv1) + 0.01;

First off, there's a few things we need to say about our 3 variables,

  • zoom: It's hard to describe what this controls, but basically it decides how 'open' or 'closed' our fractal will look. Smaller values will zoom in and tend to 'open' the fractal up, while larger values will zoom out and tend to 'close' the fractal in on itself. It can be very picky and you'll probably need to experiment with each new fractal equation to find a value that works. 1.6 happens to work well for the Julia fractal.
  • uv1: (uv-0.5) is basically making the origin (point (0,0)) at the center of the screen as opposed to the top left corner. After this the screen can be thought of as a standard graph with x and y axis, both of which run from -0.5..0.5, just like the standard graphs from trigonometry class. Multiplying this by zoom adjusts the size of the graph in the x and y dimensions, since here zoom is 1.6 our graph now runs from -0.8..0.8. This isn't really an important concept, it's just important to know that you always have to adjust uv in this way for a fractal to work.
  • c: c is another picky parameter, and the range of 'acceptable' values will be different for every fractal equation.

So what we're doing in the calculation is simple, just multiplying uv1 by itself to get z^2. The y parameter of the new float2 doesn't look like the multiplication rule above, but that's just because we've simplified it from (uv1.x*uv1.y + uv1.y*uv1.x) to (2*uv1.x*uv1.y). We can do that so long as we're multiplying a complex number by itself. Then at the end we just add c, as per the Julia equation.

All the ret line is doing is sampling sampler_main and adding 0.01. This is just the standard way of displaying the fractal if we want to use sampler_main as the source, since it makes each iteration of the fractal equation a little brighter than the last, which gives the fractal its layered look.

You can plug that simple bit of code into the warp shader and you'll get the Julia fractal. Like any fractal, we could loop the equation to be done more than once per pixel, which would give us a very different (often more complex looking) version of the same fractal. It's possible though that multiple iterations will throw off the c and zoom variables so that you'll need to tweak them to get the fractal to display nicely again.

Polar Form

There are literally hundreds of fractal equations out there and not all of them are as easy to write as the Julia equation. For instance, say we wanted to try z' = z^1.5 + c, how do we translate an exponent of 1.5 into our system of multiplying complex numbers by themselves? We can't, at least I don't know how. What we need is a different way of expressing complex numbers, and luckily one exists. It's called Polar Form, and while it looks complicated in the beginning it's actually quite a bit easier to work with. But there's a trade off, working with complex numbers in polar form usually means a lot more instruction slots, because we'll be working with sin/cos/atan functions. The 64 instruction slots in PS 2.0 just won't be enough, so polar form is pretty much limited to PS 2.X or above and will still take a toll on older video cards.

Since this is meant as an introduction to complex numbers we're not going to bother figuring out why polar form works, just how. The basic idea is simple, we turn z = (a + bi) into z = r(cos(θ) + sin(θ)i). In this form, r is called the modulus (or absolute) form of z, which is calculated with the equation: |z| = sqrt(x^2 + y^2) Since for us z is just the uv coordinates, in shader code the equation looks like this:

  float1 modulus = sqrt(uv.x*uv.x + uv.y*uv.y);

You'll notice this is just the Pythagorean theorem of a right triangle. Our x and y values are the two right sides of the triangle, and the absolute form of z is just the hypotenuse of this triangle, which is the distance from the origin (the center of the screen) to the point z. This is part of the reason why we set the origin to the center of the screen earlier with (uv-0.5). The θ (theta) is called the argument (or angle) of z. To find it we need to think of our complex number z as being on a graph:

A complex number "A" on a graph.

Trig 101 tells us that to find the angle θ we need to take the arc tangent of y/x, and the shader has a built in function for this which looks like:

  float1 theta = atan2(uv.y, uv.x);

Now we have all the variables we need to put our complex number into polar form:

  float2 uv1 = float2(modulus * cos(theta), modulus * sin(theta));

Important Note: This equation is exactly the same as (uv-0.5). That is,

  float2(uv.x-0.5, uv.y-0.5) ==  float2(modulus * cos(theta), modulus * sin(theta))

as long as we've done the math for the modulus and theta right. This is very important, because it means once Milkdrop has evaluated this equation we can use uv1 in any of the Cartesian operational rules for complex numbers. In fact you can mix the two forms all you want, as long as they're done in different steps, and as long as you compute modulus and theta in the beginning of the shader code. Important Note: Any time you change the uv1 values you need to recalculate modulus and theta before you can use it again in another polar form equation.

So why bother with polar form? Well with this we can use something called De Moivre's Theorem to raise z to any power we want without using the multiplication rule. De Moivre's equation is almost the same as the normal polar form equation: File:Dem5.gif

To see it in action, lets do the equation we were having trouble with earlier, z' = z^1.5 + c:

  uv1 = float2(pow(modulus, 1.5)*cos(theta*1.5), pow(modulus, 1.5)*sin(theta*1.5))+c;

And the Julia set equation would look like this:

  uv1 = float2(pow(modulus, 2)*cos(theta*2), pow(modulus, 2)*sin(theta*2))+c;

Which gives the exact same values as

  uv1 = float2(uv1.x*uv1.x - uv1.y*uv1.y, 2*uv1.x*uv1.y)+c;

We can use polar form to multiply two complex numbers together too, as long as both are in polar form. To multiply them, we multiply the modulus of both numbers together and add the arguments (angles). So to multiply uv1*c we'd first have to calculate the modulus and theta values for c, which we'll call mod_c and theta_c, then we do this:

  uv1 = float2(modulus*mod_c*cos(theta + theta_c), modulus*mod_c*sin(theta + theta_c));

Note this method uses a LOT more instruction slots than the Cartesian form of multiplication, and it's really a waste. I've only included it here to show that any of the Cartesian operations can be done in polar form too. Normally we'd only use polar form to do exponents, because there's no simple way to do those in Cartesian form.

As a final example lets do the equation z' = (z^0.5)*c + c. Here's the full warp shader code for this equation:

   float2 zoom = 5;
   float2 c = float2(q4, q5);
   float2 uv1 = (uv-0.5)*zoom;
   float moduz = sqrt(uv1.x*uv1.x + uv1.y*uv1.y);
   float thetaz = atan2(uv1.y, uv1.x);
   uv1 = float2((pow(moduz, 0.5)*cos(thetaz*0.5)), (pow(moduz, 0.5)*sin(thetaz*0.5)));
   uv1 = float2(uv1.x*c.x - uv1.y*c.y, uv1.y*c.x + uv1.x*c.y);
   uv1 += c;

This equation doesn't give us a fractal, it's just an example of mixing the two forms together. For (hundreds) more fractal equations take a look here

Flexi also created a nice set of tutorial presets with some examples of fractals, rotations and other uses of complex numbers. Here's the link to the .zip file: Flexi's Complex Numbers Tutorial

Complex Number to Complex Power

Ok, here's where things get a little hairy. Suppose you have an equation like F(Z) = Z^(2-Z) + C , now we're raising a complex number to a complex power. So lets step through it, like most math it's not so hard once you've done it a few times. I'm going to describe the procedure from beginning to end, to give you an idea of how to approach a problem like this.

Before we start we define the regular fractal code:

   float2 zoom = 1.5;
   float2 c = float2(q4, q5);
   float2 uv1 = (uv-0.5)*zoom;
   float moduz = sqrt(uv1.x*uv1.x + uv1.y*uv1.y);
   float thetaz = atan2(uv1.y, uv1.x);

First step is to compute the easiest bit, 2-Z. Since adding/subtracting a scalar works the same for a complex number as it does for a normal number, the operation is simply:

   float2 zp = 2-uv1;

We're going to skip all the theory behind how this next equation works, it's complicated and I don't fully understand it myself. But if you want an in depth look try [this page]. The equation for finding a complex exponent of a complex number looks like this:

   x = a + bi
   y = c + di
   r = sqrt(x.x*x.x + x.y*x.y);
   θ = atan2(x.y, x.x);
   xy = ( rc * e(-d*θ) ) * (cos(c*θ + d*log(r) + isin(c*θ + d*log(r))

First thing to note is that we only need to find modulus and theta for the base, not the exponent. Second thing, the + in front of the isin is not an addition, it's the signal that this is a complex number. The other +'s are additions. So lets see this equation in shader code for computing Zzp + c:

   float2 zp = 2-uv1;
   float mp = pow(moduz, zp.x)*exp(-zp.y*thetaz);
   uv1 = mp*float2(cos(zp.y*log(moduz) + zp.x*thetaz), sin(zp.y*log(moduz) + zp.x*thetaz))+c;    

So first thing we do is find the bit at the front of the equation, which evaluates to a float1 (a scalar). exp(x) is the built in function for the funny e and is the same as ex. Then we simply plug the right values into the rest of the equation and multiply the whole float2 by mp, which is the first bit of the equation. Notice we're not using any values from uv1 directly, just the modulus and theta we calculated from it. If we get right down to it, this is actually a combination of the polar and Cartesian forms. The base (uv1) is in polar form while the exponent (zp) is in Cartesian form.

And that's that. The key with the more complicated fractal equations is to break them down into smaller pieces, so for example for the equation F(Z) = (Z^2 + C^2)^2Z the process would go something like this:

  • Calculate Z^2 (Cartesian form square), store it in a new float2
  • Calculate C^2 (Cartesian form square), store it in C
  • Calculate 2Z (just 2*Z), store it in a new float2
  • Calculate Z^2+C^2 (just simple addition), store it in Z
  • Calculate Z^2Z (complex power of a complex number), store it in Z

The most important thing here is that we're storing the results of the first and third calculations in a temporary variable. We can't store them to Z because we need Z unaltered for later equations.

Functions for Complex Math

Stahlregen found a nice use for custom funtions that can make the shader code for complex math a lot easier to read. We're going to set up 4 custom functions that you can paste into any shader and use in the shader body. The functions will be cmul, cdiv, cpow and cexp. They'll do the same thing as their regular counterparts, except they'll work with complex numbers. Here's the code to define these functions, remember to put this before the shader_body statement:

 float2 cmul(float2 mul1, float2 mul2) {
   float2 mul = float2(mul1.x*mul2.x - mul1.y*mul2.y, mul1.y*mul2.x + mul1.x*mul2.y);
   return mul;
 float2 cdiv(float2 div1, float2 div2) {
   float2 div = float2( (div1.x*div2.x + div1.y*div2.y)/(div2.x*div2.x + div2.y*div2.y),
   (div1.y*div2.x + div1.x*div2.y)/(div2.x*div2.x + div2.y*div2.y) );
   return div;
 float2 cpow(float2 base, float ex) {
   float moduz = sqrt(base.x*base.x + base.y*base.y);
   float thetaz = atan2(base.y, base.x);
   float2 pow = float2(pow(moduz, ex)*cos(thetaz*ex), pow(moduz, ex)*sin(thetaz*ex));
   return pow;
 float2 cexp(float2 base, float2 ex) {
   float moduz = sqrt(base.x*base.x + base.y*base.y);
   float thetaz = atan2(base.y, base.x);
   float mp = pow(moduz, ex.x)*exp(-ex.y*thetaz);
   float2 sol = mp*float2(cos(ex.y*log(moduz) + ex.x*thetaz), sin(ex.y*log(moduz) + ex.x*thetaz));
   return sol;

Notice all the moduz and thetaz calculations are done in these functions, which means we don't need to worry about writing them by hand anymore. All we need to do is pass the two complex variables involved to the function and it'll take care of the rest. Now the full shader_body code for the equation above, F(Z) = Z^(2-Z) + C , would look like this:

   float2 zoom = 1.5;
   float2 c = float2(q4, q5);
   float2 uv1 = (uv-0.5)*zoom;
   float2 zp = 2-uv1;
   uv1 = cexp(uv1, zp)+c; 

Nice! As a bonus, custom functions don't count toward your instruction limit if you don't use them, so there's no performance difference between using a custom function or writing the calculation out by hand each time.

Important Note: cexp() is NOT the same as the function exp(). exp() raises the funny e to a power, cexp raises one complex number to another. So if you come across an equation that uses e^Z cexp() will not work for it. More on e^Z later.