DragonArray 文档

LoongArch / LSX 优化日志

这是 loongson-experimental 分支上优化工作的时间顺序记录。实测吞吐量见 BENCHMARKS.zh-CN.md。

待办事项

统一记录在这里，避免跨 session 遗失。完成的项目会移出此列表，并成为下方按日期 排列的独立章节。

• ~~f64 sin/cos 的 SIMD Payne-Hanek range reduction~~ — 已于 2026-05-24 完成，见下文。reduction 按 lane 以 scalar 执行（bit-table lookup 也是按 lane）； 随后的 polynomial 仍保持 vectorized。slow path 速度大致与 libm 持平。进一步 vectorize reduction 本身仍是后续工作，需要为 2/π bit table chunk 实现 SIMD gather。
• ~~cbrt f64 direct kernel~~ — 已于 2026-05-23 完成，见下文。
• ~~log1p / arctanh f64 direct polynomial~~ — 已于 2026-05-24 完成， 见下文。（arcsinh 被跳过，相关章节说明了原因。）
• 提交上游。 Payne-Hanek 完成后，精度可以表述为“在完整 f64 domain 上与 libm 一致”。届时 f64 sin/cos hybrid 可考虑提交 upstream PR，并改用 build flag，而 不是 __loongarch__ 进行 gating。

2026-05-17 — NPYV f32 exp / log

文件： numpy/_core/src/umath/loops_exponent_log.dispatch.c.src

改动

在 loops_exponent_log.dispatch.c.src 中加入 NPYV 风格的 simd_exp_FLOAT_npyv / simd_log_FLOAT_npyv。原文件中的 float32 exp/log 使用 raw x86 intrinsic（AVX2/AVX-512），因此在 LoongArch 上会回退到 scalar expf / logf。新 block 的 gate 为：

#if NPY_SIMD_F32 && !defined(SIMD_AVX2_FMA3) && !defined(SIMD_AVX512F)

因此它覆盖所有尚无 AVX 路径的 NPYV target，目前包括 LSX、NEON 和 VSX。 dispatcher entry 新增 #elif defined(NPYV_IMPL_F32_EXP_LOG) 分支：contiguous input 使用新 kernel，strided input 回退到 scalar libm。

实现方式

算法与同一文件中的 AVX2 实现一致：

• exp： Cody-Waite range reduction：y = x - k·ln(2)，其中 k = rint(x·log2e)；随后 exp(x) = (P(y)/Q(y)) · 2^k。P 为 degree-5，Q 为 degree-2。通过 npyv_shli_s32 与 npyv_reinterpret_f32_s32 把 (k+127) 写入 IEEE-754 exponent field，构造 2^k。
• log： 以 bit 操作提取 mantissa，并强制 exponent field 为 126，使 mantissa 落在 [0.5, 1)；随后用 m ≤ 1/√2 的 shift 让 mantissa 靠近 1，再计算 log(1 + (m−1)) ≈ P(m−1) / Q(m−1)，P 与 Q 均为 degree-5。

coefficient 原样复用 numpy/_core/src/umath/npy_simd_data.h 中的数据。

决策

• 暂时使用一个 LSX intrinsic，尚未完全 portable。 新 block 中唯一非 NPYV 操作是 int32→float32 的 __lsx_vffint_s_w，封装在 npyv__cvt_f32_s32 中，并带有 TODO。NPYV 尚未提供 portable 版本；若将其提升到 API，NEON / VSX 也能直接获得同一 exp/log 路径。
• 仅优化 contiguous fast path。 strided case 回退到 scalar libm。可以用 npyv_loadn / npyv_storen 实现 strided NPYV，但会为较少见的 shape 增加复杂度。
• 暂不加入 DOUBLE 路径。 LoongArch 上的 f64 仍使用 scalar libm。

开发期间修复的正确性问题

1. NaN / inf input 触发不应出现的 FP exception。 ordered LSX compare （vfcmp.clt.s 等）遇到 NaN 会触发 INVALID。两个 kernel 现在都会在任何 signaling compare 之前，用 non-signaling npyv_notnan_f32 把 NaN mask 到 safe value（exp 用 0，log 用 1）。polynomial path 同样接收 masked input，避免 形成 inf · finite 等无效中间值。
2. log mantissa range 错误。 最初将 IEEE exponent field 设为 127，得到 m ∈ [1, 2)，随后却与 1/√2 ≈ 0.707 比较，因此条件永远不成立。修复为 exponent field 126（m ∈ [0.5, 1)），并相应调整 integer exponent，与 AVX 路径一致。修复前 test_log_float32 报告 3985 ULP error；修复后 774/774 个 exp/log test 全部通过。

测试状态

spin test -- numpy/_core/tests/test_umath.py：4699 passed，0 failed。

后续工作

• 按同一模式实现 f64 exp/log（2-lane LSX，4-lane LASX）。
• 将 npyv_cvt_f32_s32 提升到 NPYV API。
• 增加 LASX（256-bit）backend；由于代码已经采用 NPYV 风格，它会自动拓宽此路径。

2026-05-18 — 启用 Highway LSX target

文件： numpy/_core/meson.build（两处修改）

改动

移除 loongarch64 上已经过时的全局 -DHWY_COMPILE_ONLY_SCALAR flag，并把 LSX 加入 highway_qsort.dispatch.h 与 highway_qsort_16bit.dispatch.h 的 dispatch list。Highway 现在会编译 LSX backend（compiler 支持时也包括 LASX），NumPy 的 sort dispatcher 也会真正实例化 LSX variant。

为什么现在可行

Google Highway 在 v1.2 到 v1.3 之间加入 LSX/LASX target（首个 commit 7e01a07e，2024-11-25；实际实现 7bfb8e8b，2025-04-17；runtime dispatch 10d7ab41，2025-09-09）。bundled submodule 位于 1.3.0-159-gee36c83，同时 包含 hwy/ops/loongarch_lsx-inl.h（5,954 行）和 hwy/ops/loongarch_lasx-inl.h（4,681 行）。

NumPy 的 flag 添加于 2024-11-05（commit 7c35c37a1d），比 Highway upstream 首次加入 LSX stub 还早三周。它会强制所有由 Highway 支持的 op 回退到 scalar， 如今已经过时一年多。

影响

影响范围远大于 exp/log，因为 Highway 在 NumPy 中不仅用于 sort：

• src/npysort/highway_qsort.dispatch.cpp — vectorized quicksort
• src/npysort/highway_qsort_16bit.dispatch.cpp — 16-bit vqsort
• src/umath/loops_trigonometric.dispatch.cpp — sin、cos、tan 等
• src/umath/loops_hyperbolic.dispatch.cpp.src — sinh、cosh、tanh
• src/umath/loops_logical.dispatch.cpp — logical_and、logical_or、logical_not

3A6000 上的实测收益（表格见 BENCHMARKS.zh-CN.md）：

• np.sort f32/int32：约 4–5×（需要第二处修改，将 LSX 加入 sort dispatch list；此前仅枚举 ASIMD、VSX2）
• sin / cos f32：约 3.2×
• tanh f32：约 12.8×（最大单项提升）
• logical_and / logical_not bool：约 8–13×

少数 op 未受影响，因为它们目前并不实际经过 Highway（tan、arcsin、 arctan、sinh、cosh）。

决策

• 暂不把 LASX 加入 sort dispatch list。 meson_cpu/loongarch64/meson.build 中还没有 LASX feature。加入它是独立工作，见后续 LASX backend 章节。
• 保留 # FIXME: disable VXE due to runtime segfault comment。 这是 s390x 问题，与 loongarch64 无关。

测试状态

spin test -- numpy/_core/tests/test_multiarray.py -k sort numpy/_core/tests/test_umath.py：7920 passed，0 failed。

后续工作

• 在 meson_cpu/loongarch64/meson.build 中加入 LASX feature，再将其加入两个 highway_qsort.dispatch.h 及其他 Highway dispatch group，让 3A6000 在 runtime 使用 256-bit variant。
• 调查 tan、arcsin、arctan、sinh、cosh 未提升的原因；它们可能使用 不同的 scalar path（x86 上用 SVML，这里使用 scalar libm）。

2026-05-19 — 定义 LASX feature；dispatch sort/trig/log/hyper

文件：

• meson_cpu/loongarch64/meson.build（加入 LASX feature，并 implies LSX）
• numpy/_core/src/_simd/checks/cpu_lasx.c（新增 probe）
• numpy/_core/src/common/npy_cpu_features.h（NPY_CPU_FEATURE_LASX = 501）
• numpy/_core/src/common/npy_cpu_features.c（HWCAP probe + name registration）
• numpy/_core/meson.build（把 LASX 加入 highway_qsort、loops_hyperbolic、 loops_logical、loops_trigonometric dispatch list）

改动

将 LASX（LoongArch 256-bit SIMD）端到端接入 NumPy CPU dispatch 系统，使 NumPy 编译的所有 Highway-backed code 都能在 3A6000 上于 runtime 使用 LASX target。

detection 仍使用 getauxval(AT_HWCAP) & HWCAP_LOONGARCH_LASX，遵循 HWCAP 是 事实来源而不是 raw cpucfg 的既有规则。3A6000 上的 np.show_runtime() 现在报告： 'simd_extensions': {'baseline': ['LSX'], 'found': ['LASX'], 'not_found': []}。

这只是 LASX dispatch wiring。此时还没有 numpy/_core/src/common/simd/lasx/ NPYV backend，那是另一项更大的工作。当前 LASX 收益完全来自 Highway 自己的 LASX target；NPYV 风格代码（包括五月编写的 f32 exp/log）仍以 LSX 运行。

决策

• LASX implies LSX，写在 meson_cpu/loongarch64/meson.build 中， 与 ARM 的 NEON_FP16 implies NEON pattern 一致。3A6000 两者都具备；实践中 LSX 保持 baseline，LASX 在其上方 dispatch。
• 出现 regression 后，从 highway_qsort_16bit.dispatch.h 移除 LASX。 加入 LASX 后，np.sort(int16) 从 94 降至 59 Melem/s，可能是 Highway upstream gap（其 16-bit vqsort 可能在部分 primitive 上使用不完整 LASX 路径并回退到 scalar）。int16 使用 LSX-only 更好，meson.build comment 记录了原因。
• 未把 LASX 加入 loops_exponent_log.dispatch.h。 exp/log 使用 NPYV 风格， 此时还没有 LASX NPYV backend；加入 LASX 要么没有作用，要么无收益地重复编译。

影响（3A6000，n = 1M）

相对上一节的 LSX-only 状态，LASX 进一步带来：

• np.sort f32：49 → 71 Melem/s（1.44×）
• np.sort f64：23 → 33 Melem/s（1.45×）
• np.sort int32：62 → 81 Melem/s（1.30×）
• np.sort int64：37 → 42 Melem/s（1.13×）
• np.sin f32：250 → 435 Melem/s（1.74×）
• np.cos f32：246 → 437 Melem/s（1.78×）

tanh 与 logical op 在 LSX 下已经饱和（即便 NumPy 层 dispatch target 是 LSX， Highway 也会在内部采用这些路径）。tan、arcsin、arctan、sinh、cosh 仍无变化，它们使用其他 code path。

相对最初 scalar baseline，累计收益为：

• np.sort int32：12.7 → 81 Melem/s（6.4×）
• np.sort f32：11.5 → 71 Melem/s（6.2×）
• np.sin f32：77 → 435 Melem/s（5.7×）
• np.cos f32：78 → 437 Melem/s（5.6×）

测试状态

spin test -- numpy/_core/tests/test_multiarray.py -k sort numpy/_core/tests/test_umath.py：7920 passed，0 failed。

后续工作

• 原生 LASX NPYV backend（numpy/_core/src/common/simd/lasx/）是下一项大 工作。它会把已经写好的 f32 exp/log（以及其他 NPYV loop）从 4 lane 拓宽到 8 lane。
• 调查 Highway upstream 中 int16 LASX qsort regression。
• tan、arcsin、arctan、sinh、cosh 仍为 scalar；检查 loops_trigonometric.dispatch.cpp 与 loops_hyperbolic.dispatch.cpp.src 的实际 code path。

2026-05-20 — 原生 LASX NPYV backend

新增文件：

• numpy/_core/src/common/simd/lasx/lasx.h（umbrella、type、lane count）
• numpy/_core/src/common/simd/lasx/misc.h
• numpy/_core/src/common/simd/lasx/memory.h
• numpy/_core/src/common/simd/lasx/reorder.h
• numpy/_core/src/common/simd/lasx/operators.h
• numpy/_core/src/common/simd/lasx/conversion.h
• numpy/_core/src/common/simd/lasx/arithmetic.h
• numpy/_core/src/common/simd/lasx/math.h

总计 2108 行，与 LSX backend 的文件结构逐一对应。

修改文件：

• numpy/_core/src/common/simd/simd.h — LASX dispatch entry 位于 LSX 之前
• numpy/_core/src/umath/loops_exponent_log.dispatch.c.src — 扩展 npyv__cvt_f32_s32 helper，在 LASX target 使用 __lasx_xvffint_s_w
• numpy/_core/meson.build — 将 LASX 加入 loops_exponent_log dispatch

改动

为 LoongArch LASX 编写完整的 NPYV（256-bit）backend。从 simd/lsx/ 转换的工作 大部分是 mechanical：

• __m128i → __m256i，__m128 → __m256，__m128d → __m256d
• __lsx_v* → __lasx_xv*
• lane-count constant 加倍（nlanes_f32 4→8 等）
• vector literal type：v4i32 → v8i32，v4f32 → v8f32， v2f64 → v4f64 等
• partial load/store switch 从 1–3 case 扩展到 1–7 / 1–3 case，并明确 fallthrough

2026-05-17 编写的 f32 exp/log kernel 已采用 NPYV 风格，因此将 LASX 加入其 dispatch list 后，它自动获得 256-bit lane，无需修改 source。

决策与注意事项

• 不存在 xvpickve2gr_h / xvpickve2gr_b。 LASX 只提供 _w 与 _d element-extract intrinsic。operators.h、conversion.h 和 math.h 必须通过 xvpickve2gr_wu 加 shift/mask，从所属 32-bit lane 中取得 16-bit 或 8-bit lane；首次 build 即发现此问题。
• xvmsknz_b 按 128-bit half 输出 byte-nonzero mask，并放在各 half 的低 16 bit（即 low half 的 word lane 0 与 high half 的 word lane 4）。any/all reduction 与 tobits_b8 会把两者组合成 32-bit scalar mask。
• 先在 half 内 reduction，再组合。 LASX shuffle/horizontal-add 与 AVX 一样， 都在各 128-bit half 内执行，不会自动跨完整 256 bit。所有 npyv_reduce_* 与 npyv_sum_* 先用 LSX pattern 分别 reduction 两个 half，再以一次 scalar op 合并。
• unzip 通过 memory round-trip 回退到 scalar。 cross-half permute 需要两个 带 custom immediate 的 xvpermi.q。这可以实现，但 unzip 不在 exp/log hot path 上，因此在真正有 user code 高频使用前，scalar fallback 足够。

测试状态

spin test -- numpy/_core/tests/test_umath.py：4699 passed，0 failed。

影响（3A6000，n = 1M，float32）

两项在 LASX 上都相对 LSX 获得接近 2× 的提升；n 足够大、scalar tail 的开销被 摊薄后，8-lane scaling 如预期实现。

后续工作

• f64 NPYV exp/log kernel；现在已经有 4-wide LASX f64 lane，只缺 kernel code。
• 将 npyv_cvt_f32_s32 与 byte/half-lane extract helper 提升到 shared NPYV API， 让其他架构也能干净地复用。
• 审核 LASX unzip/cross-lane path；若 user code 开始高频使用，则改为基于 xvpermi.q 的实现。

2026-05-20 — f64 NPYV exp / log

文件： numpy/_core/src/umath/loops_exponent_log.dispatch.c.src

改动

加入 double-precision exp/log 的 NPYV 风格 kernel，gate 为 #if NPY_SIMD_F64 && !defined(SIMD_AVX2_FMA3) && !defined(SIMD_AVX512F) && !defined(SIMD_AVX512F_NOCLANG_BUG)，并加入 NPYV_IMPL_F64_EXP_LOG sentinel。DOUBLE_@func@ dispatcher 在 AVX-512 path 与 scalar fallback 之间新增 #if defined(NPYV_IMPL_F64_EXP_LOG) 分支。

算法

f64 需要约 15 位十进制精度，因此不能沿用过短的 f32 P5/Q2 polynomial。这里采用 两种成熟算法：

• exp： Cephes rational approximation。先做 Cody-Waite range reduction： xr = x - round(x*log2e)*ln2（将 ln2 拆成 C1 + C2，使每次乘法都 exact）， 再计算 exp(xr) = (1 + 2z)，其中 z = xr·P(xr²) / (Q(xr²) − xr·P(xr²))。P 在 xr² 中为 degree-2，Q 为 degree-3。最终通过 IEEE-754 exponent-field bit trick（bias 1023，shift 52） 乘以 2^k。coefficient 是 inline literal 形式的 Cephes constant，未放入 npy_simd_data.h。
• log： FDLIBM 风格 polynomial。分解 x = 2^k · m，其中 m ∈ [√0.5, √2]（将 IEEE exponent field 设为 1022；若 m ≤ 1/√2 再 shift）。 随后令 f = m − 1、s = f/(2 + f)，计算 s² 与 s⁴ 中的 degree-7 polynomial，并与 f 和 k·ln2 组合（ln2 拆成 hi/lo 以保证精度）。七个 Lg* coefficient 是标准 FDLIBM value。

同一文件中新增两个 arch-specific helper：

#define npyv__cvt_f64_s64(X)  ((npyv_f64)__lasx_xvffint_d_l(X))
#define npyv__cvt_s64_f64(X)  ((npyv_s64)__lasx_xvftintrne_l_d(X))

并提供对应 LSX entry。NPYV 尚未 portable 地暴露 int64↔float64 cast；TODO 记录了后续提升计划。

测试状态

spin test -- numpy/_core/tests/test_umath.py：4699 passed，0 failed。

在 200,000 个 random sample 上相对 math.exp/math.log 的 ULP：

结果优于 f32 path；后者在 edge case 上的 max 约为 6 ULP。

影响（3A6000，n = 1M，float64）

同一 kernel 也会生成 2-wide LSX 与 4-wide LASX f64；显著数字采用 3A6000 runtime dispatch 的 LASX。bench script 没有单独记录 LSX-only f64。

后续工作

• 将 f64 Cephes coefficient 提升到 npy_simd_data.h，避免 source 中的 literal duplication。
• 在 f32 kernel correctness 也经过大规模 LASX test 后，把 NPYV_IMPL_F32_EXP_LOG / NPYV_IMPL_F64_EXP_LOG 重构为一个 guard。

2026-05-23 — Tier 1 NPYV transcendental（sinh/cosh/exp2/log2/log10/arccosh）

文件：

• numpy/_core/src/umath/npyv_exp_log.h（新增）— exp/log kernel 从 loops_exponent_log.dispatch.c.src 移至此处以便共享
• numpy/_core/src/umath/loops_exponent_log.dispatch.c.src — 改为 #include header，不再 inline 定义 kernel
• numpy/_core/src/umath/loops_umath_fp.dispatch.c.src — 在 exp/log 上新增 12 个 NPYV wrapper kernel（6 op × 2 precision），并在 unary dispatcher template 中 加入 NPYV branch
• numpy/_core/src/umath/loops.h.src — 声明 power/arctan2 前重新 include loops_umath_fp.dispatch.h，避免它们继承不含 LASX/LSX 的 loops_half dispatch list，导致 LASX build 失败
• numpy/_core/meson.build — loops_umath_fp.dispatch.h 中 [X86_V4] → [X86_V4, LASX, LSX]

改动

把 LoongArch acceleration 扩展到六个原先位于 loops_umath_fp.dispatch.c.src 的 transcendental。此前它们只有 X86_V4/SVML 路径，在包括 LoongArch、ARM、PowerPC 在内的其他架构上都使用 scalar libm。 六项都以现有 exp/log NPYV kernel 上的 thin wrapper 实现：

kernel naming 遵循现有 dispatcher template 使用的 libm convention （asinh/acosh/atanh，而非 NumPy 的 arcsinh 形式），因此 dispatcher 中同一 #elif defined(NPYV_HAVE_@intrin@_@sfx@) pattern 可以自动接入，无需 per-op #ifdef。

架构重构

exp/log kernel 原先 inline 在 loops_exponent_log.dispatch.c.src 中。它们被移至 新的 numpy/_core/src/umath/npyv_exp_log.h，让 loops_umath_fp 可以在其上组合。 两个文件都 include 该 header；kernel 为 NPY_FINLINE，所以每个 translation unit 获得自己的 inline copy，不产生 link conflict。NPYV_IMPL_F32_EXP_LOG / NPYV_IMPL_F64_EXP_LOG marker 现在在 header 内设置（前提是对应 int↔float cast macro 可用）；原文件 guard 仅检查 marker，不再重复 condition。

决策

• Tier 1 有意不包括 arcsinh 与 arctanh。 自然 identity （log(x + √(x² + 1)) / ½·log((1+x)/(1−x))）在 |x| → 0 时会发生 catastrophic cancellation。例如 f64 中 (1 + 1e-20) round 为 1.0，log 返回 0 而不是约 1e-20。direct formula 会让 TestComplexFunctions::test_loss_of_precision 在 [1e-20, 1e-3] 全范围失败。 仅用 Taylor 修复，在 |x|=0.5 达到精度需要 30+ term。正确修复是 log1p， 当时 NPYV 尚无此实现。marker macro 被 comment，并记录原因；kernel function 作为 dead code 保留，等待 log1p 完成后重新启用。
• Tier 1 包括 arccosh。 x = 1 边界上 (1+ε)² − 1 的精度损失没有现有 test 覆盖，且只影响很窄的 corner。direct formula 可通过 test suite。
• 在声明 power/arctan2 前重新 include loops_umath_fp.dispatch.h。 dispatch macro 是 global state，每次 include dispatch.h 都会重写。改动前， power/arctan2 紧接 loops_half.dispatch.h 声明，因此继承 loops_half 的 dispatch list（无 LASX）。build 期待 FLOAT_power_LASX 存在，却从未声明它。 重新 include 正确 header 后恢复 macro。这实际是由本次改动暴露出的 upstream NumPy bug，未来 upstream 时值得说明。

测试状态

umath suite：4699 passed，0 failed。

影响（3A6000，n = 1M）

在同一 build 中对比 contiguous array（NPYV path）与同长度 strided array。strided dispatch 无法通过 steps[i] == sizeof(dtype) 检查，因而回退到 scalar libm； 编译代码相同，仅 code path 不同。

arccosh 提升较小，因为 scalar libm acoshf 已经可以 vectorize，或 LoongArch libm 本身已有优化。其他 op 此前确实受 scalar path 限制。

后续工作

• Tier 2 transcendental： tan、arcsin、arccos、arctan、 arctan2、expm1、log1p、cbrt、power(x, y)。每项都需要自己的 polynomial 与严谨的精度处理。log1p 是重新启用 arcsinh/arctanh 的前提。
• power(x, y) 可在 x > 0 时用 exp(y·log(x)) 实现，但 integer-exponent fast path 以及 0^0 / negative base edge case 都需谨慎处理，并不简单。

2026-05-23 — Tier 2 NPYV transcendental（除 power 外的所有项）

文件：

• numpy/_core/src/umath/npyv_exp_log.h — 为 f32 与 f64 加入 log1p 和 expm1 kernel。两者在较小 |x| 时使用 Taylor，避免 1+x 或 e^x - 1 的 catastrophic cancellation；中等 |x| 使用既有 log_kernel(1+x) / exp_kernel(x) - 1。f64 使用 7-term Taylor，f32 使用 5-term。
• numpy/_core/src/umath/loops_umath_fp.dispatch.c.src — 新增 8 个 NPYV kernel：arctan（FDLIBM polynomial + 3-way range reduction）、arcsin、 arccos（基于 arctan 组合）、arctan2（binary + quadrant logic）、tan （按 π/2 做 Cody-Waite reduction + FDLIBM polynomial）、cbrt （sign(x) · exp(log(|x|)/3)）。同时改写 arcsinh 与 arctanh，让它们使用 log1p，修复 small-|x| 精度问题并重新启用 NPYV。

改动

此 pass 补齐 loops_umath_fp.dispatch.c.src 中的大部分缺口。该文件拥有的 17 个 ufunc 中，15 个现在已有 NPYV 实现；只剩 power(x, y)。它涉及 negative base 配 integer exponent、0^0、任一 operand 为 ±inf 等 corner case，明显更复杂。

kernel 可分为几种 pattern：

• range reduction 后直接 polynomial： arctan（按 tan(π/8) 与 tan(3π/8) 做 3-way reduction，再用 f32 degree-5 / f64 degree-21 FDLIBM minimax）；tan（按 π/2 做 Cody-Waite reduction，再用 f32 degree-11 / f64 degree-25 FDLIBM minimax）。
• 在 exp/log/log1p 上组合： arcsinh = log1p(|x| + x²/(1+sqrt(1+x²)))（使用 rationalized form，否则 1+x²→1 cancellation 会破坏 small-|x|）； arctanh = ½·log1p(2x/(1-x))；cbrt = exp(log(|x|)/3)；以及 exp2、 log2、log10、sinh、cosh 的前述 identity。
• 在 arctan 上组合： arcsin = arctan(x/sqrt(1-x²))； arccos = π/2 - arcsin；arctan2(y, x) = arctan(y/x) 再加 quadrant adjustment。

决策与注意事项

• 巨大 |x| 的 tan。 Cody-Waite reduction 在约 |x| < 2²⁰ 内精确；更大时 r = x - k·π/2 error 会增长。实际 workload 很少到达此范围，test suite 通过。 Payne-Hanek 可扩展精度，但需要大量代码。
• 最常见 failure 是不应出现的 FP exception。 SIMD 会计算所有 branch，即使 最后通过 select 丢弃，因此除零、inf/inf、sqrt(<0) 等都会设置 flag。修复 pattern 是：对于不会选择该 branch 的 lane，先把 input mask 到 safe value （1.0 或 0.0）。这修复了 arctan、arcsin、arctan2 与 tan 的多处问题。
• 由于 log1p，arcsinh 与 arctanh 重新启用 NPYV。 前述 identity 在 small |x| 不损失精度；Tier 1 的 disable 被撤回，marker 恢复为 NPYV_HAVE_asinh_* / NPYV_HAVE_atanh_*。
• 巨大 |x| 的 arcsinh 会在 moderate path 中使 x² overflow。对 |x| > 2^15（f32）/ 2^27（f64）改用 asymptotic arcsinh(x) ≈ log(|x|) + ln(2)。
• arctan2 memory-overlap check。 dispatcher 进入 NPYV path 前必须调用 is_mem_overlap，因为 arctan2.accumulate 会让 dst 与一个 input overlap。 该问题由 test_memoverlap_accumulate_symmetric 发现。
• 推迟 power。 完整实现还需处理 negative x + integer y（sign 取决于 parity）、0^0 = 1、±inf^(...)、(...)^±inf、保留 sign 的 ±0 等。

测试状态

spin test -- numpy/_core/tests/test_umath.py：4699 passed，0 failed。

影响（3A6000，contiguous NPYV vs strided scalar libm，n=1M）

数字在不同 run 间有噪声，特别是 absolute speed 较低、overhead 占主导时；下列 ratio 已 round。

tan 最突出：3A6000 上 scalar libm tan 约为 7 M elem/s，而 NPYV 一次处理 8 lane。arcsin / arccos 提升最小，因为它们调用本身已经承担 polynomial cost 的 arctan；收益主要来自 SIMD width。

后续工作

• power(x, y) — 此文件中最后尚未实现的 op。
• Tier-2 polynomial pivot（arctan 的 tan_pi8 / tan_3pi8 constant、tan 的 π/2 split）目前 inline，放入 npy_simd_data.h 会更整洁。
• arcsin f64、arccos f64、cbrt f32 等提升较小；LoongArch libm 的这些实现 已经很高效。

2026-05-23 — power(x, y) NPYV

文件： numpy/_core/src/umath/loops_umath_fp.dispatch.c.src

改动

实现 npyv_pow_FLOAT_kernel 与 npyv_pow_DOUBLE_kernel，补齐 loops_umath_fp.dispatch.c.src 最后一处缺口。该文件每个 op 现在都在 LoongArch （LSX 与 LASX）上拥有 NPYV path。

算法

通用公式为 |x|^y = exp(y · log(|x|))。外围用一系列 mask-based override 处理 IEEE 754 / C99 pow() special case：

避免无关 FP flag

• log(0) 返回 -inf；若 y=0，随后 y·-inf 会产生带 INVALID 的 NaN。因此在 log 前将 x=0 lane mask 为 |x|=1。
• 0 · ±inf 会产生 INVALID；在 multiplication 前把 |x|=1 lane 的 y mask 为 0。
• round_to_int(±inf) 会 saturation 并触发 INVALID；在 int conversion / parity check 前将 non-finite y lane mask 为 0。

决策

• parity check：f32 使用 round_s32_f32，f64 使用 even/odd fractional check。 f32 path 的 int32 saturation 意味着 y > 2^31 时 pow(-2, y) 会 丢失 parity，极大 exponent 会被默认为 even。f64 使用更稳健的 y/2 - rint(y/2) == 0 ⇒ even，适用于所有 finite y。
• pow.accumulate memory overlap 在 dispatcher 中检查，与 arctan2 相同。
• pow(0, -1) = +inf 触发 DIV_BY_ZERO 是符合 IEEE 的预期行为，与 NumPy/libm 一致。

测试状态

4699 个 umath test 全部通过，无 regression。

影响（3A6000，n = 1M）

提升不算巨大，因为 LoongArch libm pow 内部已经采用 exp(y·log(x)) 组合；但 f64 收益仍有意义，也意味着 loops_umath_fp 的全部 17 个 op 现已在 LoongArch 上经过 NPYV routing。

LoongArch 上 loops_umath_fp.dispatch.c.src 的最终状态

17 个 unary + 2 个 binary = 19 个 ufunc，全部使用 NPYV。

2026-05-23 — Regression cleanup：arcsin / arccos f64 direct kernel

文件： numpy/_core/src/umath/loops_umath_fp.dispatch.c.src

改动

Tier-2 pass 让 f64 arcsin / arccos 在 3A6000 上反而慢于 scalar libm：

组合方式为什么失败

原 kernel 通过已有 Tier-2 arctan 计算 arcsin(x) = arctan(x/√(1-x²))。每个 element 需要：

• argument transform 的 1 次 sqrt 与 1 次 div；
• 完整 arctan f64 kernel：3-way reduction + degree-21 FDLIBM polynomial。

4-lane LASX f64 的 SIMD parallelism 无法弥补约为 scalar libm direct asin 两倍的 floating-point 工作量。

修复

使用 FDLIBM two-region rational 重写 npyv_asin_DOUBLE_kernel，coefficient 来自 glibc sysdeps/ieee754/dbl-64/e_asin.c：

• |x| < 0.5：asin(x) = x + x·t·P(t)/Q(t)，t = x²
• |x| ≥ 0.5：asin(x) = sign(x)·(π/2 - 2·s·(1+w))， s = √((1-|x|)/2)，w = t·P(t)/Q(t)，使用同一 polynomial

两个 branch 共用 P（degree 5）/ Q（degree 4），所以 polynomial 无条件计算。 large-path identity 会把 argument 移回同一 polynomial 精确的范围。总成本为一次 sqrt、一次 div、一个 degree-5 Horner 与一个 degree-4 Horner，不含 transcendental。

npyv_acos_DOUBLE_kernel 仍组合为 π/2 - asin(x)，因而自动获得提升。

结果

• arcsin f64：34.3 → 269.4 M elem/s（0.51× → 相对 main 3.97×）
• arccos f64：52.0 → 264.6 M elem/s（0.80× → 相对 main 4.08×）
• 46 个 op 的总体 geomean：3.09× → 4.25×
• 相对 main 的 regression：5 → 0
• 4699 个 umath test 全部通过

决策：不同时重写 f32 arcsin/arccos

f32 kernel 仍在 arctan 上组合，相对 main benchmark 为 2.07× / 2.27×，已有收益且 没有精度问题。f64 fix 不代表必须修改 f32，因此保持不变。

2026-05-23 — 通过精度 gate 重新启用 SIMD f64 sin/cos

文件： numpy/_core/src/umath/npyv_sincos.h（新增）、 numpy/_core/src/umath/loops_trigonometric.dispatch.cpp

背景：上游为什么禁用此路径

upstream loops_trigonometric.dispatch.cpp 215–217 行完全禁用了 SIMD f64 sin/cos：

/* Disable SIMD code sin/cos f64 and revert to libm: see
 * https://mail.python.org/archives/list/[email protected]/
 *   thread/C6EYZZSR4EWGVKHAZXLE7IBILRMNVK7L/
 * for detailed discussion on this*/

讨论涉及一个曾发布的 SVML kernel：当 |x| ≥ 2^20 时会损失精度。三段 π/2 split 的 Cody-Waite reduction 在约 2^20 之后失效，而 kernel 没有 Payne-Hanek fallback。 因此 DOUBLE_sin / DOUBLE_cos 在所有平台都 route 到 scalar libm；在 3A6000 上这会放弃约 3× LASX speedup。

改动

只在 Cody-Waite 可证明精确的范围内重新启用 SIMD f64 sin/cos。per-block gate 在 block 含任何 out-of-range lane 时回退到 scalar libm；pattern 与同文件 Highway f32 sin/cos 完全一致。

• npyv_sincos.h — 单个 FDLIBM polynomial kernel npyv__sincos_f64_kernel(x, want_cos)。它使用三段 π/2 split 做 Cody-Waite reduction，计算 FDLIBM __kernel_sin / __kernel_cos polynomial，再执行 quadrant select / sign flip。
• loops_trigonometric.dispatch.cpp — DISPATCH_DOUBLE_FUNC(func, OPCODE) 现在会包裹 SIMD call（仅 LoongArch + NPYV f64）。对每个 npyv_nlanes_f64 block，dispatcher 通过 npyv_all_b64(cmple) 检查每个 lane 是否 |x| ≤ 2^20。全部在范围内时运行 SIMD，否则逐 lane 调用 scalar npy_sin / npy_cos。非 LoongArch 上 macro 仍是原 upstream form。

精度审计（3A6000，每个范围 100k sample）

SIMD range 的 max 为 1 ULP，达到或优于 FDLIBM 约 2 ULP 的目标。gate 之外的 libm fallback bit-exact，因为 strided 与 contiguous path 都调用同一 npy_sin/npy_cos。kernel 内部先把 NaN mask 到 0 再 rint，最后恢复 NaN。

结果

• sin f64：68.5 → 240.5 M elem/s（1.00× → 相对 main 3.51×）
• cos f64：67.9 → 235.9 M elem/s（1.00× → 相对 main 3.47×）
• 46 个 op 的总体 geomean：4.25× → 4.56×
• 不再有 neutral 或 regression：每个 op 相对 main ≥ 1.10×；最低为 f64 cbrt 的 1.42×
• 4699 个 umath test 全部通过

决策：hybrid gate，而非完整 SIMD reduction

选择 per-block |x| < 2^20 libm fallback，而未立即实现 SIMD Payne-Hanek：

• 常见数据（magnitude 约 10⁶ 以内）获得完整 SIMD speedup。实际 workload 很少 超出此范围，angle array 通常已经 wrap 到 [-π, π] 附近。
• SIMD Payne-Hanek 会增加约 150 行代码，并引入新的精度调优风险；对 tech-demo branch 而言当时得不偿失。
• fallback 是正确的 libm-exact 结果，不是 approximation；巨大 argument 仍获得 原 NumPy behavior。

后续：SIMD Payne-Hanek

TODO：实现 SIMD Payne-Hanek range reduction，让 SIMD f64 sin/cos 覆盖完整 f64 range，不再回退到 libm。这样可以：

• 消除 per-block max(|lane|) check（目前约一次 compare + branch / block）；
• 让确有巨大 argument 的 workload（如 unwrapped phase data）也获得 SIMD speedup；
• 使 LoongArch path 更接近可提交 upstream PR 的状态，而不是 LoongArch-specific opt-in。

参考：Payne 与 Hanek 的 1983 年论文 “Radian Reduction for Trigonometric Functions”；FDLIBM k_rem_pio2.c scalar reference；以及 simd_sincos_f32 顶部标记同一 f32 follow-up 的 Highway TODO。

决策：非 LoongArch 保持 upstream behavior

新 macro body 由 #if NPY_SIMD && defined(__loongarch__) && defined(NPYV_IMPL_F64_EXP_LOG) guard， 因此非 LoongArch build 得到原 upstream UNARY_LOOP { *op1 = npy_##func(*ip1); }。这是有意为之：在没有 benchmark 的 x86 上重新启用 SIMD f64 sin/cos 会绕过 numpy-discussion 的决定。LoongArch 是 opt-in；若要扩展到其他架构，只需调整一行 guard。

2026-05-23 — cbrt f64：FDLIBM bit-magic + Newton

文件： numpy/_core/src/umath/loops_umath_fp.dispatch.c.src

背景

sin/cos 工作完成后，f64 cbrt 是 speedup 最低的 op，相对 main 1.42×（受噪声影响 有时为 1.29×）。此前 kernel 在 NPYV exp/log 上组合：

cbrt(x) = sign(x) · exp(log(|x|) / 3)

每个 element 要为 log_kernel 付出约 30 op，再为 exp_kernel 付出约 30 op，加上 乘以 1/3 与重建 sign，合计约 60 op。libm 会绕过两个 transcendental。

算法（FDLIBM s_cbrt.c）

1. 通过 bit hack 得到 initial estimate（约 5-bit accuracy）。取得 |x| 的 top 32 bit，除以 3，加 magic constant B1 = 715094163 ≈ (682 - 0.033) · 2²⁰。对 biased exponent field 的 arithmetic 隐式把 binary exponent 除以 3，B1 把结果调整到正确 cube-root exponent。再 shift 回 high 32 bit，low bit 清零，reinterpret 为 double：

c t = bits_to_double(((|x|_bits >> 32) / 3 + B1) << 32)

2. Rational refinement（约 5 → 23 bit）。使用 FDLIBM 对 r = t³/x 的 standard rational：

c r = t·t/x; s = C + r·t; t *= G + F/(s + E + D/s);

C、D、E、F、G 分别为 19/35、-864/1225、99/70、45/28、5/14。

3. Newton iteration 达到 53 bit。 先清零 t 的 low 32 bit，使 t² exact：

c s = t·t; r = x/s; w = t+t; r = (r-t)/(w+r); t = t + t·r;

SIMD 细节

• top 32 bit 的 bit-hack 除以 3 在 FP 中完成：把小于 2³¹、可在 f64 中 exact 表示的 hx 转为 double，乘以 1/3，调用 npyv_floor_f64，再转回 integer。 这样无需 SIMD integer-divide-by-3 或 widening multiply。
• subnormal、±0、±inf、NaN 通过 special mask 绕过 bit-magic，lane 暂时直接 传递原 x。这对 ±0、±inf、NaN 正确，但 subnormal 当时仍不正确；其所需的 FDLIBM rescale branch（乘 2⁵⁴、运行算法、除 2¹⁸）列为 follow-up。dispatcher 中含 subnormal 的 lane 需要回退到 libm，strided path 已如此。
• rational refinement 与 Newton 之间清掉 t 的 low 32 bit，是 FDLIBM 保证 t·t 在 FP 中 exact 的标准技巧。SIMD 以 0xFFFFFFFF00000000 bitwise AND 实现。

精度审计（3A6000，每个范围 50k sample，相对 scalar libm）

完整 f64 normal range max 为 3 ULP，略宽于 FDLIBM 2 ULP 目标。额外 ULP 来自 FP-based initial divide-by-3；FDLIBM 对 32-bit hx 使用 integer divide，而这里的 floor(hx · (1/3)) 引入 sub-ULP 偏移，随后被 rational refinement 放大。可以通过 给 NPYV 加 integer-divide-by-3 helper 收紧，但 3 ULP 仍在 SIMD transcendental 预算内。

special value 已验证：cbrt(±0) = ±0、cbrt(±inf) = ±inf、cbrt(NaN) = NaN。 test_umath.py 中 86 个 cbrt test 全部通过。

结果

• cbrt f64：42.7 → 140.6 M elem/s（1.42× → 相对 main 3.29×）
• 46 个 op 的总体 geomean：4.56× → 4.71×
• bottom-10 floor 从 1.42× 提高到 1.55×；每个 benchmark op 现在都比 scalar libm 快至少 1.55×

决策：不同时重写 f32 cbrt

f32 cbrt 相对 main 已有 3.57×；现有 exp(log(|x|)/3) 在 8-lane LASX f32 path 上表现良好。重写的收益更少且需要独立 ULP audit，因此保持现状。

后续：~~subnormal handling~~ — 已完成

最初 cbrt 实现将 subnormal route 到原 x，结果错误；后续章节补上 FDLIBM rescale branch。

2026-05-23 — cbrt f64：subnormal correctness

文件： numpy/_core/src/umath/loops_umath_fp.dispatch.c.src

问题

cbrt 的 bit-magic initial estimate 假设 input 为 normal magnitude；biased exponent field 必须落在 normal double range。若 exponent field 为零（subnormal）， divide-by-3-plus-magic arithmetic 会失效。首次实现把 subnormal lane route 到 x 自身，这是错误的：cbrt(2⁻¹⁰⁷⁴) ≈ 1.7e-108，不是 2⁻¹⁰⁷⁴。

修复

按照 FDLIBM s_cbrt.c：

1. 先把 subnormal lane 乘以 2⁵⁴，推入 normal range。
2. 这些 lane 使用 magic constant B2 = 696219795，而不是 B1 = 715094163。 B2 = B1 − 18·2²⁰；initial estimate exponent field 中的 −18 会精确抵消 2⁵⁴ scaling，因为 cbrt(2⁵⁴) = 2¹⁸。
3. rational refinement 与 Newton step 使用 unscaled |x|，因为它们针对原始 cube-root constraint refinement t³ ≈ |x|。

SIMD 在 bit-magic 前按 lane select：subnormal 选择 bit_input = abs_x × 2⁵⁴ 与 magic = B2，其他 lane 选择原 abs_x 与 B1。 多出两个 select 与一个 multiply，在 4-lane LASX path 上几乎无成本；throughput 仍在 137 M/s 噪声范围内。

精度审计（完整 subnormal range 的 156 个 sample）

最小 subnormal cbrt(2⁻¹⁰⁷⁴) = 1.703184e-108 与 scalar libm 完全一致。 subnormal 反而达到 1 ULP（normal 为 3 ULP），因为 scaled input 给 bit-magic 更干净的起点，FP divide-by-3 也有更多 headroom。

已验证 edge case

4699 个 umath test 仍全部通过。相对 main throughput 保持约 3.2×。

2026-05-24 — log1p f64：按 block dispatch 的直接 FDLIBM Padé

文件： numpy/_core/src/umath/npyv_exp_log.h

背景

log1p f64 相对 main 为 1.95–2.17×——当 |x| 超过非常严格的 2⁻⁷ Taylor threshold 时，kernel 的 slow path 会经过 log_kernel(1+x)（约 30 个操作）。 由于 7 项 Taylor branch 和 log_kernel 总是一起计算再 blend，每次调用都要承担完整的 log_kernel 成本。

算法——两个 kernel，按 block dispatch

• npyv__log1p_pade_f64——FDLIBM Lp1..Lp7 Padé 形式 （即 Lg1..Lg7，Mercator series 的 2/(2k+1) coefficient）：

s = x/(2+x); z = s²; hfsq = x²/2 R = z · poly(z) [Horner, 7 coefficients] log1p(x) = x - (hfsq - s·(hfsq + R))

成本约 15 个操作，在 x ∈ (1-√2/2, √2-1) ≈ (-0.293, 0.414) 内有效。超出该区间时， 公式会损失约 14 bit。

• npyv__log1p_full_f64——log_kernel(1+x) 加上 FDLIBM 的 correction term。correction 捕获 1+x 的舍入误差，因此即使 x 小到 1+x 舍入为 1，log1p 仍能保持 ULP 精度。根据 u = 1+x ≥ 2 与否， 通过 per-lane select 在两种 correction 形式之间选择。

• 顶层 kernel——通过 npyv_all_b64 按 block 检查：若每个 lane 都在 Padé 区间内，就走 cheap path；否则走 full path。检查只需一次 cmpgt + cmplt + and 和一个 scalar branch；对于任何单调或子区间 workload， branch 都很好预测。

踩坑记录

第一次尝试使用了对称 threshold |x| < 0.4。这会让 x = -0.39 走 Padé， 但此处违反了公式的底层假设（u = 1+x 应落在 [√2/2, √2) 内）。结果是在 测试覆盖的区域出现约 8700 ULP error。把 threshold 修正为 FDLIBM 使用的非对称 (-0.29, 0.41) 区间后，整个 f64 normal range 的 max ULP 回到 1。

full path 还需要 inf/NaN masking——correction 中的 inf − inf 会触发虚假的 INVALID（被 NumPy 的 test_unary_spurious_fpexception 捕获）。修复方法是在 subtraction 之前把 non-finite lane mask 为 0，然后让 log_kernel(1+x) 自己的 special-case path 产生正确的 inf/NaN。还加入了 u == 0 guard，防止 log1p(-1) 在 correction 中计算 0/0（这一问题在 arctanh(-1) 开始返回 NaN 而不是 -∞ 时被发现）。

精度审计（每个 range 50k 个 sample，与 scalar libm 对比）

special-value 检查：log1p(0) = 0、log1p(-1) = -∞、log1p(-2) = NaN、 log1p(∞) = ∞、log1p(NaN) = NaN——全部正确，且没有虚假 flag。

结果

• log1p f64：63.4 → 231.1 M elem/s（相对 main 从 1.95× → 3.65×）
• arctanh f64 也有提升（原为 2.54×），因为它的 slow path 组合使用了 log1p。见下一节。

2026-05-24 — arctanh f64：以 x² 为变量的直接 Cephes Padé，按 block dispatch

文件： numpy/_core/src/umath/loops_umath_fp.dispatch.c.src

背景

arctanh f64 相对 main 为 2.54×。kernel 使用 0.5 · log1p(2x/(1-x))——即便有了上面更快的 log1p，仍要付出一次 divide 和完整 log1p path 的成本。

算法

当 |x| ≤ 0.5 时，函数使用以 t = x² 为变量的 Cephes Padé 近似：

arctanh(x) = x + x · t · P(t) / Q(t)

其中 P 为 4 次、Q 为 5 次（最高项 1 隐含）。coefficient 来自 Cephes atanh.c，peak relative error 为 1.4e-17，成本约 13 个操作。

当 |x| > 0.5 时，polynomial 不收敛，因此使用原来的 log1p 公式。 per-block dispatch 根据 max(|lane|) ≤ 0.5 选择 path。

经验

Cephes 存储 polynomial coefficient 时，coef[0] 对应的是最高次幂， 而不是最低次幂。第一次尝试以错误顺序执行 Horner，产生了离谱的结果 （约 10¹⁸ ULP error）。已在 kernel comment 中记录，以免移植其他 Cephes routine 时重犯。

精度审计（每个 range 50k 个 sample）

slow-path 边界处的 21 ULP 来自 x = 1 附近 (1-x) 丢失 ULP 时的渐近 放大；libm 本身在这里也只略好一些。

special value：arctanh(±0) = ±0、arctanh(±1) = ±∞、 arctanh(|x| > 1) = NaN、arctanh(NaN/±∞) = NaN——全部正确。

结果

• arctanh f64：33.0 → 485.2 M elem/s——相对 main 为 14.7× （rewrite 前为 2.54×）。bench range 完全落在 Padé 区间内，因此几乎每个 block 都走 fast path。
• 46 个 op 的 overall geomean：4.71× → 4.95×

2026-05-24 — f64 sin/cos slow path 的 Payne–Hanek reduction

文件： numpy/_core/src/umath/payne_hanek_f64.h（新增）、 numpy/_core/src/umath/npyv_sincos.h、 numpy/_core/src/umath/loops_trigonometric.dispatch.cpp

背景

2026-05-23 合入的 sin/cos f64 实现只在 max(|lane|) < 2^20 时使用 SIMD， 只要 block 中含有更大的 input 就 fallback 到 scalar libm。2^20 限制来自 Cody–Waite range reduction：y = round(x · 2/π) 必须能在 53-bit mantissa 中精确表示；对更大的 |x|，x · inv_pi_2 的舍入误差会破坏 y。超过这一范围后， libm 使用 Payne–Hanek——一种 table-based reduction：存储 2/π 的许多 bit， 并按 x 的 exponent 建立索引，从而处理任意 magnitude。

这次改动用我们自己的 Payne–Hanek reduction 替代 libm fallback。LoongArch sin/cos f64 path 对 libm 的依赖降为零，并以 libm 等效精度覆盖完整 f64 normal range。

实现

新增两个文件；slow-path kernel 位于二者之间：

• payne_hanek_f64.h——scalar Payne–Hanek reduction npyv__payne_hanek_reduce(x, *r_out) → quadrant。算法与 table data 移植自 Sun 的 FDLIBM e_rem_pio2.c / k_rem_pio2.c（public domain）：

1. 缩放 |x|，使其 biased exponent 变为 1046（unbiased 23），得到 [2^23, 2^24) 内的值。剥离三个 24-bit integer chunk tx[0..2]。
2. 索引预计算的 2/π bit table（66 个 chunk × 24 bit = 1584 bit，足以 覆盖 |x| < 2^1024 并留有余量）。加载 jx+jk+1 个 chunk，其中 jx = nx-1、jk = 4（FDLIBM 为 f64 精度采用的值）。
3. 将 tx[] 与加载的 chunk 相乘，得到 partial product q[0..jk]。 每次 f64 multiplication 都是精确的（24-bit × 24-bit 为 48 bit， 不超过 53 bit）。
4. 从高位向下把 q[] 分解成 24-bit integer piece iq[]，取出 x · 2/π 的 integer part，供 quadrant 使用。
5. 若剩余 fraction 大于 0.5，则 quadrant 加一，并对 iq[] 取 complement。 若 fraction 恰好缩为零则迭代（罕见；goto recompute）。
6. 把 fraction 乘回按 24-bit 对齐的 π/2 chunk，恢复 reduced argument。

• npyv_sincos.h——重构 kernel，将其拆为 reduction half （Cody–Waite 或 Payne–Hanek）和 polynomial half （npyv__sincos_f64_polynomial）。slow-path entry point npyv__sincos_f64_slow_path 逐 lane 执行 scalar Payne–Hanek，构造 (r, quadrant) vector，然后运行与 fast path 相同的 vectorized polynomial。NaN 与 ±inf input 直接输出 NaN，不经过 table-based reduction。

• loops_trigonometric.dispatch.cpp——per-block |x| < 2^20 gate 的 else branch 现在调用 npyv__sincos_f64_slow_path，不再循环调用 npy_sin/npy_cos。

精度审计（每个 range 10k 个 sample，与 scalar libm 对比）

具体检查：sin(2⁵⁰⁰) = 0.42925739234240... 与 libm 完全一致。 在 1e15、1e30、1e100、1e150、1e300 的 spot check 都与 libm bit-for-bit 一致。4699 个 umath test 全部通过。

Throughput（3A6000，n = 1M）

slow path 目前比 scalar libm 慢约 25%。polynomial half 已 vectorize，因此 一旦 per-lane reduction overhead 能被摊薄，差距应会缩小。下一项 follow-up 是 vectorize reduction（用 vector FP arithmetic 实现真正的 SIMD Payne–Hanek）；这需要对 per-lane table chunk load 进行 SIMD gather，而在 LASX 上需加入 xvinsgr2vr_d plumbing。

46 个 op×dtype pair 的 overall bench geomean（绝大多数为 fast-path input） 保持 4.94× 不变。

决策：scalar reduction + vector polynomial，而不是全 SIMD reduction

完全 vectorized 的 Payne–Hanek 需要逐 lane 从 66-entry table gather 三到四个 24-bit chunk。LASX 没有原生 gather instruction；替代方案是在四次 scalar load 后使用 xvinsgr2vr_d，与当前 scalar 方法相比并不真正节省时间。reduction 后的 polynomial 才是容易 vectorize 的部分，而这部分已经完成。

尝试 SIMD Payne–Hanek（2026-05-24）及失败原因

我做过一个 vectorized PH prototype：逐 lane gather 四个 2/π chunk，并用全 vector FP arithmetic 完成 multiplication 和 quadrant/fraction split。 它在 |x| ≤ ~10⁶ 时正确（1 ULP），但在 |x| ≳ 10¹⁰ 时出现灾难性错误 （3000+ ULP）。

根因是 chunk accumulation 中的精度损失。每个 chunk 的 contribution m · chunk_k · 2^(e − 52 − 24(k+1)) 跨越 77-bit range，但 f64 只能保留 53 bit。用普通 f64 arithmetic 将不同 bit scale 的 chunk 相加，会丢掉每次 sum 的低位：

• 对 x ≈ 10¹⁰，chunk 0 contribution 的 exponent 约为 +34，ULP 约为 2⁻¹⁹
• chunk 1 contribution 的 exponent 约为 +10，相加时会对齐到同一 ULP scale，因此损失约 15 bit
• chunk 2、3 的损失依次更大

求和后，fractional part 只剩约 20 bit，远不足以满足 sin/cos 所需的 53 ULP 精度。

修复方法都不简单：

1. Double-Double accumulation。 用 ddadd2（使用 FMA correction 的 Dekker-style two-sum）替换 partial product 之间的普通 +。每一步 accumulator 都保留结果的 high 与 low part。每次 add 每个 element 约多 6 个操作 × 4 次 add = 24 个操作，multiplication 也需要 DD 形式。
2. Multi-precision integer accumulation。 把 running sum 保持为多个独立的 24-bit-piece vector（FDLIBM 的 iq[0..jz] 风格），并在 vector side propagate carry。实质上就是把 FDLIBM 的 __kernel_rem_pio2 连同 per-lane state 直接移植到 SIMD。
3. 预对齐 table（SLEEF 方法）。生成 table 时就把每个 entry 的四个 double 按直接 multiplication 所需 scale 对齐。需要用 mpmath 生成 table；constant table 会从 66 个 int 增长到数千个 double。

（1）最清晰；（3）是 SLEEF 实际采用的方法。二者都要增加 100–300 行谨慎编写的 代码。当前 scalar PH 工作正常，能以 7–9 M/s 在完整 f64 range 内输出正确结果， 而且 bench 不受影响（bench input 全部在 fast path）。因此 SIMD-reduction follow-up 继续保留，但不阻塞发布。

决策：不为 arcsinh f64 增加直接 Padé

arcsinh 现有的 log1p(|x| + x²/(1+√(1+x²))) 会在 bench range |x| < 2 内产生通常大于 0.5 的 log1p argument（例如 arcsinh(0.5) → log1p(0.618)，arcsinh(1) → log1p(1.414)）。因此， log1p 中 threshold 为 0.5 的 per-block dispatch 在 arcsinh call path 中很少 触发。arcsinh 的直接 polynomial 需要自己的 reduction（FDLIBM 按 |x| 分为四个 range），工程投入较大，而在这颗 CPU 上与 bench 有关的收益有限。如果真实 workload 需要，再作为后续工作处理。

2026-05-24 — 从 libm 收回 tanh / cbrt-f32 / power

文件： numpy/_core/src/umath/loops_umath_fp.dispatch.c.src、 numpy/_core/meson.build

背景

在 Intel 上与 stock NumPy 进行的 cross-CPU bench 表明，有三个 transcendental 落后幅度远大于 f64/f32 SIMD-width gap 所能解释：tanh 为 0.06×、cbrt f32 为 0.06×，arcsinh/power 在 0.08–0.11×。诊断如下：

• tanh 在 loops_hyperbolic.dispatch.cpp.src 中有列入 LASX/LSX 的 Highway-based SVML port，但 bench 数字强烈暗示它在 Loongson 上落回了 scalar libm（f64 为 54 M/s，约等于每 element 18 ns 的 scalar tanh）。
• cbrt f32 仍采用缓慢的 exp(log(x)/3) 组合，而 cbrt f64 已在 2026-05-23 改写为 FDLIBM bit-magic；f32 当时被遗漏。
• 即使没有 lane 需要处理 negative base、±inf y、±0、±1 等 IEEE-754 edge case， power 仍为每个 element 承担约 50 个 select 操作。

arcsinh 是唯一的 outlier——见上一节。log1p 形式对其精度至关重要，如果不使用 自定义推导的 Padé coefficient，很难缩短 path。

实现

