画圆

和一点中学三角学知识还有一点 arch 常识

画圆和其他一些基础几何图形是常见的图形学需求，然而现代的图形 API 几乎只提供了绘制折线的功能（GDI 之类操作系统图形层的 API 可能会给，但是现在已经没人用了）。因此，有时候我们需要自己对圆进行插值，转化成折线再画到屏幕上。让我们先定义一些基础的结构体

rs
#[derive(Clone, Copy, Debug)]
pub struct Point {
    x: f64,
    y: f64,
}
pub struct Circle {
    radius: f64,
    origin: Point,
}

使用 f64 仅仅是因为现代 CPU 上 f32 和 f64 几乎一样快，实际的图形 API 几乎只有 f32 可以用。

写个插值并不难，如果你没把中学数学知识全忘掉的话

rs
impl Circle {
    pub fn interpolate(&self, step: usize) -> Vec<Point> {
        let mut res = vec![Point { x: 0.0, y: 0.0 }; step];
        for i in 0..step {
            let theta = 2.0 * std::f64::consts::PI / step as f64 * i as f64;
            res[i].x = self.origin.x + self.radius * theta.cos();
            res[i].y = self.origin.y + self.radius * theta.sin();
        }
        res
    }
}

然而这种写法在每次循环时都需要计算一次三角函数（一般都会有同时计算 sin 和 cos 的方法），这会带来一定的性能损 -- x87 里提供的fsincos指令非常慢，而 libc 里提供的 sincos 尽管经过了深度优化，但肯定还是会比简单的加减乘除慢。能不能去掉这次三角运算呢？

让我们观察第 n 次和第 n + 1 次的计算，可以看到连续的两次计算中，origin和radius都没有变，只有theta角发生了变化，每次固定增加delta，也就是说我们有

rs
res[i].x = res[i - 1].x + self.radius * ((theta_i_1 + delta).cos() - theta_i_1.cos());
res[i].y = res[i - 1].x + self.radius * ((theta_i_1 + delta).sin() - theta_i_1.sin());

这时，如果你还能多回忆起一点中学数学课知识的话，就能意识到我们可以用和差化积公式展开这堆三角函数

rs
res[i].x = res[i - 1].x + self.radius *
    (theta_i_1.cos() * delta.cos() - theta_i_1.sin() * delta.sin() - theta_i_1.cos());
res[i].y = res[i - 1].x + self.radius *
    (theta_i_1.sin() * delta.cos() + theta_i_1.cos() * delta.sin() - theta_i_1.sin());

观察这段代码，可以看到对theta_i的三角计算，被换成了对theta_i_1和delta的计算。再重复一遍，delta是个固定值，而theta_i_1则是前一次迭代中已经算好了的值。突然之间，每次循环所需的三角运算被消除掉了，只需要做几次加减乘除即可得到结果。这实际上相当于对这次运算做了个微分。

实现起来也没多麻烦

rs
impl Circle {
    pub fn interpolate_smarter(&self, step: usize) -> Vec<Point> {
        let mut res = vec![Point { x: 0.0, y: 0.0 }; step];
        let delta = 2.0 * std::f64::consts::PI / step as f64;
        let cos_d = delta.cos();
        let sin_d = delta.sin();
        let mut cos = 1.0;
        let mut sin = 0.0;
        for i in 0..256 {
            res[i].x = self.origin.x + self.radius * cos;
            res[i].y = self.origin.y + self.radius * sin;
            let next_cos = cos_d * cos - sin_d * sin;
            let next_sin = sin_d * cos + cos_d * sin;
            cos = next_cos;
            sin = next_sin;
        }
        res
    }
}

跑个分看看

test tests::normal   ... bench:       1,886.72 ns/iter (+/- 4.30)
test tests::better   ... bench:         498.40 ns/iter (+/- 0.73)

快了三倍。

这也告诉了我们要如何战胜编译器的自动优化 -- 通过编译器没掌握的领域知识。

然而这个优化能不能适用于对于直线的线性插值呢？（先别问为什么不用 GPU 插值）让我们也试试看

rs
pub struct Line {
    origin: Point,
    orient: Vector2,
}
impl Line {
    pub fn interpolate(&self, step: usize) -> Vec<Point> {
        let mut res = vec![Point { x: 0.0, y: 0.0 }; step];
        for i in 0..step {
            res[i].x = self.origin.x + self.orient.x * i as f64;
            res[i].y = self.origin.y + self.orient.y * i as f64;
        }
        res
    }
    pub fn interpolate_smarter(&self, step: usize) -> Vec<Point> {
        let mut res = vec![Point { x: 0.0, y: 0.0 }; step];
        let mut p = self.origin;
        for i in 0..step {
            res[i] = p;
            p.x += self.orient.x;
            p.y += self.orient.y;
        }
        res
    }
}

跑分结果几乎没有区别

test tests::normal_l ... bench:         117.82 ns/iter (+/- 1.50)
test tests::better_l ... bench:         119.24 ns/iter (+/- 0.29)

这是因为现代 CPU 都实现了乱序执行，如果两条指令之间没有依赖关系，而 CPU 内又有多余的计算单元，那么两条指令就可以并发执行，这叫做指令级并发（ILP）。在普通版本中，两次循环之间没有任何关系，因此可以并发执行；而“优化”过的版本中，后一次循环依赖前一次的结果。因此，即使后者比前者少掉一个乘法（SSE 这样的 SIMD 指令可以一次计算两个 f64 乘法），前者却能靠更高的指令并发度扳回性能劣势。

这也警示我们，任何性能优化的基础都是观测和跑分，不进行实际测量而空谈性能是毫无意义的。