1. 从数学公式到可部署算子的挑战

在前两篇文章中，我们分别介绍了：

- 如何用切比雪夫多项式和查找表逼近菲涅尔余弦积分$C(x)$ 在中小$x$段的值；

- 如何实现高性能的 `sin/cos` 函数来支持渐近段计算。

但要把这些技术组装成一个真正可用的 Ascend C 算子，还需要解决很多工程问题：

- 多段函数的统一处理：如何在一个内核中优雅地处理近零段、主区间段、渐近段，以及负数对称？

- 掩码驱动的执行流：如何在 SIMD/SIMT 架构上避免分支，让所有线程都保持高利用率？

- 缓冲区资源管理：本地内存有限，如何规划几十个中间变量的生命周期，避免爆内存？

- 精度与性能的平衡：各个段之间的衔接是否平滑？整体误差是否可控？

本文将详细介绍这些工程实践，以及在开发过程中遇到的一些关键决策点。

---

2. 整体分段策略的演进

2.1 最初的想法：全用渐近公式

如第一篇所述，最开始的想法很简单：既然菲涅尔积分有渐近展开，那全区间都用渐近公式不就行了？

C(x) \approx \frac{1}{2} + \frac{1}{\pi x}\,\sin\!\left(\frac{\pi}{2}x^2\right).

实现了一版之后，发现问题很明显：

- 在$x < 2$时，误差大到完全不可用（相对误差 > 10%）；

- 即使用更高阶的渐近展开（加入 $1/x^3, 1/x^5$ 项），在中小 $x$段仍然不够准。

2.2 转折点：李俊昌（LIGC）公式

后来接触到 LIGC 公式（李俊昌等人提出的有理逼近公式），这类公式通过引入修正因子，能够在更大的$x$范围内保持高精度。

经过实测，在 $x \ge 10.2$ 的区间内，LIGC 公式能够满足精度要求（绝对误差 $< 10^{-5}$）。

这个发现让整个方案有了明确的分水岭：

- 大 $x$ 段（$x \ge 10.2$）：用 LIGC 公式（或简化版的渐近公式）；

- 中小 $x$ 段（$0 < x < 10.2$）：需要另找方法。

2.3 切比雪夫方案填补中小 x 段

对于$x \in (0, 10.2)$，最终选择了「压缩 + 切比雪夫 + LUT」的方案（见第一篇文章）。这个方案的优点是：

- 在每个子区间上，误差均匀分布；

- 实现结构简单，适合向量化；

- 离线优化好系数后，运行时只需要查表 + 递推。

2.4 完整的分段策略

综合起来，最终的分段策略是：

| 区间 | 方法 | 实现函数 |

|------|------|----------|

| $\|x\| \le \epsilon$ | $C(x) \approx x$ | 直接 `Select` |

| $\epsilon < \|x\| < 10.2$ | 切比雪夫多项式 + LUT | `CalculateChebyshev` |

| $\|x\| \ge 10.2$ | LIGC 公式或渐近公式 | `CalculateAsymptotic` |

| $x < 0$| 利用奇函数性质 $C(-x) = -C(x)$ | 最后统一处理符号 |

其中 $\epsilon$ 通常取 $10^{-8}$ 左右，保证数值稳定性。

图1 三代方案的对比

---

3. 掩码驱动的无分支实现

3.1 掩码构造

首先对输入 $x$构造各种掩码：

```cpp

AscendC::LocalTensor<float> x_abs = B_inner1.Get<float>();

AscendC::Abs(x_abs, x, length);

// 负数掩码

AscendC::LocalTensor<bool> mask_neg = B_mask1.Get<bool>();

AscendC::Compare(mask_neg, x, 0.0f, CMPMODE::LT, length);

// 近零掩码

AscendC::LocalTensor<bool> mask_near_zero = B_mask2.Get<bool>();

AscendC::Compare(mask_near_zero, x_abs, EPS, CMPMODE::LE, length);

// 渐近段掩码

AscendC::LocalTensor<bool> mask_asym = B_mask3.Get<bool>();

AscendC::Compare(mask_asym, x_abs, BOUNDARY_19, CMPMODE::GE, length);

// 主区间掩码 = 非近零 且 非渐近

AscendC::LocalTensor<bool> mask_main = B_mask4.Get<bool>();

AscendC::Not(mask_main, mask_near_zero, length);

AscendC::LocalTensor<bool> mask_temp = B_mask_temp.Get<bool>();

AscendC::Not(mask_temp, mask_asym, length);

AscendC::And(mask_main, mask_main, mask_temp, length);

```

其中 `BOUNDARY_19 = 10.2f`，`EPS = 1e-8f`。