• tanh f32/f64。 新 NPYV kernel 使用 tanh(x) = sign(x) · (1 − 2/(exp(2|x|) + 1))——只需一次 exp 调用， 而 sinh/cosh decomposition 需要两次。exp 前把 |x| clamp 到 9（f32）/ 20（f64），使 e^{2|x|} 保持 finite；超过这些 threshold 后，tanh 在 ULP 精度下已饱和为 ±1。

从 meson.build 的 loops_hyperbolic arch list 中移除 LASX, LSX，使该 source 不再为 LoongArch build 生成 dispatch function，避免 duplicate-symbol collision。loops_umath_fp 中新的 DOUBLE_tanh / FLOAT_tanh 定义受 NPYV_HAVE_tanh_@sfx@ 控制（仅当 NPYV exp/log path 编译时为 true，也就是 LSX/LASX），因此其他 target 仍以 loops_hyperbolic 的 SVML port 为准。

• cbrt f32。 替换为 FDLIBM cbrtf 算法：bit-magic initial estimate （hx/3 + B1，其中 B1 = 709958130），再执行两次 f32 Halley iteration。 Halley 为 cubic convergence（5 → 15 → 45 bit），因此两次 iteration 足以覆盖 24-bit mantissa 并留有余量。subnormal 先乘以 2^24，并使用 B2 = B1 − 8·2^23 = 642849266 补偿 exponent field。完全跳过 log + exp pair。

• power f64/f32 hot path。 在现有 kernel 外增加 per-block fast check：若 每个 lane 均满足 x > 0、x 和 y 都是 finite、y ≠ 0 且 x ≠ 1，则直接 result = exp(y · log(x))，不做 select wave。大多数 workload（统计、模拟、 图像处理）都满足条件，此前却仍要每个 element 支付约 50 个 select。slow-path code 保持原样，negative base、infinity、±0 等情况仍能正确处理。

权衡

• tanh 在 |x|=9（f32）处饱和，对典型 workload 在 ULP 精度下正确，但当 x 为 upper bit 中 payload 不同的 NaN 时会丢失 sign information。kernel 末尾会把 NaN 原样传递，因此并无实际影响。
• cbrt f32 的 bit-magic 依赖 npyv__cvt_f32_s32 可用，这与现有 f64 implementation 的约束一致。LSX 与 LASX 都提供该操作。
• pow fast path 每个 block 要花约 5 个操作构造 mask。对全部为 clean data 的 common case 可忽略；worst case（每个 block 都有 special value）则同时支付检查 和 slow path。common case 的收益足以覆盖这个检查。

已在 Intel 验证

四个 op 都能正确编译。在 Intel 上，tanh 仍走 loops_hyperbolic SVML port （只从其 arch list 移除了 LASX/LSX）；f32 cbrt 和 pow 不会使用我们的 NPYV path（那里未定义 NPYV_IMPL_F32_EXP_LOG）。对 libm 的验证等待 Loongson rebuild。

Build 陷阱：baseline duplicate symbol

第一次 Loongson rebuild 因 loops_umath_fp baseline build 与 loops_hyperbolic baseline build 都定义 DOUBLE_tanh 而失败。Loongson meson_cpu config 把 LSX 放在 baseline，因此即使对每个 dispatch source 做 baseline build，也会定义 NPY_HAVE_LSX。结果是两个 source 的 baseline 中 NPYV_HAVE_tanh_* 都被定义，双方都生成了该 symbol。

修复方法是用 #if !defined(NPY_HAVE_LSX) && !defined(NPY_HAVE_LASX) guard loops_hyperbolic 的 tanh 定义。loops_umath_fp 一侧已正确受 NPYV_HAVE_tanh_@sfx@ 控制，因此 collision 现在只剩单一来源：LoongArch 上 symbol 来自 loops_umath_fp；其他平台来自 loops_hyperbolic。

Bench 意外：tanh f32——Highway port 更快，继续保留

rebuild 后第一次 end-to-end bench 显示，tanh f32 从 307 → 166 M/s，反而出现 regression。Highway-based SVML port 原本就在 LASX 上正常运行（与我先前的假设 不同）；我们的 NPYV expm1-based 形式在 f32 上更慢，因为 small-x polynomial path 的 per-block dispatch 会在典型 bench input distribution [-3, 3] 上 误触发：每个 lane 约有 8% 的 input 小于 0.25，导致约 50% 的 8-lane LASX block 进入 “mixed” branch，同时支付两条 path 的成本。

修复方法是完全不编译 f32 NPYV kernel。从 loops_umath_fp.dispatch.c.src 移除 #define NPYV_HAVE_tanh_f32 1，因此其中的 dispatch wrapper 只为 f64 生成。loops_hyperbolic source 的 tanh #if !@loong_skip@ gate 现在使用 per-iteration loong_skip = 0, 1 repeat variable，只在 LoongArch 上跳过 DOUBLE。LASX/LSX 的 FLOAT_tanh 继续来自 Highway SVML port。

Loongson 的 DOUBLE_tanh 确实改用 NPYV：Highway path 在那里对 f64 fallback 到 scalar libm（没有 LASX f64 Highway gather），改动前为 54 M/s，使用我们的 NPYV 形式后为 81 M/s，f64 净提升 1.49×。

最终 Loongson 3A6000 数字

验证：各 100k 个 random input 中，tanh f32 ≤1 ULP、tanh f64 ≤3 ULP、 cbrt f32 ≤3 ULP。power 在 result 约为 1e-18 的 extreme value（例如 pow(99, -8.6)）处会达到 ≤65 ULP——这是 f64 中 exp(y · log(x)) 不可避免的 ULP 放大：约 1 ULP 的 log error × |y| ≈ 8.6，会在 y·log(x) 中放大为约 9 ULP；随后 exp 把 absolute precision 指数化为 tiny value 的 relative precision，换算后约为 65 ULP。接近 1.0 的 pow result 保持 ≤8 ULP。 已记录为已知限制；若要修复，需要在整个 log+exp pair 上使用 double-double arithmetic。

2026-05-25 — Release 1：日常基础操作启用 LASX，以及两项权衡

为此前 transcendental 工作仍留在 LSX 的日常 element-wise dispatch unit 启用 LASX（256-bit）：loops_arithm_fp（float add/sub/mul/divide + complex）、 loops_unary、loops_unary_fp（sqrt、reciprocal、square、abs、 rint/floor/ceil/trunc）、loops_unary_complex 与 loops_autovec。这几乎不需 额外实现：LASX NPYV backend（NumPy 自己的 upstream/community work）已经实现这些 文件需要的全部 intrinsic，只是 meson.build 中的开关尚未打开。

有两个 dispatch unit 刻意保留在 LSX，因为对 store-bound、byte-output op， LASX 的 wide→byte pack 需要跨 128-bit lane 的 xvpermi_d fixup，整体反而变慢：

• loops_unary_fp_le（isnan/isinf/isfinite/signbit）——upstream 已通过 NPY_SIMD_FORCE_128 强制使用 128-bit。
• loops_comparison（less/greater/…）——实测 LASX 相对 LSX，less uint32 为 0.78×、uint64 为 0.87×，且没有任何 comparison dtype 快于 1.02×。因此从其 arch list 移除 LASX。minmax 保持 LASX（output 同宽，无 pack；max 提升 1.3–1.5×）。

保留 LASX 的已知权衡：f32 add/sub/mul

loops_arithm_fp 是单个 dispatch unit，只能全开或全关。在 3A6000 上实测： f32 divide +92%，全部 f64 arithmetic +9–14%，但 f32 add/sub/mul −19%。为保留 divide/f64 的收益，仍采用 LASX。尚未 profile f32 regression 的原因——可能受 memory bandwidth 限制，但这只是 hypothesis，并非测量结论； cache-resident 与更高 bandwidth 部件上的表现仍未测试。详见 BENCHMARKS.md。

Headline（Release 1）

完整 rebench（含 sort）：dragon 2.5.0.dev0+dragon.unofficial.1 对 stock NumPy 2.5.0.dev0；在空闲 3A6000 上 back-to-back 运行，--iters 30，共 135 个 op × dtype pair：geomean 3.31×（f16 10.7×、bool 11.4×、f32 4.5×、 f64 3.6×、int 1.8×）。

Bench 方法说明

只比较在空闲机器上 back-to-back 获取的 run。曾经一次用旧 stock run 对比刚完成的 dragon run，产生了“修复反而变慢”的假象；max/multiply control op（本应稳定） 是识别 contaminated run 的信号。