45个发布版本

0.9.2 2023年9月3日
0.9.1 2023年3月19日
0.9.0 2022年3月3日
0.8.1 2021年5月21日
0.3.3 2019年11月29日

#26 in 数学

Download history 3557/week @ 2024-05-01 3243/week @ 2024-05-08 3112/week @ 2024-05-15 3242/week @ 2024-05-22 3502/week @ 2024-05-29 3176/week @ 2024-06-05 2954/week @ 2024-06-12 3064/week @ 2024-06-19 2850/week @ 2024-06-26 1199/week @ 2024-07-03 3042/week @ 2024-07-10 3048/week @ 2024-07-17 3219/week @ 2024-07-24 2673/week @ 2024-07-31 4074/week @ 2024-08-07 2583/week @ 2024-08-14

12,951 每月下载量
用于 89 个crate(24直接使用)

MIT OR Apache-2.0 OR Zlib

440KB
9K SLoC

crates.io docs.rs

紫外线

这是一个用于计算机图形和游戏相关的线性代数和几何代数的crate,但非常快速,无论是从生产力还是从运行时性能来看。

从生产力的角度来看,紫外线没有使用泛型,并设计为尽可能直观的接口,从而实现快速的编译时间和清晰的代码。此外,没有泛型和Rust类型系统“黑客”技术,导致错误清晰简洁,用户易于解析和修复。

从运行时性能来看,紫外线从一开始就注重性能。为此,我们为每种类型提供两种不同的类型,每种类型都具有几乎相同的功能,一种使用常规的f32标量值,另一种使用SIMD f32x4向量。这种设计意图明确,并且允许代码充分利用SIMD。

“宽”类型使用“SoA”(数组结构)架构,这样每个宽数据结构实际上包含与其关联的数据类型的4个或8个数据,并将同时执行任何操作的所有SIMD 'lanes'。例如,一个 Vec3x8 等同于捆绑在一起的一个数据结构中的8个 Vec3

这样做可能比标准的“AoS”(结构体数组)布局快得多(速度快10倍),但这取决于您的工作负载和算法需求。算法必须精心设计才能充分利用这一点,而且做到这一点可能比说的要难,特别是如果您的算法涉及大量分支。

ultraviolet 是第一个以 "AoSoA" 方式设计的 Rust 数学库,尽管现在 nalgebra 也支持其几个数据结构。

基准测试

请查看 mathbench-rs 获取最新的基准测试结果(可能不是完全与 git master 一致)。

Cargo 功能

为了进一步提高构建时间,ultraviolet 将各种功能置于功能标志之下。例如,二维和三维射影几何代数以及 f64 和整数类型默认禁用。为了启用它们,请在您的 Cargo.toml 中启用相应的 crate 功能标志。例如

[dependencies]
ultraviolet = { version = "0.9", features = [ "f64", "int" ] }

将启用 f64int 功能。以下是可用的功能列表

  • f64 – 启用 f64 位宽的浮点数支持。命名约定为 D[Type],例如 DVec3x4 将是一个包含 4 个三维向量,每个向量具有 f64 精度的集合。
  • int – 启用整数向量类型。
  • bytemuck – 启用将许多类型转换为字节数组的转换,用于与图形 API 一起使用。
  • mint – 通过 mint 接口启用与其他数学 crate 的互操作。
  • num-traits – 启用与其他数学 crate 互操作的身份特征。
  • serde – 为许多标量类型启用 SerializeDeserialize 实现功能。

Crate 功能

该 crate 目前正在我的光线追踪器 rayn 中进行测试,并被各种独立的 Rust 游戏开发者用于各种项目。它完成了用户目前所需的所有功能。

这个库有几个相对独特/新颖的功能,其中最重要的是几何代数的使用。

我们不是实现复数代数(用于二维旋转)和四元数代数(用于三维旋转),而是使用来自几何代数的概念——旋量,来表示二维和三维旋转。

对程序员来说,这意味着您将使用 Rotor3 类型来代替四元数,尽管您可以期望它基本上能做四元数能做的所有事情。实际上,四元数与旋量直接同构(这意味着它们在本质上是一样的,只是表达方式不同)。做出这个决定有两个原因:首先,数学推导实际上非常简单易懂。在这个库中实现的所有 Rotor 结构的数学推导都写在 GitHub 仓库的 derivations 文件夹中;我手动推导了它们作为实现的一部分。

