关于自制、快速 C# 数学库的漫谈
4.64/5 (11投票s)
为什么不自己构建 C# 数学库，而不是使用 DirectX 或 XNA 等 API 中包含的库？
引言
夜深了。非常晚。事实上，现在已经非常晚了，我应该说差不多是早晨了……但我沉浸在一项我拖延了太久的任务中：扩展我自己的数学库，直接用 C# 编写。
首先，我必须说，这篇文章无意成为任何严肃的东西，只是关于数学、C# 和 .Net 的一些漫谈。一些好想法，也有一些坏主意……;)
背景
为什么要写我自己的库？
因为我认为尽量少依赖图形 API 是件好事。为什么？
- 因为 API 会不时被弃用（就像托管 DirectX 那样），或者演进、改变，或者以某种形式重新出现（就像托管 DirectX 的某些部分似乎将通过 Windows 7 Code Pack 重新出现那样），或者你就是决定使用另一个 API，或者你希望你的代码同时兼容 **XNA** 和 **DirectX**。
- 因为我希望能够调试我数学的内部细节。
- 因为编写你自己的库（即使你复制了大部分内容），你也会学到很多东西。
- 因为看到事物的内部运作，你就会意识到某些任务需要付出多少努力，你会三思而后行……
还需要更多理由吗？
等一下。更快的代码，不使用 API？
对于 CPU 相关的任务，是的，这是可能的。
XNA（或 DirectX）并没有什么神奇之处。它们只是像你一样进行计算。但当然，你必须为每个任务选择最佳算法。如果你这样做，你的代码速度将与 API 中的代码一样快，但关键是，当你了解事物的内部工作原理时，你将对事物的成本有更详细的了解，并优化你的代码。此外，你还可以根据你的特定需求对事物进行一些微调，进一步提高速度。
但是，请记住，像微软这样的公司有数百甚至数千名顶尖工程师在他们的产品上工作，所以如果有一个更好的方法来做某事，API 就会拥有它。你也需要那么优秀！
让我们来看一个例子，看看你需要什么样的优化才能使你自己的数学库与 API 中包含的数学库一样快（甚至更快）。我们将包含一些矩阵运算，并对不同的算法进行一些速度测试。
关于仿射变换矩阵的说明
在 3D 图形应用程序中，变换通常由基本矩阵（这些基本矩阵是：平移、旋转、缩放和单位矩阵）的组合构成。这样构建的矩阵总是仿射变换（矩阵的右列为 0,0,0,1 – 如果使用行表示法）。
本文中的一些算法优化仅适用于仿射变换矩阵。你应该避免将它们用于非仿射矩阵。透视投影矩阵就是这种情况，它们不属于仿射变换组。
关于 JIT 编译和 .NET 数学性能的说明
这是一个非常有趣的问题：局部变量是否比成员变量更快？
起初，你可能会认为它们不是，因为它们每次调用方法时都需要重新创建。但是，正如这篇有趣的博文指出的那样，JIT 编译器可以将局部变量寄存化，但不能对成员字段做同样的处理。所以，当我们谈论数学计算时，我们应该说，是的，它们是。正如你在测试中看到的，这会带来显著的差异。所以
|
如果你在 .Net 中开发并且有一个性能关键的数学计算方法，你应该使用局部变量副本以加快计算速度。 |
事实上，如果你看一下 XNA 代码，`Matrix` 类中每个非静态的、性能关键的方法都会为每个成员字段创建一个局部副本。
第一部分：矩阵行列式计算
通用方法（基于 Graphics Gems I）
下面是一个计算矩阵行列式的相当直接的算法。它主要用于非实时应用程序，因为它相当慢且使用双精度。它基于 3x3 方法计算 4x4 行列式，然后是 2x2 等等。它基于 Graphics Gems I 中用于矩阵求逆的算法。
public double Determinant4x4()
{
double ans;
double a1, a2, a3, a4, b1, b2, b3, b4, c1, c2, c3, c4, d1, d2, d3, d4;
a1 = this.M11; b1 = this.M12;
c1 = this.M13; d1 = this.M14;
a2 = this.M21; b2 = this.M22;
c2 = this.M23; d2 = this.M24;
a3 = this.M31; b3 = this.M32;
c3 = this.M33; d3 = this.M34;
a4 = this.M41; b4 = this.M42;
c4 = this.M43; d4 = this.M44;
ans = a1 * Determinant3x3(b2, b3, b4, c2, c3, c4, d2, d3, d4)
- b1 * Determinant3x3(a2, a3, a4, c2, c3, c4, d2, d3, d4)
+ c1 * Determinant3x3(a2, a3, a4, b2, b3, b4, d2, d3, d4)
- d1 * Determinant3x3(a2, a3, a4, b2, b3, b4, c2, c3, c4);
return ans;
}
public double Determinant3x3(double a1, double a2, double a3, double b1, double b2, double b3, double c1, double c2, double c3)
{
double ans;
ans = a1 * Determinant2x2(b2, b3, c2, c3)
- b1 * Determinant2x2(a2, a3, c2, c3)
+ c1 * Determinant2x2(a2, a3, b2, b3);
return ans;
}
public double Determinant2x2(double a, double b, double c, double d)
{
double ans;
ans = a * d - b * c;
return ans;
}
通用方法（实时应用程序优化）
第二个版本也是一个完整的通用方法，但针对实时应用程序进行了优化。它快得多，因为它进行了计算，并且通过将重复的乘法分组到 numXX 变量中，节省了乘法和加法。
public float Determinant()
{
float num18 = (M33 * M44) - (M34 * M43);
float num17 = (M32 * M44) - (M34 * M42);
float num16 = (M32 * M43) - (M33 * M42);
float num15 = (M31 * M44) - (M34 * M41);
float num14 = (M31 * M43) - (M33 * M41);
float num13 = (M31 * M42) - (M32 * M41);
return ((((M11 * (((M22 * num18) - (M23 * num17)) + (M24 * num16))) - (M12 * (((M21 * num18) -
(M23 * num15)) + (M24 * num14)))) + (M13 * (((M21 * num17) - (M22 * num15)) +
(M24 * num13)))) - (M14 * (((M21 * num16) - (M22 * num14)) + (M23 * num13))));
}
这是 XNA 或 SlimDX 中使用的版本。它是事实上的标准。
优化 I（更快的内存访问）
如引言中所述，关于局部变量与成员字段的性能，我们将使用此概念对上述方法进行优化。
public float Determinant()
{
float L11 = this.M11;
float L12 = this.M12;
float L13 = this.M13;
float L14 = this.M14;
float L21 = this.M21;
float L22 = this.M22;
float L23 = this.M23;
float L24 = this.M24;
float L31 = this.M31;
float L32 = this.M32;
float L33 = this.M33;
float L34 = this.M34;
float L41 = this.M41;
float L42 = this.M42;
float L43 = this.M43;
float L44 = this.M44;
float C1 = (L33 * L44) - (L34 * L43);
float C2 = (L32 * L44) - (L34 * L42);
float C3 = (L32 * L43) - (L33 * L42);
float C4 = (L31 * L44) - (L34 * L41);
float C5 = (L31 * L43) - (L33 * L41);
float C6 = (L31 * L42) - (L32 * L41);
return ((((L11 * (((L22 * C1) - (L23 * C2)) + (L24 * C3))) - (L12 * (((L21 * C1) -
(L23 * C4)) + (L24 * C5)))) + (L13 * (((L21 * C2) - (L22 * C4)) + (L24 * C6)))) -
(L14 * (((L21 * C3) - (L22 * C5)) + (L23 * C6))));
}
性能结果请参见测试。
优化 II（基于 Graphics Gems II，仅适用于仿射变换）
以下算法是 Graphics Gems II 中仿射矩阵求逆算法代码的直接翻译（完整算法见下篇文章）。它只进行 12 次乘法和 7 次加法。请记住，此方法仅适用于仿射变换（有关更多信息，请参见引言）。
public float Determinant()
{
return (this.M11 * this.M22 * this.M33) +
(this.M12 * this.M23 * this.M31) +
(this.M13 * this.M21 * this.M32) -
(this.M13 * this.M22 * this.M31) -
(this.M12 * this.M21 * this.M33) -
(this.M11 * this.M23 * this.M32);
}
测试
好了，是时候衡量这些小家伙们了。
性能
测试在 Intel Core 2 Quad Q9450 2.66Ghz、3 Gb 内存（非多线程）上，计算 4x4 仿射矩阵行列式 1 亿次。
| 通用方法（Graphics Gems I，非实时） |
1 分钟：0.9 秒 |
| 通用方法（实时） |
2.6832 秒 |
| 优化 I（内存访问） |
2.496 秒 |
| 优化 II（仿射变换） |
1.404 秒 |
可靠性
可靠性测试使用了两个不同的矩阵。这里没有意外，所有算法都运行良好。以下是用于测试的矩阵（仿射变换）以及每个算法计算出的结果（行列式）。
矩阵 1
| 0.7071068 | -0.235702276 | 0.6666667 | 0.0 |
| 0.0 | 0.9428091 | 0.33333334 | 0.0 |
| -0.7071068 | -0.235702276 | 0.6666667 | 0.0 |
| 0.0 | 0.0 | -15.0 | 1 |
-
| 通用方法（Graphics Gems I） |
1.0000001641611829 |
| 通用方法（实时） |
1.00000012 |
| 优化 I（内存访问） |
1.00000012 |
| 优化 II（仿射变换） |
1.00000012 |
矩阵 2
| 9.848019 | 0.10061159 | 0.0 | 0.0 |
| -3.50578356 | 0.282625765 | 0.0 | 0.0 |
| 0.0 | 0.0 | 94.22 | 0.0 |
| 0.0 | 0.0 | 0.0 | 1 |
-
| 通用方法（Graphics Gems I） |
295.47639778127973 |
| 通用方法（实时） |
295.4764 |
| 优化 I（内存访问） |
295.4764 |
| 优化 II（仿射变换） |
295.4764 |
第二部分：矩阵求逆计算
最昂贵的矩阵运算之一（在常用运算中）是计算矩阵的逆。我们将从通用算法开始，稍后研究几个优化。
通用方法
这个算法（这里解释）大致如下：
public Matrix Inverse()
{
float determinant = this.Determinant();
if (determinant == 0)
throw new System.InvalidOperationException("Matrix has zero determinant (singular matrix). Inverse doesn´t exist");
Matrix b = new Matrix();
b.M11 = M22 * M33 * M44 + M23 * M34 * M42 + M24 * M32 * M43 - M22 * M34 * M43 - M23 * M32 * M44 - M24 * M33 * M42;
b.M12 = M12 * M34 * M43 + M13 * M32 * M44 + M14 * M33 * M42 - M12 * M33 * M44 - M13 * M34 * M42 - M14 * M32 * M43;
b.M13 = M12 * M23 * M44 + M13 * M24 * M42 + M14 * M22 * M43 - M12 * M24 * M43 - M13 * M22 * M44 - M14 * M23 * M42;
b.M14 = M12 * M24 * M33 + M13 * M22 * M34 + M14 * M23 * M32 - M12 * M23 * M34 - M13 * M24 * M32 - M14 * M22 * M33;
b.M21 = M21 * M34 * M43 + M23 * M31 * M44 + M24 * M33 * M41 - M21 * M33 * M44 - M23 * M34 * M41 - M24 * M31 * M43;
b.M22 = M11 * M33 * M44 + M13 * M34 * M41 + M14 * M31 * M43 - M11 * M34 * M43 - M13 * M31 * M44 - M14 * M33 * M41;
b.M23 = M11 * M24 * M43 + M13 * M21 * M44 + M14 * M23 * M41 - M11 * M23 * M44 - M13 * M24 * M41 - M14 * M21 * M43;
b.M24 = M11 * M23 * M34 + M13 * M24 * M31 + M14 * M21 * M33 - M11 * M24 * M33 - M13 * M21 * M34 - M14 * M23 * M31;
b.M31 = M21 * M32 * M44 + M22 * M34 * M41 + M24 * M31 * M42 - M21 * M34 * M42 - M22 * M31 * M44 - M24 * M32 * M41;
b.M32 = M11 * M34 * M42 + M12 * M31 * M44 + M14 * M32 * M41 - M11 * M32 * M44 - M12 * M34 * M41 - M14 * M31 * M42;
b.M33 = M11 * M22 * M44 + M12 * M24 * M41 + M14 * M21 * M42 - M11 * M24 * M42 - M12 * M21 * M44 - M14 * M22 * M41;
b.M34 = M11 * M24 * M32 + M12 * M21 * M34 + M14 * M22 * M31 - M11 * M22 * M34 - M12 * M24 * M31 - M14 * M21 * M32;
b.M41 = M21 * M33 * M42 + M22 * M31 * M43 + M23 * M32 * M41 - M21 * M32 * M43 - M22 * M33 * M41 - M23 * M31 * M42;
b.M42 = M11 * M32 * M43 + M12 * M33 * M41 + M13 * M31 * M42 - M11 * M33 * M42 - M12 * M31 * M43 - M13 * M32 * M41;
b.M43 = M11 * M23 * M42 + M12 * M21 * M43 + M13 * M22 * M41 - M11 * M22 * M43 - M12 * M23 * M41 - M13 * M21 * M42;
b.M44 = M11 * M22 * M33 + M12 * M23 * M31 + M13 * M21 * M32 - M11 * M23 * M32 - M12 * M21 * M33 - M13 * M22 * M31;
// This last operation involves a float * Matrix multiplication
return (1f / determinant) * b;
}
你可以使用本文第一部分中提供的任何方法来计算行列式。
优化 I（内存访问和部分运算节省）
此方法应用了相同的概念，但使用了局部变量并将重复的计算也打包到局部变量中。
public static Matrix InvertXNALocalValues(Matrix matrix)
{
Matrix outMatrix = new Matrix();
float L1 = matrix.M11;
float L2 = matrix.M12;
float L3 = matrix.M13;
float L4 = matrix.M14;
float L5 = matrix.M21;
float L6 = matrix.M22;
float L7 = matrix.M23;
float L8 = matrix.M24;
float L9 = matrix.M31;
float L10 = matrix.M32;
float L11 = matrix.M33;
float L12 = matrix.M34;
float L13 = matrix.M41;
float L14 = matrix.M42;
float L15 = matrix.M43;
float L16 = matrix.M44;
float L17 = (L11 * L16) - (L12 * L15);
float L18 = (L10 * L16) - (L12 * L14);
float L19 = (L10 * L15) - (L11 * L14);
float L20 = (L9 * L16) - (L12 * L13);
float L21 = (L9 * L15) - (L11 * L13);
float L22 = (L9 * L14) - (L10 * L13);
float L23 = ((L6 * L17) - (L7 * L18)) + (L8 * L19);
float L24 = -(((L5 * L17) - (L7 * L20)) + (L8 * L21));
float L25 = ((L5 * L18) - (L6 * L20)) + (L8 * L22);
float L26 = -(((L5 * L19) - (L6 * L21)) + (L7 * L22));
float L27 = 1f / ((((L1 * L23) + (L2 * L24)) + (L3 * L25)) + (L4 * L26));
float L28 = (L7 * L16) - (L8 * L15);
float L29 = (L6 * L16) - (L8 * L14);
float L30 = (L6 * L15) - (L7 * L14);
float L31 = (L5 * L16) - (L8 * L13);
float L32 = (L5 * L15) - (L7 * L13);
float L33 = (L5 * L14) - (L6 * L13);
float L34 = (L7 * L12) - (L8 * L11);
float L35 = (L6 * L12) - (L8 * L10);
float L36 = (L6 * L11) - (L7 * L10);
float L37 = (L5 * L12) - (L8 * L9);
float L38 = (L5 * L11) - (L7 * L9);
float L39 = (L5 * L10) - (L6 * L9);
outMatrix.M11 = L23 * L27;
outMatrix.M21 = L24 * L27;
outMatrix.M31 = L25 * L27;
outMatrix.M41 = L26 * L27;
outMatrix.M12 = -(((L2 * L17) - (L3 * L18)) + (L4 * L19)) * L27;
outMatrix.M22 = (((L1 * L17) - (L3 * L20)) + (L4 * L21)) * L27;
outMatrix.M32 = -(((L1 * L18) - (L2 * L20)) + (L4 * L22)) * L27;
outMatrix.M42 = (((L1 * L19) - (L2 * L21)) + (L3 * L22)) * L27;
outMatrix.M13 = (((L2 * L28) - (L3 * L29)) + (L4 * L30)) * L27;
outMatrix.M23 = -(((L1 * L28) - (L3 * L31)) + (L4 * L32)) * L27;
outMatrix.M33 = (((L1 * L29) - (L2 * L31)) + (L4 * L33)) * L27;
outMatrix.M43 = -(((L1 * L30) - (L2 * L32)) + (L3 * L33)) * L27;
outMatrix.M14 = -(((L2 * L34) - (L3 * L35)) + (L4 * L36)) * L27;
outMatrix.M24 = (((L1 * L34) - (L3 * L37)) + (L4 * L38)) * L27;
outMatrix.M34 = -(((L1 * L35) - (L2 * L37)) + (L4 * L39)) * L27;
outMatrix.M44 = (((L1 * L36) - (L2 * L38)) + (L3 * L39)) * L27;
return outMatrix;
}
优化 II（仅适用于仿射变换 - 基于 Kevin Wu 的文章，Graphics Gems II）
《Graphics Gems》丛书是任何图形或游戏开发人员书架上的必备之物。你可以在此处购买我正在谈论的书，或者此处获取有关《Graphics Gems》系列的更多信息。Google 在以下链接中提供了一些本书的摘录，但每两页就有一页缺失：Graphics gems II - Página 342。
因此，Wu 先生提出的算法如下（移植到 C#，原为 C++）。
public static bool Inverse(Matrix min, Matrix mout)
{
float det_1;
float pos, neg, temp;
double PRECISION_LIMIT = 1.0e-15;
// Calculate determinand
pos = neg = 0.0f;
temp = min.M11 * min.M22 * min.M33;
if (temp >= 0.0)
pos += temp;
else
neg += temp;
temp = min.M12 * min.M23 * min.M31;
if (temp >= 0.0)
pos += temp;
else
neg += temp;
temp = min.M13 * min.M21 * min.M32;
if (temp >= 0.0)
pos += temp;
else
neg += temp;
temp = -min.M13 * min.M22 * min.M31;
if (temp >= 0.0)
pos += temp;
else
neg += temp;
temp = -min.M12 * min.M21 * min.M33;
if (temp >= 0.0)
pos += temp;
else
neg += temp;
temp = -min.M11 * min.M23 * min.M32;
if (temp >= 0.0)
pos += temp;
else
neg += temp;
det_1 = pos + neg;
// Is the submatrix A singular? then Matrix M has no inverse
if ((det_1 == 0.0) || (Math.Abs(det_1 / (pos - neg)) < PRECISION_LIMIT))
throw new System.InvalidOperationException("Matrix has zero determinant (singular matrix). Inverse doesn´t exist");
// Calculate inverse(A) = adj(A) / det(A)
det_1 = 1.0f / det_1;
mout.M11 = (min.M22 * min.M33 - min.M23 * min.M32) * det_1;
mout.M21 = -(min.M21 * min.M33 - min.M23 * min.M31) * det_1;
mout.M31 = (min.M21 * min.M32 - min.M22 * min.M31) * det_1;
mout.M12 = -(min.M12 * min.M33 - min.M13 * min.M32) * det_1;
mout.M22 = (min.M11 * min.M33 - min.M13 * min.M31) * det_1;
mout.M32 = -(min.M11 * min.M32 - min.M12 * min.M31) * det_1;
mout.M13 = (min.M12 * min.M23 - min.M13 * min.M22) * det_1;
mout.M23 = -(min.M11 * min.M23 - min.M13 * min.M21) * det_1;
mout.M33 = (min.M11 * min.M22 - min.M12 * min.M21) * det_1;
// Calculate -C * inverse(A)
mout.M41 = -(min.M41 * mout.M11 + min.M42 * mout.M21 + min.M43 * mout.M31);
mout.M42 = -(min.M41 * mout.M12 + min.M42 * mout.M22 + min.M43 * mout.M32);
mout.M43 = -(min.M41 * mout.M13 + min.M42 * mout.M23 + min.M43 * mout.M33);
// Fill in last column
mout.M14 = mout.M24 = mout.M34 = 0.0f;
mout.M44 = 1.0f;
return true;
}
请记住，此算法仅适用于仿射变换。当然，可以使用第一部分中解释的局部变量进行一些内存访问优化。
优化 III（仅适用于由旋转和平移组成的矩阵）
如果你知道你的矩阵只由旋转和平移组成，你可以应用这种（甚至更快）的方法，它避免了计算行列式，并利用了这类矩阵的一些特性。你可以在这里阅读所有信息。
public Matrix Inverse()
{
Matrix r = new Matrix();
r.M11 = this.M11;
r.M12 = this.M21;
r.M13 = this.M31;
r.M14 = 0f;
r.M21 = this.M12;
r.M22 = this.M22;
r.M23 = this.M32;
r.M24 = 0f;
r.M31 = this.M13;
r.M32 = this.M23;
r.M33 = this.M33;
r.M34 = 0f;
r.M41 = -(r.M11 * this.M41 + r.M21 * this.M42 + r.M31 * this.M43);
r.M42 = -(r.M12 * this.M41 + r.M22 * this.M42 + r.M32 * this.M43);
r.M43 = -(r.M13 * this.M41 + r.M23 * this.M42 + r.M33 * this.M43);
r.M44 = 1f;
return r;
}
附加优化
如果需要，可以在这些算法中包含其他优化，检测特殊的矩阵，这些矩阵的逆可以更快地计算。例如对角矩阵（只有对角线需要用 1/value 进行求逆）或正交矩阵（其逆等于其转置）。
测试
我们做了一个非常简单的测试，在 Intel Core 2 Quad Q9450 2.66Ghz、3 Gb 内存（非多线程）上，进行了 50,000,000（5000 万）次仿射变换的逆计算。
| 通用方法 |
11.523 秒 |
| 优化 I（内存访问和运算节省） |
5.4912 秒 |
| 优化 II（仅适用于仿射变换） |
2.7144 秒 |
| 优化 III（仅适用于旋转和平移） |
1.326 秒 |
结论
稍微搜索一下，确实可以达到图形 API 的性能，正如上面提到的，拥有（并构建）自己的库的优势是数不胜数的。
参考文献
有关上述算法或书籍的更多信息，您可以查看以下链接：
- http://www.wintellect.com/CS/blogs/jeffreyr/archive/2004/11/13/field-access-versus-local-variable-access-performance.aspx
- http://www.acm.org/tog/GraphicsGems/
- http://www.cvl.iis.u-tokyo.ac.jp/~miyazaki/tech/teche23.html
- http://www.amazon.com/exec/obidos/ASIN/0120644819?tag=realtimerenderin
- http://books.google.es/books?id=tQn2pILgt9wC&pg=PA342&dq=graphics+gems+fast+inverse+matrix
- http://gpwiki.org/index.php/MathGem:Fast_Matrix_Inversion
保重！