贝塞尔三次方的公式，涉及到4个点。如p1,p2,p3,p4，其中p1是起点，p4是重点， p2，和p3是控制点。公式如下：
B(t) = p1 * (1 - t) * (1 - t) * (1 - t) +
P2 * 3t * (1 - t) * (1 - t) +
p3 * 3t * t * (1 - t) * +
p3 * t * t * t;
那个指数不好打字，所以我把平方和立方都乘开了。t 在 [0 , 1]的范围内表示的是时间，B函数是在 t 时间的时刻，贝塞尔三次方曲线上的一个点的值。在一段连续的时间内，利用这个公式，我们就能计算出一系列的点，连接起来就是贝塞尔三次方曲线。P代表的点(x, y)，所以我带入方程的时候，需要分别令p为x和y分别计算。
为了优化，我们现在要展开这个公式，就是展开每个括号。 以后变成以 t 为系数的方程。如下
B(t) = t * t * t * (-p1 + 3p2 - 3p3 + p4) +
3t * t * (p1 - 2p2 + p3) +
3t * (-p1 + p2) +
p1;
我们可以看到，p1234点都是已知的值，t 是变量， t 的变化就能计算出点的位置。如果我们假设贝塞尔上次曲线有10个点构成，那么每个点的时间就是 t = 1 / 10, 并且这10个点的间隔是相等的0.1， 0.2， 0.3 … 这样。 于是我们只要代入 t 的值，我们就能计算出10个点的所有值。
=======================================================
我们观察到，时间 t 是倍数递增的。比如上面就是0 .1， 0.2， 0.3 直到 1.0， 步长无论取值多少都是 1 / step 倍数递增。还可以看到点p1234都是定值，只有 t 是变化的。所以我们，令
D = (-p1 + 3p2 - 3p3 + p4)
c = 3 * (p1 - 2p2 + p3)
B = 3 * (-p1 + p2)
A = p1
B(t) = A + tB + ttC + tttD
这里有一种算法，把乘法转换成加法。首先我们看一个函数如下：
f(t) = tx;
f(1) = x;
f(2) = 2x;
f(3) = 3x;
f(4) = 4x;
// 我们可以看到这个f(t) = tx函数就可以用加法来计算;
// 有几个t就循环几次迭代出 f(t) 的值
while(t--) {
x++;
}
// 那么 f(t) = t * t * x 这个函数如何计算：
f(1) = x;
3x;
f(2) = 4x; 2x;
5x;
f(3) = 9x; 2x;
7X;
f(4) = 16x;
// 中间的数值是相邻的差值, 我们看到 2x 是一个定值，于是有：
fx dfx ddfx
0 x 2x
f(1) x 3x
f(2) 4x 5x
f(3) 9x 7x
f(4) 16x
// 以上规律不言自明， fx 存储了最终f(x) = t * t * x 的值
// 我们只要设定初始的 dfx, ddfx 即可
fx = 0;
dfx = x;
ddfx = 2x;
while(t--) {
fx += dfx;
dfx += ddfx;
}
// 我们再来看看 f(t) = t * t * t * x 这个函数如何转化，原理和上面一样
f(1) = x;
7x;
f(2) = 8x; 12x;
19x; 6x;
f(3) = 27x; 18x;
37x; 6x;
f(4) = 64x; 24x;
61x;
f(5) = 125x;
// 中间的数值是相邻的差值, 我们看到 2x 是一个定值，于是有：
fx dfx ddfx dddfx
0 x 6x 6x
f(1) x 7x 12x
f(2) 8x 19x 18x
f(3) 27x 37x 24x
// fx 存储了最后 f(t) = t * t * t * x的数值，代码如下：
fx = 0;
dfx = x;
ddfx = 6x;
dddfx = 6x;
while(t--) {
fx += dfx;
dfx += ddfx;
ddfx += dddfx;
}
有了以上的算法，我们就可以B(t) = A + tB + ttC + tttD这个函数看成是 B(t) = A + f(B) + f(C) + f(D) 的新函数， 我们就可以应用把乘法转换成加法的运算了。
=======================================================
C语言实现如下：
// 定义： p1(x1, y1) p2(x2, y2) p3(x3, y3) p4(x4, y4)
// 曲线点的个数
int pointNum = 10;
float t = 1 / pointNum;
float t2 = t * t;
float t3 = t * t2;
// 预计算
float pre_3t = 3 * t;
float pre_3t2 = 3 * t2;
float pre_6t2 = 6 * t2;
float pre_6t3 = 6 * t3;
// (p1 - 2p2 + p3)
float tmp1x = x1 - x2 * 2.0 + x3;
float tmp1y = y1 - y2 * 2.0 + y3;
// (-p1 + 3p2 - 3p3 + p4)
float tmp2x = (x2 - x3) * 3.0 - x1 + x4;
float tmp2y = (y2 - y3) * 3.0 - y1 + y4;
// p1
float fx = x1;
float fy = y1;
// (p2 - p1) * 3t +
// (p1 - 2p2 + p3) * 3t2 +
// (-p1 + 3p2 - 3p3 + p4) * t3
float dfx = (x2 - x1) * pre_3t + tmp1x * pre_3t2 + tmp2x * t3;
float dfy = (y2 - y1) * pre_3t + tmp1y * pre_3t2 + tmp2y * t3;
// (p1 - 2p2 + p3) * 6t2 +
// (-p1 + 3p2 - 3p3 + p4) * 6t3
float ddfx = tmp1x * pre_6t2 + tmp2x * pre_6t3;
float ddfy = tmp1y * pre_6t2 + tmp2y * pre_6t3;
// (-p1 + 3p2 - 3p3 + p4) * 6t3
float dddfx = tmp2x * pre_6t3;
float dddfy = tmp2y * pre_6t3;
// fx, fy 就是t值对应计算出的曲线值
while(pointNum--) {
fx += dfx;
fy += dfy;
dfx += ddfx;
dfy += ddfy;
ddfx += dddfx;
ddfy += dddfy;
}