3.2 分段计算与结果合并

初始化结果为 0，然后根据掩码逐段计算：

```cpp

AscendC::LocalTensor<float> y = B_result1.Get<float>();

AscendC::Duplicate(y, 0.0f, length);

// 近零段：C(x) ≈ x

AscendC::Select(y, mask_near_zero, x, y, length);

// 主区间段：切比雪夫

AscendC::LocalTensor<float> y_main = B_inner2.Get<float>();

CalculateChebyshev(y_main, x_abs, length);

AscendC::Select(y, mask_main, y_main, y, length);

// 渐近段：LIGC 公式

AscendC::LocalTensor<float> y_asym = B_inner3.Get<float>();

CalculateAsymptotic(y_asym, x_abs, length);

AscendC::Select(y, mask_asym, y_asym, y, length);

// 处理负数

AscendC::LocalTensor<float> y_neg = B_inner4.Get<float>();

AscendC::Muls(y_neg, y, -1.0f, length);

AscendC::Select(y, mask_neg, y_neg, y, length);

```

这样整个逻辑流是：

1. 所有元素并行计算各自的掩码；

2. 根据掩码分别调用对应的计算函数；

3. 用 `Select` 把各段结果合并到最终输出。

3.3 掩码操作的性能考量

虽然掩码操作本身是向量化的，但仍然有一些优化空间：

- **掩码合并**：尽量减少 `Compare` 和逻辑运算的次数，比如用一次 `And` 代替多次嵌套 `Select`；

- **短路评估**：在某些情况下，可以通过统计掩码中的非零元素数量，提前跳过某些计算分支（当前实现中为了代码简洁性没有做这个优化）；

- **掩码复用**：`mask_temp` 在多处被复用，避免分配过多的布尔缓冲区。

在实际测试中，掩码操作的开销占总计算时间的比例不到 5%，主要瓶颈还是在切比雪夫递推和三角函数计算上。

图2 掩码机制分析

---

4. 缓冲区生命周期规划

4.1 问题：本地内存是稀缺资源

在 Ascend AI Core 上，每个核的本地内存（UB）是有限的。如果算子内部有几十个中间变量，每个都单独分配缓冲区，很容易爆内存，导致：

- Tiling 尺寸被迫变小，影响吞吐；

- 甚至无法编译通过（编译器报错「本地内存不足」）。

4.2 解决方案：缓冲区复用

核心思路是：不同函数中的临时变量，如果生命周期不重叠，可以共享同一块缓冲区。

在代码中，预先分配了一批「功能性缓冲区」：

```cpp

// 主计算缓冲区

AscendC::GlobalTensor<float> B_inner1..B_inner6;  // 用于切比雪夫递推

AscendC::GlobalTensor<float> B_tri1, B_tri2;      // 用于三角函数

AscendC::GlobalTensor<float> B_result1;           // 最终结果

// 三角函数专用

AscendC::GlobalTensor<float> B_c1..B_c7;          // 常数、角度、幂次

// 掩码缓冲区

AscendC::GlobalTensor<bool> B_mask1..B_mask4, B_mask_temp;

```

然后在每个函数中，明确规定「谁用哪个缓冲区」：

| 函数 | 使用的缓冲区 | 用途 |

|------|-------------|------|

| `CalculateChebyshev` | `B_inner1..6` | $t, y, T_0, T_1, T_n, \text{tmp}$|

| `SinWithReduction` | `B_c1..c7, B_tri1, B_tri2` | $\theta_0, \phi, \text{sign}, \phi^{2..7}$ |

| `CalculateAsymptotic` | `B_inner1..3, B_c1..c3` | 中间项、角度、临时变量 |

| `Calculate` (主函数) | `B_inner1..4, B_result1` | $x_{\text{abs}}, y_{\text{main}}, y_{\text{asym}}, y_{\text{neg}}$ |

4.3 生命周期图与冲突检测

在开发过程中，利用一个简化的「缓冲区生命周期图」来辅助规划：

```

时间线 →

────────────────────────────────────────────

Calculate():

  B_inner1 (x_abs)     [====================]

  B_inner2 (y_main)         [======]

  B_inner3 (y_asym)              [======]

  B_inner4 (y_neg)                    [===]

CalculateChebyshev():

  B_inner1 (t)              [=====]

  B_inner3 (T0)             [=====]

  B_inner4 (T1)             [=====]

  B_inner5 (Tn)             [=====]

  B_inner6 (tmp)            [=====]

SinWithReduction():

  B_c1 (theta0)                  [====]

  B_c2 (phi2)                    [====]

  B_tri1 (phi)                   [====]

  B_tri2 (sign)                  [====]

```

通过这个图可以清楚地看到：

- `B_inner1` 在主函数中存储 `x_abs`，生命周期贯穿始终，不能被复用；

- `B_inner3` 在 `CalculateChebyshev` 中作为 `T0`，在主函数中作为 `y_asym`，两者不会同时活跃，可以复用；

- `B_c1..c7` 专门服务于三角函数，与切比雪夫的缓冲区完全隔离。

4.4 调试经验：缓冲区冲突

在开发过程中，曾经多次遇到「缓冲区冲突」的 bug，表现为：

- 某个中间结果被意外覆盖，导致最终输出错误；

- 误差突然在某些特定$x$ 值附近变大。

排查方法是：

1. 在每个函数入口和出口加日志，打印关键变量的前几个元素；

2. 检查是否有两个变量在同一时间段内都需要有效值，但被分配到了同一块缓冲区；

3. 逐步调整分配方案，直到所有冲突消除。

最终在代码注释中也记录了这些调整（例如「修复 B_inner3 与 B_tri1 的冲突」）。

通过这种「按需分配、用完就复用」的策略，把总缓冲区数量控制在了 20 个左右（包括掩码），相比 naive 实现（每个变量一个缓冲区）节省了 50% 以上的本地内存。

---

5. 渐近公式的实现与分段衔接

5.1 LIGC 公式的数学形式

李俊昌公式的一种形式是：

C(x) \approx \frac{1}{2} + f(x) \cdot \frac{\sin\!\left(\frac{\pi}{2}x^2\right)}{\pi x},

其中 $f(x)$ 是一个修正因子，常见的形式包括：

f(x) = 1 - \alpha\,e^{-\beta(x - x_0)},

或者更复杂的有理函数。这个修正因子的作用是：

- 在 $x$接近渐近段起点（例如$x \approx 10$）时，通过$e^{-\beta(x-x_0)}$ 平滑过渡；

- 在 $x$ 很大时，$f(x) \to 1$，退化为标准的渐近公式。

5.2 AscendC 实现

在 `CalculateAsymptotic` 函数中，实现大致如下：

```cpp

void CalculateAsymptotic(

    AscendC::LocalTensor<float> result,

    AscendC::LocalTensor<float> x_abs,

    uint32_t length

) {

    // 1. 计算 θ = (π/2) * x^2

    AscendC::LocalTensor<float> theta = B_c1.Get<float>();

    AscendC::Mul(theta, x_abs, x_abs, length);          // x^2

    AscendC::Muls(theta, theta, PI / 2.0f, length);     // (π/2) * x^2

    // 2. 计算 sin(θ)

    AscendC::LocalTensor<float> sin_val = B_c2.Get<float>();

    SinWithReduction(sin_val, theta, length);

    // 3. 计算 1 / (πx)

    AscendC::LocalTensor<float> inv_pix = B_inner1.Get<float>();

    AscendC::Muls(inv_pix, x_abs, PI, length);          // πx

    AscendC::Reciprocal(inv_pix, inv_pix, length);      // 1 / (πx)

    // 4. 修正因子 f(x) = 1 - α*exp(-β(x - x0))

    // 当前版本中简化为 f(x) ≈ 1 （或者用更简单的多项式修正）

    AscendC::LocalTensor<float> correction = B_c3.Get<float>();

    AscendC::Duplicate(correction, 1.0f, length);

   

    // 可选：加入修正项

    // AscendC::Subs(tmp, x_abs, X0, length);           // x - x0

    // AscendC::Muls(tmp, tmp, -BETA, length);          // -β(x - x0)

    // AscendC::Exp(tmp, tmp, length);                  // exp(...)

    // AscendC::Muls(tmp, tmp, ALPHA, length);          // α*exp(...)

    // AscendC::Subs(correction, 1.0f, tmp, length);    // 1 - α*exp(...)

    // 5. 组合结果：0.5 + f(x) * sin(θ) / (πx)

    AscendC::Mul(result, sin_val, inv_pix, length);     // sin(θ) / (πx)

    AscendC::Mul(result, result, correction, length);   // f(x) * ...

    AscendC::Adds(result, result, 0.5f, length);        // 0.5 + ...

}

```

5.3 分段边界的选择与衔接

选择 $x = 10.2$ 作为渐近段起点并不是随意的，而是通过离线扫参数确定的：

1. 精度测试：在 $x \in [8, 15]$范围内密集采样，分别计算切比雪夫方案和 LIGC 方案的误差；

2. 交叉点寻找：找到一个 $x_0$，使得在 $x < x_0$时切比雪夫更准，$x > x_0$时 LIGC 更准；

