问题定义

设线段 AA 由顶点 (Asx,Asy)(A_{sx}, A_{sy}) 和顶点 (Aex,Aey)(A_{ex}, A_{ey}) 组成，记为 (Asx,Asy)  (Aex,Aey)‾\overline{(A_{sx}, A_{sy})\;(A_{ex}, A_{ey})}。

设线段集 BB 包含 nn 条线段：

B={(Bsx[0],Bsy[0])  (Bex[0],Bey[0])‾,(Bsx[1],Bsy[1])  (Bex[1],Bey[1])‾,…,(Bsx[n−1],Bsy[n−1])  (Bex[n−1],Bey[n−1])‾}B = \left\{ \overline{(B_{sx}[0], B_{sy}[0])\;(B_{ex}[0], B_{ey}[0])}, \overline{(B_{sx}[1], B_{sy}[1])\;(B_{ex}[1], B_{ey}[1])}, \dots, \overline{(B_{sx}[n-1], B_{sy}[n-1])\;(B_{ex}[n-1], B_{ey}[n-1])} \right\}

目标：线段 AA 与线段集 BB 中的每条线段求交，得到以下信息：

  1. 交点坐标数组：{(Rx[0],Ry[0]),…,(Rx[n−1],Ry[n−1])}\{(R_x[0], R_y[0]), \dots, (R_x[n-1], R_y[n-1])\}
  2. 直线相交判定数组：Rintersection[n]R_{intersection}[n]，为 true\text{true} 表示 AA 所在直线与 BiB_i 所在直线存在交点，且交点落在两条线段上
  3. 线段重叠判定数组：Roverlap[n]R_{overlap}[n]，为 true\text{true} 表示 AA 与 BiB_i 共线且有重叠部分
  4. 线段集大小 size=n\text{size} = n

本文从纯标量实现出发，逐步引入 AVX2（256-bit，一次处理 4 个 double）和 AVX-512F（512-bit，一次处理 8 个 double）指令进行矢量化优化，并介绍两项关键优化——无除法相交判定和重叠计算跳过——最终在 AMD Zen4 平台（WSL2 环境）上完成性能验证。

数学基础

直线参数方程

线段上的任意一点可表示为参数方程。对于线段 AA：

