Calculating Bezier curves fast
via forward differencing

Calculating bezier curves is quite processor intensive task, but this
algorithm will make calculating at least 10 times faster. Especially with
bezier surfaces of bezier spaces, you will need the speed.

All polynomical functions can be calculated via forward differencing.
The key idea is to start at some point of the function, move forwards
at constant step and use Taylor series to calculate the next value. With
polynomical functions, Nth derivate equals 0. Thus the Taylor series can
be cut off there and get rid of the error term. Here's the Tailor series
(Rn is the error term).

Now, what we want is to calculate the next point f(x+t), which simplifies
the formula a little:

How does this help us? For a practical example, lets derive everything
step by step for the cubic Bezier formula, for which

Since the function is 3rd order, it's 4th derivate equals 0.

Thus, we can rewrite the Taylor series as:

(1)

It calculates the exact value of f(x+t) (since all the other terms in
the Taylor series equal 0). Now all we need to know now is f(x), f '(x),
f ''(x) and f '''(x). Since bezier curves are drawn from x=0 to x=1, we
should calculate the values of these functions at x=0

Now, with these precalculations and formula (1) above (where x = 0 and
t = the point we want), we could easily calculate any point in the bezier
curve and we'd only need a couple of multiplications to get the desired
point. But lets not stop here, because we can get rid of all the slow
multiplications.

We need to have a constant step by which we move through the bezier curve.
This is t. Now look at the formula to calculate next step:

To get rid of multiplications, we need a variable for f(x), a variable
for f '(x)*t, a variable for f ''(x)*t^2 / 2 and a variable for f '''(x)*t^3
/ 6. Let's call these f, fd, fdd_per_2 and fddd_per_6, respectively.

The assignment to evaluate the next point f(x+t) becomes:

f = f + fd + fdd_per_2 + fddd_per_6

Now we can calculate the next point without multiplications in the 'loop'.
Hooray! But the problem we now have is that to calculate the next point
with this same formula, we need to update fd, fdd_per_2 and fddd_per_6
too, since they're supposed to represent f '(x)*t, f ''(x)*t^2 / 2 and
f '''(x)*t^3 / 6, and x just increased by t. Back to Taylor then.
First we need to fiqure out how to update fd. Directly from the series
follows this:

Now lets just multiply both sides by t

That can be read as an assignment:
fd = fd + ? + ?

Time to introduce new variables.

fd = fd + fdd + fddd_per_2

To calculate fdd's next value, we need:

That can be read as an assignment:
fdd = fdd + fddd

We still need to calculate fdd_per_2. Here's the non-surprising result:

That can be read as an assignment:
fdd_per_2 = fdd_per_2 + fddd_per_2

Now, there's still variables fddd, fddd_per_2 and fddd_per_6 that we
haven't updated. Well, the next derivate equals zero from there, so these
variables are in fact constants. All done, then.

Here's a C code that wraps it up pretty nicely, and you can see how all
the variables are first initialized and then updated in the loop.

The initialization needs 11 multiplications, but after that you can get
all the points you need with very low cost.

The results

To produce a 9x9 bicubic bezier surface, you'll need 143 multiplications
and 923 adds. Compared to this
Gamasutra article which introduced an algorithm with 1488 muls (!)
and 1506 adds. Thanks to that article I could do this comparsion table
for 9x9 bicubic bezier surfaces:

Algorithm name

Muls

Adds

This algorithm

143

923

Central Differencing

1488

1506

De Casteljau's algorithm

3672

3672

Direct evaluation of original formula

8586

3978

And note that this new algo produces the exact same results as the original
formula, so it is not an approximation.

Further reading

The algorithm I presented here may not be the best one, but it should
be plenty fast for many applications. If you need more, read these (of
which I've read none):