另一方面,四元数通常被视为我们程序员用来进行旋转的黑盒,因为它们具有一些不错的特性,但我们并不真正理解。您也可以用同样的方式使用旋转器,但您可以轻松地理解它们。第二,从某种意义上说,它们可以被视为比四元数“更正确”。具体来说,它们有助于更正确地理解旋转是发生在平面内而不是围绕发生的事情,正如通常所想的那样。最后,旋转器也可以推广到4维甚至更高维,如果有人想的话,他们可以实施一个Roter4,它保留Rotor3/四元数的所有特性,但在4维中进行旋转,这是四元数根本无法做到的。

如果您发现缺少了您需要的功能,请在GitHub问题跟踪器和/或Rust社区discord服务器(我在那里是Fusha)上向我报告,我会尝试为您添加,如果我认为它符合lib的愿景的话 :)

示例

欧拉积分

欧拉积分是一种用于数值求解常微分方程的方法。如果这听起来很复杂,不用担心!如果您不打算实施任何类型的物理模拟,这个方法在游戏中很常见。继续阅读下面的代码!

关键是,如果您在多个浮点值上执行相同的数学运算而没有条件(没有if),则将它们移植到宽数据类型和并行处理相当简单。

以下是欧拉积分的标量示例

fn integrate(
    pos: &mut [uv::Vec3],
    vel: &mut [uv::Vec3],
    acc: &[uv::Vec3],
    dt: f32,
) {
    for ((position, velocity), acceleration) in pos.iter_mut().zip(vel).zip(acc) {
        *velocity = *velocity + *acceleration * dt;
        *position = *position + *velocity * dt;
    }
}

代码遍历每对相应的位置、速度和加速度向量。它首先通过已过去的时间调整速度,然后通过速度乘以过去的时间调整位置。

这些都是需要以相同方式应用到所有相关变量中的乘法、加法和赋值运算符。

要将此函数移植到宽数据类型和并行处理,我们只需更改数据类型即可!新的函数看起来像这样

fn integrate_x8(
    pos: &mut [uv::Vec3x8],
    vel: &mut [uv::Vec3x8],
    acc: &[uv::Vec3x8],
    dt: f32x8,
) {
    for ((position, velocity), acceleration) in pos.iter_mut().zip(vel).zip(acc) {
        *velocity = *velocity + *acceleration * dt;
        *position = *position + *velocity * dt;
    }
}

此函数现在并行处理8组向量,带来了显著的加速!

唯一的缺点是创建向量切片的调用代码需要修改,以便用8组值而不是一组值填充这些宽数据类型。该标量代码如下所示

let mut pos: Vec<uv::Vec3> = Vec::with_capacity(100);
let mut vel: Vec<uv::Vec3> = Vec::with_capacity(100);
let mut acc: Vec<uv::Vec3> = Vec::with_capacity(100);

// You would probably write these constant values in-line but
// they are here for illustrative purposes
let pos_x = 1.0f32;
let pos_y = 2.0f32;
let pos_z = 3.0f32;

let vel_x = 4.0f32;
let vel_y = 5.0f32;
let vel_z = 6.0f32;

let acc_x = 7.0f32;
let acc_y = 8.0f32;
let acc_z = 9.0f32;

for ((position, velocity), acceleration) in pos.iter_mut().zip(vel).zip(acc) {
    pos.push(uv::Vec3::new(pos_x, pos_y, pos_z));
    vel.push(uv::Vec3::new(vel_x, vel_y, vel_z));
    acc.push(uv::Vec3::new(acc_x, acc_y, acc_z));
}

而要为8车道宽的Vec3x8数据类型填充相同的内容,代码可能看起来像这样

let mut pos: Vec<uv::Vec3x8> = Vec::with_capacity(100 / 8 + 1);
let mut vel: Vec<uv::Vec3x8> = Vec::with_capacity(100 / 8 + 1);
let mut acc: Vec<uv::Vec3x8> = Vec::with_capacity(100 / 8 + 1);

let pos_x = uv::f32x8::splat(1.0f32);
let pos_y = uv::f32x8::splat(2.0f32);
let pos_z = uv::f32x8::splat(3.0f32);

let vel_x = uv::f32x8::splat(4.0f32);
let vel_y = uv::f32x8::splat(5.0f32);
let vel_z = uv::f32x8::splat(6.0f32);