{cx=Asx+A⃗sex⋅scy=Asy+A⃗sey⋅s\begin{cases} c_x = A_{sx} + \vec{A}_{sex} \cdot s \\ c_y = A_{sy} + \vec{A}_{sey} \cdot s \end{cases}

其中 A⃗sex=Aex−Asx\vec{A}_{sex} = A_{ex} - A_{sx}，A⃗sey=Aey−Asy\vec{A}_{sey} = A_{ey} - A_{sy}，参数 s∈[0,1]s \in [0, 1]。

同理，线段 BiB_i 上的任意一点：

{cx=Bsx[i]+B⃗sex[i]⋅tcy=Bsy[i]+B⃗sey[i]⋅t\begin{cases} c_x = B_{sx}[i] + \vec{B}_{sex}[i] \cdot t \\ c_y = B_{sy}[i] + \vec{B}_{sey}[i] \cdot t \end{cases}

其中 B⃗sex[i]=Bex[i]−Bsx[i]\vec{B}_{sex}[i] = B_{ex}[i] - B_{sx}[i]，B⃗sey[i]=Bey[i]−Bsy[i]\vec{B}_{sey}[i] = B_{ey}[i] - B_{sy}[i]，参数 t∈[0,1]t \in [0, 1]。

求交点参数

联立两参数方程：

{Asx+A⃗sex⋅s=Bsx[i]+B⃗sex[i]⋅tAsy+A⃗sey⋅s=Bsy[i]+B⃗sey[i]⋅t\begin{cases} A_{sx} + \vec{A}_{sex} \cdot s = B_{sx}[i] + \vec{B}_{sex}[i] \cdot t \\ A_{sy} + \vec{A}_{sey} \cdot s = B_{sy}[i] + \vec{B}_{sey}[i] \cdot t \end{cases}

用 Cramer 法则求解。行列式：

Δ=A⃗sey⋅B⃗sex−A⃗sex⋅B⃗sey\Delta = \vec{A}_{sey} \cdot \vec{B}_{sex} - \vec{A}_{sex} \cdot \vec{B}_{sey}

当 Δ≠0\Delta \neq 0 时两条直线不平行，存在唯一交点。参数 ss 的分子：

s1=Asx⋅B⃗sey+Bsx⋅(Asy−Bey)+Bex⋅(Bsy−Asy)s_1 = A_{sx} \cdot \vec{B}_{sey} + B_{sx} \cdot (A_{sy} - B_{ey}) + B_{ex} \cdot (B_{sy} - A_{sy})

参数 tt 的分子：

t1=Asx⋅(Aey−Bsy)+Aex⋅(Bsy−Asy)−Bsx⋅A⃗seyt_1 = A_{sx} \cdot (A_{ey} - B_{sy}) + A_{ex} \cdot (B_{sy} - A_{sy}) - B_{sx} \cdot \vec{A}_{sey}

则 s=s1/Δs = s_1 / \Delta，t=t1/Δt = t_1 / \Delta。

相交与重叠判定

  • 直线相交：Δ≠0  ∧  s∈[0,1]  ∧  t∈[0,1]\Delta \neq 0 \;\land\; s \in [0,1] \;\land\; t \in [0,1]
  • 线段重叠：Δ=0\Delta = 0（两直线平行），且至少一个端点落在对方线段上

其中"点是否在线段上"的判定函数 contain：已知线段起点和方向向量，通过反解参数 tx=(cx−startx)/dir⃗xt_x = (c_x - \text{start}_x) / \vec{\text{dir}}_x，检查 t∈[0,1]t \in [0,1] 且两个坐标方向解出的 tt 一致。

关键洞察：无除法判定 s∈[0,1]s \in [0,1]

这是本文最核心的优化。s=s1/Δs = s_1 / \Delta，要判定 s≥0∧s≤1s \geq 0 \land s \leq 1，可以完全避免除法：

Δ>0:s≥0  ⟺  s1≥0,s≤1  ⟺  s1≤ΔΔ<0:s≥0  ⟺  s1≤0,s≤1  ⟺  s1≥Δ\begin{aligned} \Delta > 0 &: \quad s \geq 0 \iff s_1 \geq 0,\qquad s \leq 1 \iff s_1 \leq \Delta \\ \Delta < 0 &: \quad s \geq 0 \iff s_1 \leq 0,\qquad s \leq 1 \iff s_1 \geq \Delta \end{aligned}

同理适用于 tt。这一变换将除法（延迟 ~14–20 CPU 周期）替换为比较（延迟 ~3 周期），对于绝大多数不相交的线段对（实际测试中 ≥99.99%\ge 99.99\%），完全消除了除法的开销。

完整代码

头文件与辅助函数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
#include <immintrin.h>   // Intel SIMD intrinsics (AVX2 + AVX-512F)
#include <math.h>         // 浮点数学函数
#include <stdio.h>        // 文件 I/O
#include <stdlib.h>       // 通用工具
#include <sys/time.h>     // 计时

#include <cstring>        // memcpy
#include <functional>     // std::function
#include <iostream>       // std::cout
#include <vector>         // std::vector
using namespace std;

// 判定点 (cx, cy) 是否在线段 A 上
// 线段 A: 起点 (Asx, Asy)，方向向量 (Asex, Asey)
bool contain(double Asx, double Asy, double Asex, double Asey,
             double cx,  double cy)
{
    // 线段 A 上的点可表示为: cx = Asx + Asex * t,  cy = Asy + Asey * t
    // 分别从 x 和 y 坐标反解出参数 t
    double tx = (cx - Asx) / Asex;  // 从 x 坐标反解参数 t
    double ty = (cy - Asy) / Asey;  // 从 y 坐标反解参数 t
    // t 必须两方向一致（说明点在直线上）且落在 [0,1] 内（说明点在线段上）
    return tx == ty && 0 <= tx && tx <= 1;
}

串行版本（含无除法判定优化）

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
// =====================================================================
// 串行线段求交（标量版本，含无除法判定优化）
//
// 核心变换:
//   s = s1 / Δ ∈ [0,1]  ⇔
//     (Δ > 0  && 0 <= s1 && s1 <= Δ) ||
//     (Δ < 0  && Δ <= s1 && s1 <= 0)
//
// 这一变换消除了 99.996% 线段对上的除法运算，
// 只有真正相交的线段对才执行除法计算交点坐标。
// =====================================================================
void intersectVector(
    double Asx,           // 线段 A 起点 x
    double Asy,           // 线段 A 起点 y
    double Aex,           // 线段 A 终点 x
    double Aey,           // 线段 A 终点 y
    double *Barrsx,       // 线段集 B 起点 x 数组
    double *Barrsy,       // 线段集 B 起点 y 数组
    double *Barrex,       // 线段集 B 终点 x 数组
    double *Barrey,       // 线段集 B 终点 y 数组
    double *Rarrx,        // [输出] 交点 x 坐标数组
    double *Rarry,        // [输出] 交点 y 坐标数组
    bool *Rintersection,  // [输出] 是否直线相交
    bool *Roverlap,       // [输出] 是否共线重叠
    int size)             // 线段集 B 的大小
{
    // 线段 A 的方向向量（起点 → 终点）
    double Asex = Aex - Asx;
    double Asey = Aey - Asy;

    for (int i = 0; i < size; i++)
    {
        // ---- 1. 加载线段 B_i 的起点、终点 ----
        double Bsx = Barrsx[i], Bsy = Barrsy[i];
        double Bex = Barrex[i], Bey = Barrey[i];

        // ---- 2. 线段 B_i 的方向向量 ----
        double Bsex = Bex - Bsx;
        double Bsey = Bey - Bsy;

        // ---- 3. 求行列式（分母）—— Δ = A_sey·B_sex − A_sex·B_sey ----
        // 若 Δ == 0，则 A 与 B_i 平行或共线
        double denom = Asey * Bsex - Asex * Bsey;

        if (denom == 0)
        {
            // ---- 平行/共线分支 ----
            // 不可能有普通交点
            Rintersection[i] = false;
            Rarrx[i] = 0; Rarry[i] = 0;

            // 判定是否重叠：检查四个端点是否至少有一个落在对方线段上
            Roverlap[i] =
                (contain(Asx, Asy, Aex, Aey, Bsx, Bsy) ||  // B起点是否在A上
                 contain(Asx, Asy, Aex, Aey, Bex, Bey) ||  // B终点是否在A上
                 contain(Bsx, Bsy, Bex, Bey, Asx, Asy) ||  // A起点是否在B上
                 contain(Bsx, Bsy, Bex, Bey, Aex, Aey));   // A终点是否在B上
        }
        else
        {
            // ---- 不平行分支：无除法判定相交 ----
            Roverlap[i] = false;  // Δ != 0 不可能重叠

            // 计算参数 s 的分子 s₁ 和参数 t 的分子 t₁
            // s₁ = A_sx·B_sey + B_sx·(A_sy−B_ey) + B_ex·(B_sy−A_sy)
            double s1 = Asx * Bsey + Bsx * (Asy - Bey) + Bex * (Bsy - Asy);

            // t₁ = A_sx·(A_ey−B_sy) + A_ex·(B_sy−A_sy) − B_sx·A_sey
            double t1 = Asx * (Aey - Bsy) + Aex * (Bsy - Asy) - Bsx * Asey;

            // ★ 核心优化: 用乘法+比较替代除法判定 s ∈ [0,1] 和 t ∈ [0,1]
            //   s = s₁/Δ, t = t₁/Δ
            //   若 Δ > 0: s ∈ [0,1] ⇔ 0 <= s₁ <= Δ
            //   若 Δ < 0: s ∈ [0,1] ⇔ Δ <= s₁ <= 0
            bool s_ok, t_ok;
            if (denom > 0)
            {
                s_ok = (s1 >= 0 && s1 <= denom);
                t_ok = (t1 >= 0 && t1 <= denom);
            }
            else  // Δ < 0（已知 Δ != 0）
            {
                s_ok = (s1 <= 0 && s1 >= denom);
                t_ok = (t1 <= 0 && t1 >= denom);
            }

            if (s_ok && t_ok)
            {
                // ★ 只有真正相交时才执行除法计算交点坐标
                double s = s1 / denom;
                Rintersection[i] = true;
                Rarrx[i] = Asx + s * Asex;   // 交点 x = A_sx + A_sex·s
                Rarry[i] = Asy + s * Asey;   // 交点 y = A_sy + A_sey·s
            }
            else
            {
                Rintersection[i] = false;
                Rarrx[i] = 0; Rarry[i] = 0;
            }
        }
    }
}

AVX2 版本（4-wide SIMD，含无除法判定优化）

AVX2 使用 256-bit 的 YMM 寄存器，每个寄存器可容纳 4 个双精度浮点数（__m256d）。循环步长由标量的 i++ 变为 i += 4，每次迭代同时处理 4 条线段的求交。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
// =====================================================================
// AVX2 矢量化线段求交 (一次处理 4 条线段)
//
// AVX2 关键 Intrinsics 说明:
//   _mm256_set1_pd(x)       将标量 x 广播到 4 个 lane
//   _mm256_load_pd(ptr)     从对齐内存加载 4 个 double
//   _mm256_store_pd(ptr, v)  存储 4 个 double 到对齐内存
//   _mm256_add_pd / _mm256_sub_pd / _mm256_mul_pd / _mm256_div_pd
//                            逐元素加/减/乘/除
//   _mm256_cmp_pd(a, b, imm) 逐元素比较，返回掩码向量（全1或全0）
//   _mm256_and_pd / _mm256_or_pd  逐元素按位与/或
//   _mm256_movemask_pd(v)    提取每个 lane 的最高位组成 4-bit 整数
// =====================================================================

// 广播常量：全 0 和全 1 的向量，供比较指令使用
const __m256d zero256 = _mm256_set1_pd(0);
const __m256d one256  = _mm256_set1_pd(1);

// 矢量化 contain：一次判定 4 个点是否分别在 4 条线段上
__m256d containSIMD(
    __m256d vAsx, __m256d vAsy,    // 4 条线段的起点（广播自同一条线段 A）
    __m256d vAsex, __m256d vAsey,  // 4 条线段的方向向量
    __m256d vcx, __m256d vcy)      // 4 个待判定点
{
    // 从 x 和 y 坐标分别反解参数 t
    __m256d tx = _mm256_div_pd(          // 从 x 坐标: t_x = (c_x − A_sx) / A_sex
        _mm256_sub_pd(vcx, vAsx), vAsex);
    __m256d ty = _mm256_div_pd(          // 从 y 坐标: t_y = (c_y − A_sy) / A_sey
        _mm256_sub_pd(vcy, vAsy), vAsey);

    // 三条判定条件逐元素求与: t_x == t_y && t_x >= 0 && t_x <= 1
    __m256d c1 = _mm256_cmp_pd(tx, ty,     _CMP_EQ_OS);  // tx == ty
    __m256d c2 = _mm256_cmp_pd(tx, zero256, _CMP_GE_OS); // tx >= 0
    __m256d c3 = _mm256_cmp_pd(tx, one256,  _CMP_LE_OS); // tx <= 1

    return _mm256_and_pd(c1, _mm256_and_pd(c2, c3));
}

void intersectSIMD(
    double Asx,           // 线段 A 起点 x
    double Asy,           // 线段 A 起点 y
    double Aex,           // 线段 A 终点 x
    double Aey,           // 线段 A 终点 y
    double *Barrsx,       // 线段集 B 起点 x 数组（32 字节对齐）
    double *Barrsy,       // 线段集 B 起点 y 数组（32 字节对齐）
    double *Barrex,       // 线段集 B 终点 x 数组（32 字节对齐）
    double *Barrey,       // 线段集 B 终点 y 数组（32 字节对齐）
    double *Rarrx,        // [输出] 交点 x 坐标数组（32 字节对齐）
    double *Rarry,        // [输出] 交点 y 坐标数组（32 字节对齐）
    bool *Rintersection,  // [输出] 是否直线相交
    bool *Roverlap,       // [输出] 是否共线重叠
    int size)             // 线段集 B 的大小（应为 4 的倍数）
{
    // ---- 线段 A 广播到 4 个 lane ----
    __m256d vAsx = _mm256_set1_pd(Asx);
    __m256d vAsy = _mm256_set1_pd(Asy);
    __m256d vAex = _mm256_set1_pd(Aex);
    __m256d vAey = _mm256_set1_pd(Aey);

    // 线段 A 的方向向量（各分量在 4 个 lane 中相同）
    __m256d vAsex = _mm256_sub_pd(vAex, vAsx);  // A_sex = A_ex − A_sx
    __m256d vAsey = _mm256_sub_pd(vAey, vAsy);  // A_sey = A_ey − A_sy

    // 每次迭代处理 4 条 B 线段
    for (int i = 0; i < size; i += 4)
    {
        // ---- 加载 4 条 B 线段的起终点 ----
        __m256d vBsx = _mm256_load_pd(Barrsx + i);
        __m256d vBsy = _mm256_load_pd(Barrsy + i);
        __m256d vBex = _mm256_load_pd(Barrex + i);
        __m256d vBey = _mm256_load_pd(Barrey + i);

        // ---- 计算 B 线段的方向向量 ----
        __m256d vBsex = _mm256_sub_pd(vBex, vBsx);
        __m256d vBsey = _mm256_sub_pd(vBey, vBsy);

        // ---- 行列式 Δ = A_sey·B_sex − A_sex·B_sey ----
        __m256d denom = _mm256_sub_pd(
            _mm256_mul_pd(vAsey, vBsex),
            _mm256_mul_pd(vAsex, vBsey));

        // ================================================================
        // 分支 1: Δ == 0 → 平行/共线，判定 overlap
        // ================================================================
        // 比较得到掩码向量（Δ==0 的 lane 为全 1，否则全 0）
        __m256d denom_eq_zero = _mm256_cmp_pd(denom, zero256, _CMP_EQ_OS);

        // movemask: 提取每个 lane 的高位组成 4-bit 整数，
        // 若非零则说明至少有一个 lane 的 Δ == 0
        int denom_zero_mask = _mm256_movemask_pd(denom_eq_zero);

        __m256d overlap;
        if (denom_zero_mask != 0)
        {
            // ★ 至少有一条线段与 A 平行，需要执行重叠判定
            // 检查 4 条 B 线段的 4 个端点是否落在 A 上，或 A 的端点落在 B 上
            overlap = _mm256_and_pd(
                denom_eq_zero,  // 只保留 Δ==0 的 lane
                _mm256_or_pd(   // 四种端点落在线段上的情况求或
                    _mm256_or_pd(
                        // B起点是否在A上
                        containSIMD(vAsx, vAsy, vAsex, vAsey, vBsx, vBsy),
                        // B终点是否在A上
                        containSIMD(vAsx, vAsy, vAsex, vAsey, vBex, vBey)),
                    _mm256_or_pd(
                        // A起点是否在B上（交换角色）
                        containSIMD(vBsx, vBsy, vBsex, vBsey, vAsx, vAsy),
                        // A终点是否在B上
                        containSIMD(vBsx, vBsy, vBsex, vBsey, vAex, vAey))));
        }
        else
        {
            // ★ 没有平行线段，完全跳过重叠判定（常见路径，性能关键）
            overlap = zero256;
        }

        // ================================================================
        // 分支 2: Δ != 0 → 不平行，判定直线相交
        // ================================================================
        __m256d denom_ne_zero = _mm256_cmp_pd(denom, zero256, _CMP_NEQ_OS);
        int denom_ne_mask = _mm256_movemask_pd(denom_ne_zero);

        __m256d intersect;
        if (denom_ne_mask != 0)
        {
            // ---- 计算参数 s 的分子和 t 的分子 ----
            // s₁ = A_sx·B_sey + B_sx·(A_sy−B_ey) + B_ex·(B_sy−A_sy)
            __m256d s1 = _mm256_add_pd(
                _mm256_add_pd(
                    _mm256_mul_pd(vAsx, vBsey),
                    _mm256_mul_pd(vBsx, _mm256_sub_pd(vAsy, vBey))),
                _mm256_mul_pd(vBex, _mm256_sub_pd(vBsy, vAsy)));

            // t₁ = A_sx·(A_ey−B_sy) + A_ex·(B_sy−A_sy) − B_sx·A_sey
            __m256d t1 = _mm256_sub_pd(
                _mm256_add_pd(
                    _mm256_mul_pd(vAsx, _mm256_sub_pd(vAey, vBsy)),
                    _mm256_mul_pd(vAex, _mm256_sub_pd(vBsy, vAsy))),
                _mm256_mul_pd(vBsx, vAsey));

            // ---- ★ 无除法判定 s ∈ [0,1] ----
            // 根据 Δ 的正负号分别处理
            __m256d denom_gt = _mm256_cmp_pd(denom, zero256, _CMP_GT_OS);
            __m256d denom_lt = _mm256_cmp_pd(denom, zero256, _CMP_LT_OS);

            // Δ > 0 时: 0 <= s₁ <= Δ
            // Δ < 0 时: Δ <= s₁ <= 0
            __m256d s_ok = _mm256_or_pd(
                _mm256_and_pd(denom_gt,
                    _mm256_and_pd(
                        _mm256_cmp_pd(s1, zero256, _CMP_GE_OS),  // s₁ >= 0
                        _mm256_cmp_pd(s1, denom,   _CMP_LE_OS))), // s₁ <= Δ
                _mm256_and_pd(denom_lt,
                    _mm256_and_pd(
                        _mm256_cmp_pd(s1, zero256, _CMP_LE_OS),  // s₁ <= 0
                        _mm256_cmp_pd(s1, denom,   _CMP_GE_OS))));// s₁ >= Δ

            // 同理判定 t ∈ [0,1]: Δ>0 → 0<=t₁<=Δ; Δ<0 → Δ<=t₁<=0
            __m256d t_ok = _mm256_or_pd(
                _mm256_and_pd(denom_gt,
                    _mm256_and_pd(
                        _mm256_cmp_pd(t1, zero256, _CMP_GE_OS),
                        _mm256_cmp_pd(t1, denom,   _CMP_LE_OS))),
                _mm256_and_pd(denom_lt,
                    _mm256_and_pd(
                        _mm256_cmp_pd(t1, zero256, _CMP_LE_OS),
                        _mm256_cmp_pd(t1, denom,   _CMP_GE_OS))));

            // 综合判定: Δ != 0 且 s_ok 且 t_ok
            intersect = _mm256_and_pd(
                denom_ne_zero, _mm256_and_pd(s_ok, t_ok));

            // ★ 检查是否有 lane 真正相交
            int intersect_mask = _mm256_movemask_pd(intersect);
            if (intersect_mask != 0)
            {
                // 至少有一个交点，计算 s 和交点坐标（此处才执行除法）
                __m256d s   = _mm256_div_pd(s1, denom);
                __m256d vRx = _mm256_add_pd(vAsx, _mm256_mul_pd(s, vAsex));
                __m256d vRy = _mm256_add_pd(vAsy, _mm256_mul_pd(s, vAsey));
                _mm256_store_pd(Rarrx + i, vRx);
                _mm256_store_pd(Rarry + i, vRy);
            }
            else
            {
                // 4 个 lane 都不相交，无需除法，直接存储零
                _mm256_store_pd(Rarrx + i, zero256);
                _mm256_store_pd(Rarry + i, zero256);
            }
        }
        else
        {
            // 所有 4 条线段都平行于 A，不可能有普通交点
            intersect = zero256;
            _mm256_store_pd(Rarrx + i, zero256);
            _mm256_store_pd(Rarry + i, zero256);
        }

        // ---- 从向量掩码中提取标量 bool 结果 ----
        for (int j = 0; j < 4; j++)
        {
            // AVX2 比较返回全 1（-1.0 ≈ NaN）或全 0（0.0）
            Rintersection[i + j] = (intersect[j] != 0.0);
            Roverlap[i + j]      = (overlap[j]   != 0.0);
        }
    }
}

AVX-512F 版本（8-wide SIMD，含无除法判定优化）

AVX-512F 是 Intel 在 2013 年引入的 512-bit SIMD 指令集。相比 AVX2，它的关键改进包括：

  • 512-bit ZMM 寄存器：每寄存器容纳 8 个 double（__m512d）
  • 掩码寄存器 __mmask8：比较指令直接返回 8-bit 掩码，用位运算组合，比向量掩码更高效
  • 32 个向量寄存器（AVX2 仅 16 个），减少寄存器溢出压力

注意：AMD Zen4 架构通过"双泵"（double-pumping）两个 256-bit 单元实现 AVX-512F，吞吐量与 AVX2 相近。但掩码优化和循环次数减半仍带来约 19% 的额外提升。Intel 原生 512-bit 实现可获得更高收益。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
const __m512d zero512 = _mm512_set1_pd(0);
const __m512d one512  = _mm512_set1_pd(1);

// AVX-512F 版 contain：返回 __mmask8 掩码而非向量
// __mmask8 是 8-bit 整数，每个 bit 对应一个 lane 的比较结果
__mmask8 containAVX512F(
    __m512d vAsx, __m512d vAsy,
    __m512d vAsex, __m512d vAsey,
    __m512d vcx,  __m512d vcy)
{
    __m512d tx = _mm512_div_pd(_mm512_sub_pd(vcx, vAsx), vAsex);
    __m512d ty = _mm512_div_pd(_mm512_sub_pd(vcy, vAsy), vAsey);

    // _mm512_cmp_pd_mask: 比较返回 __mmask8，三条条件用 & 组合
    return _mm512_cmp_pd_mask(tx, ty,     _CMP_EQ_OS) &
           _mm512_cmp_pd_mask(tx, zero512, _CMP_GE_OS) &
           _mm512_cmp_pd_mask(tx, one512,  _CMP_LE_OS);
}

void intersectAVX512F(
    double Asx, double Asy, double Aex, double Aey,
    double *Barrsx, double *Barrsy, double *Barrex, double *Barrey,
    double *Rarrx, double *Rarry,
    bool *Rintersection, bool *Roverlap,
    int size)
{
    // 线段 A 广播到 8 个 lane
    __m512d vAsx = _mm512_set1_pd(Asx);
    __m512d vAsy = _mm512_set1_pd(Asy);
    __m512d vAex = _mm512_set1_pd(Aex);
    __m512d vAey = _mm512_set1_pd(Aey);
    __m512d vAsex = _mm512_sub_pd(vAex, vAsx);
    __m512d vAsey = _mm512_sub_pd(vAey, vAsy);

    // 每次迭代处理 8 条 B 线段
    for (int i = 0; i < size; i += 8)
    {
        // 加载 8 条 B 线段的起终点（64 字节对齐）
        __m512d vBsx = _mm512_load_pd(Barrsx + i);
        __m512d vBsy = _mm512_load_pd(Barrsy + i);
        __m512d vBex = _mm512_load_pd(Barrex + i);
        __m512d vBey = _mm512_load_pd(Barrey + i);
        __m512d vBsex = _mm512_sub_pd(vBex, vBsx);
        __m512d vBsey = _mm512_sub_pd(vBey, vBsy);

        // 行列式
        __m512d denom = _mm512_sub_pd(
            _mm512_mul_pd(vAsey, vBsex),
            _mm512_mul_pd(vAsex, vBsey));

        // ---- 平行/共线分支 ----
        __mmask8 denom_zero = _mm512_cmp_pd_mask(denom, zero512, _CMP_EQ_OS);
        __mmask8 overlap;
        if (denom_zero)
        {
            // 至少有一条 B 线段与 A 平行
            overlap = denom_zero &
                      (containAVX512F(vAsx, vAsy, vAsex, vAsey, vBsx, vBsy) |
                       containAVX512F(vAsx, vAsy, vAsex, vAsey, vBex, vBey) |
                       containAVX512F(vBsx, vBsy, vBsex, vBsey, vAsx, vAsy) |
                       containAVX512F(vBsx, vBsy, vBsex, vBsey, vAex, vAey));
        }
        else
        {
            overlap = 0;  // 没有平行线段，跳过重叠计算（常见路径）
        }

        // ---- 不平行分支：无除法判定相交 ----
        __mmask8 denom_ne = _mm512_cmp_pd_mask(denom, zero512, _CMP_NEQ_OS);
        __mmask8 intersect;
        if (denom_ne)
        {
            // 计算 s 和 t 的分子
            __m512d s1 = _mm512_add_pd(
                _mm512_add_pd(
                    _mm512_mul_pd(vAsx, vBsey),
                    _mm512_mul_pd(vBsx, _mm512_sub_pd(vAsy, vBey))),
                _mm512_mul_pd(vBex, _mm512_sub_pd(vBsy, vAsy)));

            __m512d t1 = _mm512_sub_pd(
                _mm512_add_pd(
                    _mm512_mul_pd(vAsx, _mm512_sub_pd(vAey, vBsy)),
                    _mm512_mul_pd(vAex, _mm512_sub_pd(vBsy, vAsy))),
                _mm512_mul_pd(vBsx, vAsey));

            // 根据 Δ 正负号做无除法判定
            __mmask8 denom_gt = _mm512_cmp_pd_mask(denom, zero512, _CMP_GT_OS);
            __mmask8 denom_lt = _mm512_cmp_pd_mask(denom, zero512, _CMP_LT_OS);

            // Δ>0: 0<=s₁<=Δ; Δ<0: Δ<=s₁<=0
            __mmask8 s_ok =
                (denom_gt &
                 _mm512_cmp_pd_mask(s1, zero512, _CMP_GE_OS) &
                 _mm512_cmp_pd_mask(s1, denom,   _CMP_LE_OS)) |
                (denom_lt &
                 _mm512_cmp_pd_mask(s1, zero512, _CMP_LE_OS) &
                 _mm512_cmp_pd_mask(s1, denom,   _CMP_GE_OS));

            // 同样判定 t ∈ [0,1]
            __mmask8 t_ok =
                (denom_gt &
                 _mm512_cmp_pd_mask(t1, zero512, _CMP_GE_OS) &
                 _mm512_cmp_pd_mask(t1, denom,   _CMP_LE_OS)) |
                (denom_lt &
                 _mm512_cmp_pd_mask(t1, zero512, _CMP_LE_OS) &
                 _mm512_cmp_pd_mask(t1, denom,   _CMP_GE_OS));

            intersect = denom_ne & s_ok & t_ok;

            if (intersect)
            {
                // 至少有一个交点，计算除法得到精确坐标
                __m512d s   = _mm512_div_pd(s1, denom);
                __m512d vRx = _mm512_add_pd(vAsx, _mm512_mul_pd(s, vAsex));
                __m512d vRy = _mm512_add_pd(vAsy, _mm512_mul_pd(s, vAsey));
                _mm512_store_pd(Rarrx + i, vRx);
                _mm512_store_pd(Rarry + i, vRy);
            }
            else
            {
                _mm512_store_pd(Rarrx + i, zero512);
                _mm512_store_pd(Rarry + i, zero512);
            }
        }
        else
        {
            intersect = 0;
            _mm512_store_pd(Rarrx + i, zero512);
            _mm512_store_pd(Rarry + i, zero512);
        }

        // 从 __mmask8 掩码中提取 bool 结果
        for (int j = 0; j < 8; j++)
        {
            Rintersection[i + j] = (intersect >> j) & 1;
            Roverlap[i + j]      = (overlap   >> j) & 1;
        }
    }
}

测试框架

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
// ===========================
// 功能正确性测试
// ===========================
void test()
{
    const int size = 7, pad8 = 8;  // pad8 用于 SIMD 对齐（避免越界访问）
    double Asx = 1, Asy = 1, Aex = 3, Aey = 3;

    // 7 个测试用例，覆盖: 普通相交, 不相交, 共线重叠
    // (1,3)-(3,1) x (1,1)-(3,3) = (2,2)
    // (1,3)-(3,2) x (1,1)-(3,3) = (2.33,2.33)
    // (1,2)-(2,1) x (1,1)-(3,3) = (1.5,1.5)
    // (1,2)-(3,10) x (1,1)-(3,3) = 不相交
    // (1,2)-(5,8) x (1,1)-(3,3)  = 不相交
    // (1,1)-(3,3) x (1,1)-(3,3)  = 重叠
    // (0,0)-(2,2) x (1,1)-(3,3)  = 重叠
    double varrsxs[8] = {1, 1, 1, 1, 1, 1, 0, 0};
    double varrsys[8] = {3, 3, 2, 2, 2, 1, 0, 0};
    double varrexs[8] = {3, 3, 2, 3, 5, 3, 2, 0};
    double varreys[8] = {1, 2, 1, 10, 8, 3, 2, 0};

    // 使用 _mm_malloc 分配 64 字节对齐内存（满足 AVX-512F 要求）
    double *Barrsx = (double *)_mm_malloc(pad8 * sizeof(double), 64);
    double *Barrsy = (double *)_mm_malloc(pad8 * sizeof(double), 64);
    double *Barrex = (double *)_mm_malloc(pad8 * sizeof(double), 64);
    double *Barrey = (double *)_mm_malloc(pad8 * sizeof(double), 64);
    memcpy(Barrsx, varrsxs, size * sizeof(double));
    memcpy(Barrsy, varrsys, size * sizeof(double));
    memcpy(Barrex, varrexs, size * sizeof(double));
    memcpy(Barrey, varreys, size * sizeof(double));

    double *Rx = (double *)_mm_malloc(pad8 * sizeof(double), 64);
    double *Ry = (double *)_mm_malloc(pad8 * sizeof(double), 64);
    bool Rintersection[size], Roverlap[size];

    // 依次用三种方法测试，验证结果一致性
    auto printResults = [&](const char *tag) {
        for (int i = 0; i < size; i++) {
            if (Rintersection[i])
                printf("%s intersect at: (%.2f,%.2f)\n", tag, Rx[i], Ry[i]);
            else {
                printf("no intersection");
                if (Roverlap[i]) printf(" but overlap");
                printf("\n");
            }
        }
        printf("\n");
    };

    intersectVector(Asx, Asy, Aex, Aey, Barrsx, Barrsy, Barrex, Barrey,
                    Rx, Ry, Rintersection, Roverlap, size);
    printResults("seq");

    intersectSIMD(Asx, Asy, Aex, Aey, Barrsx, Barrsy, Barrex, Barrey,
                  Rx, Ry, Rintersection, Roverlap, size);
    printResults("avx2");

    intersectAVX512F(Asx, Asy, Aex, Aey, Barrsx, Barrsy, Barrex, Barrey,
                     Rx, Ry, Rintersection, Roverlap, size);
    printResults("avx512");

    _mm_free(Barrsx); _mm_free(Barrsy); _mm_free(Barrex); _mm_free(Barrey);
    _mm_free(Rx); _mm_free(Ry);
}

// ===========================
// 内存批量性能测试
// 场景: 一条线段对 40,000 条线段求交，循环 50 次
// 总计 200 万次求交
// ===========================
void performanceTest()
{
    const int size = 8 * 5000;   // 40,000（确保被 8 整除）
    const int count = 50;        // 迭代次数

    double sx = 0, sy = 1, ex = 10, ey = 2;  // 线段 A

    // 分配对齐内存
    double *vsx = (double *)_mm_malloc(size * sizeof(double), 64);
    double *vsy = (double *)_mm_malloc(size * sizeof(double), 64);
    double *vex = (double *)_mm_malloc(size * sizeof(double), 64);
    double *vey = (double *)_mm_malloc(size * sizeof(double), 64);

    // 构造测试数据: 所有 B 线段与 A 平行（Δ == 0）
    // 这用于测试最坏情况（必须执行 overlap 计算）
    for (int i = 0; i < size; i++) {
        vsx[i] = -1; vsy[i] = 0;
        vex[i] = 10 * ((double)i / (double)size);
        vey[i] = 3;
    }

    double *vpix = (double *)_mm_malloc(size * sizeof(double), 64);
    double *vpiy = (double *)_mm_malloc(size * sizeof(double), 64);
    bool *intersectSuccess = new bool[size];
    bool *overlap = new bool[size];

    clock_t t0;
    double t_seq, t_avx2, t_avx512;

    // Sequential
    t0 = clock();
    for (int i = 0; i < count; i++)
        intersectVector(sx, sy, ex, ey, vsx, vsy, vsx, vsy,
                        vpix, vpiy, intersectSuccess, overlap, size);
    t_seq = (double)(clock() - t0) / CLOCKS_PER_SEC;
    printf("  [Sequential ]  %f s\n", t_seq);

    // AVX2
    t0 = clock();
    for (int i = 0; i < count; i++)
        intersectSIMD(sx, sy, ex, ey, vsx, vsy, vsx, vsy,
                      vpix, vpiy, intersectSuccess, overlap, size);
    t_avx2 = (double)(clock() - t0) / CLOCKS_PER_SEC;
    printf("  [AVX2       ]  %f s  (加速 %.1f%% vs seq)\n",
           t_avx2, (t_seq - t_avx2) / t_seq * 100);

    // AVX-512F
    t0 = clock();
    for (int i = 0; i < count; i++)
        intersectAVX512F(sx, sy, ex, ey, vsx, vsy, vsx, vsy,
                         vpix, vpiy, intersectSuccess, overlap, size);
    t_avx512 = (double)(clock() - t0) / CLOCKS_PER_SEC;
    printf("  [AVX-512F   ]  %f s  (加速 %.1f%% vs seq, %.1f%% vs avx2)\n",
           t_avx512, (t_seq - t_avx512) / t_seq * 100,
           (t_avx2 - t_avx512) / t_avx2 * 100);

    _mm_free(vsx); _mm_free(vsy); _mm_free(vex); _mm_free(vey);
    _mm_free(vpix); _mm_free(vpiy);
    delete[] intersectSuccess; delete[] overlap;
}

// ===========================
// 文件 I/O: 读取/写入多边形顶点
// ===========================
// 文件格式:
//   第一行: 顶点数量 N
//   后续 N 行: x,y (每行一个顶点)
void read_vertices(const char *filename, int &size,
                   double *&x_list, double *&y_list)
{
    FILE *fp = fopen(filename, "r");
    if (!fp) { printf("Failed to open %s\n", filename); exit(1); }
    fscanf(fp, "%d", &size);

    // 分配对齐内存（pad 到 8 的倍数，确保 SIMD 不会越界）
    int size_pad = ((size + 7) / 8) * 8;
    x_list = (double *)_mm_malloc(size_pad * sizeof(double), 64);
    y_list = (double *)_mm_malloc(size_pad * sizeof(double), 64);

    for (int i = 0; i < size; i++)
        fscanf(fp, "%lf,%lf", x_list + i, y_list + i);
    printf("%s: read %d vertices\n", filename, size);
    fclose(fp);
}

void write_vertices(const char *filename,
                    vector<pair<double, double>> &vec)
{
    FILE *fp = fopen(filename, "w");
    fprintf(fp, "%d\n", (int)vec.size());
    for (size_t i = 0; i < vec.size(); i++)
        fprintf(fp, "%f,%f\n", vec[i].first, vec[i].second);
    fclose(fp);
}

// ===========================
// 多边形线段求交基准测试
//
// 读取两个多边形的顶点文件，将多边形边作为线段集：
//   source: 顶点数较少的文件 (inputLess) 的每条边 = 线段 A
//   target: 顶点数较多的文件 (inputMore) 的每条边 = 线段集 B
//
// 对于每个 source 线段，调用 method 与所有 target 线段求交，
// 收集全部交点。
// ===========================
using IntersectMethod = void (*)(
    double, double, double, double,
    double*, double*, double*, double*,
    double*, double*, bool*, bool*, int);

double polyTest(
    const char *inputMore,   // 顶点数较多的多边形文件
    const char *inputLess,   // 顶点数较少的多边形文件
    const char *output,      // 交点输出文件
    IntersectMethod method)  // 求交算法（sequential/avx2/avx512f）
{
    clock_t t0 = clock();

    // 读取两个多边形的顶点
    int s_size, c_size;
    double *s_x, *s_y, *c_x, *c_y;
    read_vertices(inputMore, s_size, s_x, s_y);
    read_vertices(inputLess, c_size, c_x, c_y);

    int numTgt = s_size - 1;                    // 目标线段数
    int numTgt_pad = ((numTgt + 7) / 8) * 8;    // SIMD 对齐

    // 目标线段集 B: 将相邻顶点对组成线段
    double *Bsx = s_x, *Bsy = s_y;               // B_i 起点 = 顶点 i
    double *Bex = (double *)_mm_malloc(numTgt_pad * sizeof(double), 64);
    double *Bey = (double *)_mm_malloc(numTgt_pad * sizeof(double), 64);
    memcpy(Bex, s_x + 1, numTgt * sizeof(double));  // B_i 终点 = 顶点 i+1
    memcpy(Bey, s_y + 1, numTgt * sizeof(double));

    vector<pair<double, double>> intersections;

    // 临时数组（每次调用复用）
    double *vpix = (double *)_mm_malloc(numTgt_pad * sizeof(double), 64);
    double *vpiy = (double *)_mm_malloc(numTgt_pad * sizeof(double), 64);
    bool *intersectSuccess = new bool[numTgt_pad];
    bool *overlap = new bool[numTgt_pad];

    // 遍历 source 的每条边（线段 A）
    for (int i = 0; i < c_size - 1; i++)
    {
        // 调用求交算法: source 第 i 条边 vs 全部 target 边
        method(c_x[i], c_y[i], c_x[i + 1], c_y[i + 1],
               Bsx, Bsy, Bex, Bey,
               vpix, vpiy, intersectSuccess, overlap, numTgt);

        // 收集交点
        for (int j = 0; j < numTgt; j++)
            if (intersectSuccess[j])
                intersections.push_back(make_pair(vpix[j], vpiy[j]));
    }

    // 释放资源
    _mm_free(Bex); _mm_free(Bey);
    _mm_free(vpix); _mm_free(vpiy);
    delete[] intersectSuccess; delete[] overlap;

    write_vertices(output, intersections);

    double elapsed = (double)(clock() - t0) / CLOCKS_PER_SEC;
    printf("  耗时: %f s,  交点数量: %d\n",
           elapsed, (int)intersections.size());
    return elapsed;
}

// ===========================
// 主函数
// ===========================
int main()
{
    cout << "====== 功能测试 ======" << endl;
    test();

    cout << "\n====== 内存批量测试 (40,000 线段 × 50 次) ======" << endl;
    performanceTest();

    const char *inputMore = "poly/s.txt";
    const char *inputLess = "poly/c.txt";

    cout << "\n====== 多边形求交测试 (~1 亿线段对) ======" << endl;

    cout << "--- Sequential ---" << endl;
    double t_seq = polyTest(inputMore, inputLess,
                            "poly/seq_output.txt", intersectVector);

    cout << "--- AVX2 ---" << endl;
    double t_avx2 = polyTest(inputMore, inputLess,
                             "poly/avx2_output.txt", intersectSIMD);

    cout << "--- AVX-512F ---" << endl;
    double t_avx512 = polyTest(inputMore, inputLess,
                               "poly/avx512_output.txt", intersectAVX512F);

    printf("\n====== 最终加速效果 ======\n");
    printf("  Sequential:  %f s  (基线)\n", t_seq);
    printf("  AVX2:        %f s  (加速 %.1f%%)\n",
           t_avx2, (t_seq - t_avx2) / t_seq * 100);
    printf("  AVX-512F:    %f s  (加速 %.1f%% vs seq, %.1f%% vs avx2)\n",
           t_avx512, (t_seq - t_avx512) / t_seq * 100,
           (t_avx2 - t_avx512) / t_avx2 * 100);

    return 0;
}

数据生成脚本

使用 Python 生成两个"多边形"的顶点数据：

  • s.txt（较大，20000 顶点）：参数方程 r(θ)=100+50sin⁡(0.02θ)r(\theta) = 100 + 50\sin(0.02\theta) 的螺旋线，确保与 c.txt 产生丰富交点
  • c.txt（较小，5000 顶点）：锯齿形折线 y=80sin⁡(0.15x)cos⁡(0.04x)y = 80\sin(0.15x)\cos(0.04x)，从左到右横向穿过螺旋线
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
import math
import random
random.seed(42)

# s.txt: 螺旋多边形（20000 个顶点）
n_s = 20000
with open("poly/s.txt", "w") as f:
    f.write(f"{n_s}\n")
    for i in range(n_s):
        angle = i * 0.0005 * math.pi     # 角度逐步增大
        r = 100 + 50 * math.sin(i * 0.02)  # 半径周期性摆动
        x = r * math.cos(angle)
        y = r * math.sin(angle)
        f.write(f"{x:.6f},{y:.6f}\n")

# c.txt: 锯齿折线横穿（5000 个顶点）
n_c = 5000
with open("poly/c.txt", "w") as f:
    f.write(f"{n_c}\n")
    for i in range(n_c):
        x = -150 + 300 * i / (n_c - 1)      # 从左到右均匀分布
        y = 80 * math.sin(i * 0.15) * math.cos(i * 0.04)  # 波形振荡
        f.write(f"{x:.6f},{y:.6f}\n")

print(f"s.txt: {n_s} 个顶点")
print(f"c.txt: {n_c} 个顶点")
print(f"总线段对: {(n_c-1)} × {(n_s-1)} = {(n_c-1)*(n_s-1):,}")

将脚本保存为 gen_data.py，运行前先创建目录：

1
2
mkdir -p poly
python3 gen_data.py

编译与运行

编译环境要求

  • 编译器：GCC 9+ 或 Clang 10+（需支持 AVX2 和 AVX-512F intrinsics）
  • CPU：Intel Haswell+ (AVX2) / Skylake-X+ (AVX-512F)，或 AMD Zen4+ (AVX2 + AVX-512F)
  • 操作系统：Linux / WSL2

编译命令

1
2
3
4
5
6
7
8
# 编译（需同时开启 AVX2 和 AVX-512F）
g++ -std=c++11 -O2 -mavx2 -mavx512f main.cpp -o main

# 或使用 Clang
clang++ -std=c++11 -O2 -mavx2 -mavx512f main.cpp -o main

# 若仅支持 AVX2（编译含 AVX-512F 的代码但不执行）
g++ -std=c++11 -O2 -mavx2 main.cpp -o main_avx2_only

运行

1
2
3
4
5
# 完整测试（功能验证 + 性能基准）
./main

# 仅功能测试
./main 2>&1 | head -30

性能结果与分析

测试环境

项目 配置
CPU AMD Ryzen 7 7840HS (Zen4, 8C/16T, 3.8GHz 基频)
L1 Cache 256KB (data) + 256KB (instruction) per core
L2 Cache 1MB per core (8MB total)
L3 Cache 16MB (shared)
内存 LPDDR5-6400
环境 WSL2 (Hyper-V), GCC 11.4

注意：AMD Zen4 的 AVX-512F 通过双泵（double-pumping）两个 256-bit 单元实现，吞吐量 ≈\approx AVX2。在 Intel 原生 512-bit CPU 上，AVX-512F 的加速比预期更高。

测试结果

1. 内存批量测试（40,000 线段 ×\times 50 次 == 200 万次求交）

在此测试中，所有 BB 线段与 AA 平行（Δ≡0\Delta \equiv 0），必须执行 overlap 计算。数据完全在 L2 缓存内（~640KB）。

方法 耗时 加速比 (vs seq) 加速比 (vs avx2)
Sequential 0.0165s 1.00×1.00\times（基线） —
AVX2 (4-wide) 0.0045s 3.67×3.67\times —
AVX-512F (8-wide) 0.0041s 4.05×4.05\times 1.10×1.10\times

2. 多边形求交测试（19,999×4,999≈119{,}999 \times 4{,}999 \approx 1 亿对）

使用螺旋线与折线数据进行真实多边形求交。交点数量 3,780 个，相交率仅 0.0038%0.0038\%——99.996%99.996\% 的除法被无除法判定优化消除。

方法 耗时 加速比 (vs seq) 加速比 (vs avx2)
Sequential 0.2211s 1.00×1.00\times（基线） —
AVX2 (4-wide) 0.1914s 1.13×1.13\times —
AVX-512F (8-wide) 0.1550s 1.30×1.30\times 1.19×1.19\times

为什么加速比达不到 4×4\times / 8×8\times？

这是计算特征决定的硬件限制，与实现质量无关。

1. 计算密集度极低（内存瓶颈）

每次 SIMD 迭代加载 44 个数组（Bsx,Bsy,Bex,BeyB_{sx}, B_{sy}, B_{ex}, B_{ey}）== 1616 或 3232 个 double，运算仅 ∼15\sim 15 条向量指令。算术强度（FLOPS/Byte）极低，瓶颈在内存带宽而非计算单元。SIMD 增加计算宽度但无法加速数据加载。

2. 标量乱序执行有天然优势

各线段对之间完全独立，无数据依赖。现代 CPU 的乱序执行窗口（∼200\sim 200+ μ\muop）可以同时执行多个线段对的标量计算，部分重叠延迟。SIMD 将元素"捆绑"在一起，反而限制了调度灵活性。

3. 多边形测试的额外限制

  • 数据量 ∼1\sim 1 MB（L2 缓存边界），缓存未命中比内存测试多
  • 无除法判定消除了除法瓶颈，但增加的比较指令使总指令数变多
  • 99.996%99.996\% 不相交意味着大部分计算在"快速路径"（无除法）上，瓶颈转移到乘加和比较

4. AMD Zen4 的 AVX-512 实现

Zen4 用两个 256-bit 单元拼接执行 512-bit 指令，吞吐量与 AVX2 同。AVX-512F 的 19%19\% 提升主要来自：掩码寄存器（减少向量 ↔\leftrightarrow 标量转换）、循环次数减半（更少的分支和迭代开销）、32 个寄存器（减轻溢出压力）。

在 Intel 原生 512-bit CPU（如 Xeon Platinum）上，AVX-512F 预期可达 2×2\times 于 AVX2。

WSL2 影响

WSL2 环境下 clock() 计时存在显著波动（同条件下最大偏差 ∼29%\sim 29\%），原因是 Hyper-V 虚拟化下 vCPU 由 Windows 宿主导调度。对纯计算吞吐量影响较小（<5%< 5\%），但影响了计时精度。建议在裸机 Linux 上做更精确的 benchmark。

关键优化总结

优化项 原理 效果
无除法判定 Δ>0  ⟹  (s∈[0,1]  ⟺  0≤s1≤Δ)\Delta > 0 \implies \big(s \in [0,1] \iff 0 \leq s_1 \leq \Delta\big) 消除 99.996%99.996\% 的除法
重叠计算跳过 _mm256_movemask_pd 检查 Δ=0\Delta = 0，全非零则跳过 多边形数据中几乎永不触发
交点坐标延迟计算 只有 intersect == true 时才做除法和坐标计算 除法从 2/2/对 →∼7.6×10−5/\to \sim 7.6 \times 10^{-5}/对
64 字节对齐 _mm_malloc(..., 64) 满足 AVX-512F 对齐要求
数组填充到 88 的倍数 分配时 pad 到 88 的倍数 避免 SIMD 越界读取

总结

本文以线段组求交为切入点，展示了从纯标量 →\to AVX2 →\to AVX-512F 的 SIMD 优化全流程。两项关键优化——无除法相交判定和重叠计算跳过——将除法从每条线段对的必执行路径上消除，使得在真实多边形数据上的加速比从最初的负值（SIMD 反而更慢）扭转为 AVX2 +13%+13\%、AVX-512F +30%+30\%。

加速比未达 4×4\times / 8×8\times 的根本原因在于该算法的算术强度太低（数据加载开销与计算量相当），这是硬件限制，并非实现不足。在 Δ≡0\Delta \equiv 0 的"重计算"场景下，AVX2 可获得 3.67×3.67\times 加速，AVX-512F 可获得 4.05×4.05\times 加速，验证了 SIMD 对计算密集型代码的巨大价值。

如需进一步提升多边形求交的整体吞吐量，方向在于增大单次处理的数据量（将多个源线段一起 SIMD 处理）以提升计算密度，或使用空间索引（如 R-tree）减少需要求交的线段对数量。