3. 平滑性检查：在 $x = x_0$ 附近，检查两个方案的输出是否连续（误差跳变 $< 10^{-7}$）。

最终确定 $x_0 = 10.2$时，两个方案在该点附近的输出差异在$10^{-7}$ 量级，满足平滑衔接的要求。

在实现中，通过掩码 `mask_asym` 和 `mask_main` 的互斥性，保证每个元素只会被其中一个方案处理，不会出现「两个方案都算、然后插值」的情况（那样反而会引入额外误差）。

5.4 边界情况的特殊处理

在$x$非常大时（例如$x > 100$），$\sin(\frac{\pi}{2}x^2)$的角度会非常大，虽然 `SinWithReduction` 理论上能处理，但浮点精度损失会变得明显。

一个可选的优化是：在$x > x_{\max}$时，直接输出 $C(x) = 0.5$（因为振荡项的振幅已经衰减到机器精度以下）。在当前实现中，$x_{\max}$被隐式设为约 50，超出部分由 LIGC 公式自然处理。


6. 可扩展性与通用化思路

这套实现的核心思想——「分段 + 多项式逼近 + 查找表 + 掩码驱动」——不仅适用于菲涅尔积分，也可以推广到其它特殊函数。

例如，要在 Ascend C 上实现贝塞尔函数 $J_0(x)$，可以采用类似的流程：

1. **分段策略**：

   - 小$x$段：用泰勒展开；

   - 中等$x$段：用切比雪夫多项式；

   - 大$x$段：用渐近展开（涉及 $\sin$和 $\cos$，可以复用现有实现）。

2. **离线优化**：

   - 确定分段边界和每段的阶数；

   - 拟合切比雪夫系数；

   - 验证精度和衔接平滑性。

3. **运行时实现**：

   - 构造掩码判断每个元素属于哪个段；

   - 调用对应的计算函数；

   - 用 `Select` 合并结果。

这种「模块化 + 查表驱动」的架构，可以大幅降低实现新特殊函数的工作量。

---

7. 总结与展望

本文详细介绍了菲涅尔余弦积分算子在 Ascend C 上的工程实践，重点讲解了：

1. 整体架构：从最初的纯渐近公式，演进到「LIGC + 切比雪夫 + LUT」的混合方案；

2. 掩码驱动：通过布尔掩码和 `Select` 实现无分支的多段函数处理；

3. 缓冲区管理：通过生命周期规划和复用，把本地内存占用降低 50%；

4. 分段衔接：精心选择边界点，保证各段误差均衡且平滑过渡；

这三篇文章构成了一个完整的技术分享系列，从数学原理、算法设计到工程实现，全面展现了在 Ascend C 上实现复杂数值算子的完整过程。希望这些经验能为其他开发者提供参考。

文章一：[文章一：基于切比雪夫多项式与查找表的菲涅尔余弦积分逼近：从渐近公式到分段优化的演进 —— Ascend C实现与优化](文章一：基于切比雪夫多项式与查找表的菲涅尔余弦积分逼近：从渐近公式到分段优化的演进 —— Ascend C实现与优化-CSDN博客)

文章二：[文章二：面向 Ascend C 的高性能 sin/cos 实现：从角度规约到向量化多项式逼近](文章二：面向 Ascend C 的高性能 sin/cos 实现：从角度规约到向量化多项式逼近-CSDN博客)

---


 

# c语言 # 开发语言
Logo
昇腾开源生态专区

昇腾计算产业是基于昇腾系列（HUAWEI Ascend）处理器和基础软件构建的全栈 AI计算基础设施、行业应用及服务，https://devpress.csdn.net/organization/setting/general/146749包括昇腾系列处理器、系列硬件、CANN、AI计算框架、应用使能、开发工具链、管理运维工具、行业应用及服务等全产业链

更多推荐

cover

HarmonyOS 6 自定义人脸识别模型10：基于MindSpore Lite框架的自定义人脸识别功能实现

avatar 昇腾开源生态专区
cover

DeepSeek 崩了 13 小时，不是故障，是 V4 在换引擎

avatar 昇腾开源生态专区

体系结构论文（107）：AscendOptimizer: Episodic Agent for Ascend NPU Operator Optimization

摘要：本文提出AscendOptimizer系统，针对华为Ascend NPU的AscendC算子优化难题，通过双阶段交替优化方法实现性能提升。系统将算子拆分为host侧tiling program和device侧kernel program：Stage I采用进化搜索优化tiling策略，利用硬件反馈探索可行解空间；Stage II通过"优化回退"机制从优质kernel反向构

avatar 昇腾开源生态专区