let acc_x = uv::f32x8::splat(7.0f32);
let acc_y = uv::f32x8::splat(8.0f32);
let acc_z = uv::f32x8::splat(9.0f32);

for ((position, velocity), acceleration) in pos.iter_mut().zip(vel).zip(acc) {
    pos.push(uv::Vec3x8::new(pos_x, pos_y, pos_z));
    vel.push(uv::Vec3x8::new(vel_x, vel_y, vel_z));
    acc.push(uv::Vec3x8::new(acc_x, acc_y, acc_z));
}

请注意,数学意义上的100 / 8相当于12.5,但我们不能方便地有一个半大小的Vec3x8

处理这些“余数”向量的方法有很多。您可以选择退回到标量代码,或者逐步退回到更窄的宽类型,例如Vec3x4,或者您可以考虑计算一些额外向量(您不会使用)的成本是否值得在代码中添加复杂性。

光线-球体相交

一次操作一个值的标量代码需要重构以利用SIMD和4-/8宽数据类型。

以下是使用Vec3进行点矢量操作的标量光线-球体相交代码示例

fn ray_sphere_intersect(
    ray_o: uv::Vec3,
    ray_d: uv::Vec3,
    sphere_o: uv::Vec3,
    sphere_r_sq: f32,
) -> f32 {
    let oc = ray_o - sphere_o;
    let b = oc.dot(ray_d);
    let c = oc.mag_sq() - sphere_r_sq;
    let descrim = b * b - c;

    if descrim > 0.0 {
        let desc_sqrt = descrim.sqrt();

        let t1 = -b - desc_sqrt;
        if t1 > 0.0 {
            t1
        } else {
            let t2 = -b + desc_sqrt;
            if t2 > 0.0 {
                t2
            } else {
                f32::MAX
            }
        }
    } else {
        f32::MAX
    }
}

本移植指南将不会讨论算法的细节,而是关注如何将代码转换为应用宽数据类型的并行SIMD操作。

首先,需要将参数和返回类型从标量 Vec3 转换为宽 Vec3x8f32x8

