The taylor series is a mechanism for approximating a function by summing polynomial terms.

First, some definitions:

f(x)	the function to approximate
fn(x)	the nth derivative of f (note that f0(x) = f(x))
n!	n factorial, n! = n * (n - 1) * (n - 2) * ... * 2

The taylor series expansion of f(x) around a point (c) is then written:

        n -> infinity
        +----
         \      fn(c) * (x - c) n
  f(x) = +      -----------------
         /             n!
        +----
        n = 0
 

As an example, consider the case where:

f(x)	=>	exp(x)
c	=>	0

Since the derivative of exp(x) is exp(x) we know that fⁿ(x) = exp(x) for all n.

Here are the first few terms of the series:

value of n	term	value of term
0	exp(0) * (x - 0)0 / 0!	1
1	exp(0) * (x - 0)1 / 1!	x
2	exp(0) * (x - 0)2 / 2!	x2 / 2
3	exp(0) * (x - 0)3 / 3!	x3 / 6
4	exp(0) * (x - 0)4 / 4!	x4 / 24

We can write the expansion of exp(x) around c = 0 more concisely as follows:

        n -> infinity
        +----
         \      xn
  f(x) = +      --
         /      n!
        +----
        n = 0
 

Now lets look at how accurate the approximation is as we add more terms. In this case, lets use x = 3.

exp(3)	n	sum	approximation	difference
20.08553692318766774092	0	1	1	19.08553692318766774092
20.08553692318766774092	1	1 + x	4	16.08553692318766774092
20.08553692318766774092	2	1 + x + x2/2	8.5	11.58553692318766774092
20.08553692318766774092	3	1 + x + x2/2 + x3/6	13	7.08553692318766774092
20.08553692318766774092	4	1 + x + x2/2 + x3/6 + x4/24	16.375	3.71053692318766774092

As you can see, the approximation converges pretty rapidly on the actual value of exp(x).

Ok, so how do we apply this to generating a fractal? Firstly, all the statements about 'x' above apply equally well to 'z' where z is the complex number (x + iy).

In our case, assume 'current' is the point 'z' for which we wish to create an approximation of exp(z). Our algorithm (in pseudo code) looks like this:

 set actualValue   = exp(current)
 set maxIterations = 20
 set tolerance     = .00001
 set count         = 0
 set zSum          = [0, 0]
 set denominator   = 1

 while count is less than maxIterations and
       the magnitude of the difference between actualValue and zSum is
       greater than tolerance
 do
    zSum  = zSum + currentcount / denominator

    if count is greater than 0
       denominator = denominator * (count + 1)
    endif

    count = count + 1
 done
 

In this example, our expansion point (c) is zero, we can readily generalize our algorithm to handle other values of (c)

 set actualValue   = exp(current)
 set maxIterations = 20
 set tolerance     = .00001
 set count         = 0
 set zSum          = [0, 0]
 set denominator   = 1
 set c             = [2, 3]
 set cValue        = exp(c)

 while count is less than maxIterations and
       the magnitude of the difference between actualValue and zSum is
       greater than tolerance
 do
    zSum  = zSum + cValue * (current - c)count / denominator

    if count is greater than 0
       denominator = denominator * (count + 1)
    endif

    count = count + 1
 done
 

Finally, here is a full-on fractorama formula file demonstrating this technique

fractal
{
   mapping {
      (-2.33672316384180778215,
       3.41009999999999990905,
       1.18621468926553697010,
       8.00000000000000000000) => (150, 195)
   }


   formula
   {
      // The value we wish to approximate

      zvalue = exp(current);

      $xavg = 0;    
      $yavg = 0;    
      $mavg = 0;    

      $xmax = 0;    
      $ymax = 0;    
      $mmax = 0;    

      $xtot = 0;    
      $ytot = 0;    
      $mtot = 0;    
 
      $maxIterations = 15;
      $tolerance     = 1e-6;

      zSum           = 0;
      $denominator   = 1;

      // The point around which we'll create our expansion.  We use
      // exp(current) here but we could just as well have used:
      //
      // point = [2, 3]
      // point = sin(current)
      // etc.
      //
      // There is no special significance to the fact that we're using
      // exp(current) and we're attempting to approximate exp.

      point = exp(current);

      // Loop until we've been through the loop '$maxIterations' times
      // or we deem that zSum is close enough to zvalue

      while($count < $maxIterations && ssq(zvalue - zSum) > $tolerance)
      {
         zSum += exp(point) * ((current - point) ^ $count) / $denominator;

         if($count > 0) { $denominator *= ($count + 1); }

         value = cot(zSum);

         $x = mag(value ^ .1);
         $y = mag(value ^ .2);
         $m = mag(value ^ .4);

         $xavg = ($x + $xavg * $count) / ($count + 1);
         $yavg = ($y + $yavg * $count) / ($count + 1);
         $mavg = ($m + $mavg * $count) / ($count + 1);

         $xdev = sqrt(($x - $xavg) * ($x - $xavg));
         $ydev = sqrt(($y - $yavg) * ($y - $yavg));
         $mdev = sqrt(($m - $mavg) * ($m - $mavg));

         $xmax = max($xmax, $xdev);
         $ymax = max($ymax, $ydev);
         $mmax = max($mmax, $mdev);

         $xtot += $xdev;
         $ytot += $ydev;
         $mtot += $mdev;
      }


      $r = 0;
      $g = 0;
      $b = 0;

      if($count > 0)
      {
         $r = ($xtot/$count) / $xmax;
         $g = ($ytot/$count) / $ymax;
         $b = ($mtot/$count) / $mmax;

         $r = get_sin_color($r, 255, 1);
         $g = get_sin_color($g, 255, 1);
         $b = get_sin_color($b, 255, 1);
      }

      set_color($r, $g, $b);
   }
}

And here is the resulting image: