命名标记
上面定义的宏，只在匹配NumPy的数组对象定义了正确的名称时才有效。在上面的代码中，数组被命名为py_m和py_r。为了在不同的方法中使用相同的宏，NumPy数组的名称需要保持一致。
计算力
特别是与上面五行的Python代码相比，计算力数组的方法显得颇为繁琐。
static inline void compute_F(npy_int64 N, PyArrayObject *py_m, PyArrayObject *py_r, PyArrayObject *py_F) { npy_int64 i, j; npy_float64 sx, sy, Fx, Fy, s3, tmp; // Set all forces to zero. for(i = 0; i < N; ++i) { F(i, 0) = F(i, 1) = 0; } // Compute forces between pairs of bodies. for(i = 0; i < N; ++i) { for(j = i + 1; j < N; ++j) { sx = r(j, 0) - r(i, 0); sy = r(j, 1) - r(i, 1); s3 = sqrt(sx*sx + sy*sy); s3 *= s3 * s3; tmp = m(i) * m(j) / s3; Fx = tmp * sx; Fy = tmp * sy; F(i, 0) += Fx; F(i, 1) += Fy; F(j, 0) -= Fx; F(j, 1) -= Fy; } } } |
请注意，我们使用牛顿第三定律（成对出现的力大小相等且方向相反）来降低内环范围。不幸的是，它的复杂度仍然为O（N^2）。
演化函数
该文件中的最后一个函数是导出的演化方法。
static PyObject *evolve(PyObject *self, PyObject *args) { // Declare variables. npy_int64 N, threads, steps, step, i; npy_float64 dt; PyArrayObject *py_m, *py_r, *py_v, *py_F; // Parse arguments. if (!PyArg_ParseTuple(args, "ldllO!O!O!O!", &threads, &dt, &steps, &N, &PyArray_Type, &py_m, &PyArray_Type, &py_r, &PyArray_Type, &py_v, &PyArray_Type, &py_F)) { return NULL; } // Evolve the world. for(step = 0; step< steps; ++step) { compute_F(N, py_m, py_r, py_F); for(i = 0; i < N; ++i) { v(i, 0) += F(i, 0) * dt / m(i); v(i, 1) += F(i, 1) * dt / m(i); r(i, 0) += v(i, 0) * dt; r(i, 1) += v(i, 1) * dt; } } Py_RETURN_NONE; } |
在这里，我们看到了Python参数如何被解析。在该函数底部的时间步长循环中，我们看到的速度和位置向量的x和y分量的显式计算。
性能
C版本的演化方法比Python版本更快，这应该不足为奇。在上面提到的相同的i5台式机中，C实现的演化方法能够实现每秒17972个时间步长。相比Python实现，这方面有70倍的提升。
观察
注意，C代码一直保持尽可能的简单。输入参数和输出矩阵可以进行类型检查，并分配一个Python装饰器函数。删除分配，不仅能加快处理，而且消除了由Python对象不正确的引用计数造成的内存泄露（或更糟）。
下一部分
在本系列文章的下一部分，我们将通过发挥C-相邻NumPy矩阵的优势来提升这种实现的性能。之后，我们来看看使用英特尔的SIMD指令和OpenMP来进一步推进。