fn ray_sphere_intersect_x8(
    ray_o: uv::Vec3x8,
    ray_d: uv::Vec3x8,
    sphere_o: uv::Vec3x8,
    sphere_r_sq: uv::f32x8,
) -> uv::f32x8 {

每次函数调用将并行处理8个光线-球体交点。函数的前四行保持不变

    let oc = ray_o - sphere_o;
    let b = oc.dot(ray_d);
    let c = oc.mag_sq() - sphere_r_sq;
    let descrim = b * b - c;

尽管这段代码相同,但8条光线和球体的计算将同时进行!

下一行标量代码测试 descrim 的值,看它是否大于 0.0。当一次操作8个值时,代码不能沿两条不同的路径分支,因为8个值中的每个 descrim 的值可能会导致分支到不同的操作集合。为了支持这一点,我们需要将代码转换回标量代码,然后我们会失去并行处理的所有性能优势。

那么,我们应该如何转换呢?我们需要考虑一个权衡,这取决于发散的频率,也就是说,取决于分支沿一条或多条路径的频率。如果对于给定数据和算法,大多数分支将沿一条路径进行是非常可能的,那么我们可以检查是否所有通道都沿该路径进行,然后根据这个结果进行分支。这种倾向于一个分支路径的情况相对较少,对于这个算法来说,分支两种情况都很常见,所以这种方法会产生较慢的代码。

另一种方法是计算所有8个通道的分支结果,然后使用掩码过滤最终从可能的结果中选择正确的值。

为了创建具有 0.0 的8个 descrim 值的掩码

    let desc_pos = descrim.cmp_gt(uv::f32x8::splat(0.0));

在原始标量版本的真正情况下,我们然后有更多的算术运算,当我们在矢量化版本上执行时,它们看起来完全相同

    let desc_sqrt = descrim.sqrt();

    let t1 = -b - desc_sqrt;

现在在标量代码中,我们基于 t1 > 0.0 的另一个分支,所以我们应用相同的技巧,但稍微多了一点

    let t1_valid = t1.cmp_gt(uv::f32x8::splat(0.0)) & desc_pos;

末尾的 & desc_pos 执行位与操作,以组合表示每个 t1 > 0.0 通道是真是假的掩码,以及每个 descrim > 0.0 通道是真是假的掩码,如果两个通道都为真,则 t1_mask 中的该通道掩码值将为真,否则该通道的值将为 false。这是组合嵌套逻辑。

t1 > 0.0 条件的真值仅返回 t1,但假值有一些额外的计算和分支,可以以类似的方式进行移植

    let t2 = -b + desc_sqrt;
    let t2_valid = t2.cmp_gt(uv::f32x8::splat(0.0)) & desc_pos;

这听起来可能比标量代码慢,因为将此算法应用于宽数据类型会为所有分支进行计算,而不管哪个分支是正确的,你是对的!

这种方法确实是一种权衡,取决于分支方向的可能性和分支计算的成本。然而,即使是像我们在分析中的光线-球体相交这样的特别分支密集型算法,在实践中,能够同时计算多个数据的好处往往会使整体收益增加!就像所有的优化一样,测量才是真理。

到目前为止,我们几乎已经移植了整个算法。对于每个8个通道,我们都有了t1t2的值。在t1_valid中,我们有掩码值,表示每个通道是否满足descrim > 0.0 && t1 > 0.0。同时,我们还有t2_valid,其值精确指示descrim > 0.0 && t2 > 0.0。当标量代码不返回t1t2时,它返回f32::MAX。我们如何为每个通道选择正确的返回值呢?

ultraviolet在掩码类型上有一个blend函数,它使用每个通道的真实或假值来从计算的真实和假值中选择。所以如果a是计算分支真实情况的宽向量值,而b是假情况,使用掩码m,我们可以通过调用m.blend(a, b)ab中选择,根据m得到期望的输出值!

让我们尝试通过查看其逻辑控制流来将该技术应用于标量代码。

    if descrim > 0.0 {
        if t1 > 0.0 {
            t1
        } else {
            if t2 > 0.0 {
                t2
            } else {
                f32::MAX
            }
        }
    } else {
        f32::MAX
    }

所以如果我们取最外层的if条件...

   let t = t1_valid.blend(t1, ???);

descrim > 0.0 && t1 > 0.0测试的假值是什么?有两种可能性 - 要么是descrim <= 0.0,这是descrim > 0.0条件的假值,要么是descrim > 0.0 && t1 <= 0.0,这是我们处理t2的else情况。这看起来很复杂。让我们看看标量代码中的descrim > 0.0 && t2 > 0.0情况,并尝试使用blend

    let t = t2_valid.blend(t2, uv::f32x8::splat(std::f32::MAX));

因此,当descrim > 0.0 && t2 > 0.0时,存在两种错误情况,要么是descrim <= 0.0,我们希望返回f32::MAX,要么是descrim > 0.0 && t2 <= 0.0,我们希望返回f32::MAX,因此我们可以使用blend来选择正确的值以覆盖descrim > 0.0条件下的错误情况,以及t1 > 0.0条件下的错误情况,这样就只剩下解决t1 > 0.0条件的正确情况...

这正是t1_valid.blend(t1, ???)所要选择的!因此,我们可以将这两个混合操作结合起来,如下所示

    let t = t2_valid.blend(t2, uv::f32x8::splat(std::f32::MAX));
    let t = t1_valid.blend(t1, t);

t现在包含t1t2f32::MAX,具体取决于每个通道的情况!我们已经完成了将标量算法代码移植到利用8通道宽数据类型的SIMD操作,以并行计算8个光线-球体交点!

下面是使用宽Vec3x8类型实现的相同光线-球体交点算法的完整示例

fn ray_sphere_intersect_x8(
    sphere_o: uv::Vec3x8,
    sphere_r_sq: uv::f32x8,
    ray_o: uv::Vec3x8,
    ray_d: uv::Vec3x8,
) -> uv::f32x8 {
    let oc = ray_o - sphere_o;
    let b = oc.dot(ray_d);
    let c = oc.mag_sq() - sphere_r_sq;
    let descrim = b * b - c;

    let desc_pos = descrim.cmp_gt(uv::f32x8::splat(0.0));

    let desc_sqrt = descrim.sqrt();

    let t1 = -b - desc_sqrt;
    let t1_valid = t1.cmp_gt(uv::f32x8::splat(0.0)) & desc_pos;

    let t2 = -b + desc_sqrt;
    let t2_valid = t2.cmp_gt(uv::f32x8::splat(0.0)) & desc_pos;

    let t = t2_valid.blend(t2, uv::f32x8::splat(std::f32::MAX));
    let t = t1_valid.blend(t1, t);

    t
}

依赖关系

~1.5MB
~27K SLoC