np.add(a, b)

の旅：NumPy 内部構造の探求

NumPy の開発に携わり、「配列を加算する」というシンプルな操作が、Python の呼び出しから SIMD カーネルの実行まで、どのような仕組みで処理されているのかを追跡しました。

本記事では、NumPy 2.5.2 を基準として、

np.add(a, b)

注意: 本記事は NumPy 2.5.2 の内部構造に基づいています。バージョン間での変更がある場合は、該当するタグで独自に検証（spelunking）を行ってください。C の読解は必須ですが、事前の NumPy 知識はありません。

地図：呼び出しの階層

追跡の全体像を表した「宝地図」です。上位から下位へ分解されます。

• np.add(a, b)

  • (Python: エントリーポイント)

• ▼

• ufunc_generic_fastcall

  • (C: 引数の解析・初期化)

• ▼

• __array_ufunc__ オーバーライドチェック

  • (サードパーティライブラリへの分岐判定)

• ▼

• プロモーションとディスパッチ

  • (float64 ループの探索、キャッシュ化)

• ▼

• 単純なループ

  または

  NpyIter

  • (反復戦略：高速パスかイテレーター生成か)

• ▼

• DOUBLE_add

  • (実際の内部ループ、SIMD 対応)

1.

np.add

はオブジェクトである

np.add

numpy.ufunc

>>> type(np.add)
<class 'numpy.ufunc'>
>>> np.add.nin, np.add.nout
(2, 1)          # 入力数:2, 出力数:1
>>> len(np.add.types)
22              # 登録された型シグネチャの数
>>> np.add.types[11:14]
['ee->e', 'ff->f', 'dd->d']

• ufunc（通用関数）の本質:
  • 対応する型シグネチャごとに 1 つの C 関数。
  • メタデータ（入出力数など）を保持。
• 追跡対象:

  dd->d

  (double 入力 → double 出力)。

2. C への潜入：ufunc_generic_fastcall

Python が

np.add(a, b)

ufunc_generic_fastcall

static PyObject *
ufunc_generic_fastcall(PyUFuncObject *ufunc, ...) {
    // 1. 入力・出力・キーワード引数を抽出
    
    // 2. オーバーライドチェック (`__array_ufunc__` など)
    //    サードパーティライブラリ（Dask, CuPy）のフック判定
    
    // 3. 引数→配列に変換、DTyes（データ型）を抽出
    
    // 4. ディスパッチとプロモーションの実行
    PyArrayMethodObject *ufuncimpl = promote_and_get_ufuncimpl(...);
    
    // 5. デスクリプタの解決
    
    // 6. 最終処理：内部ループへの実行
    errval = PyUFunc_GenericFunctionInternal(ufunc, ufuncimpl, ...);
    
    // 7. 結果をラップして返す
}

3. オーバーライドチェックとサードパーティの介入

NumPy は実際に計算を行う前に、引数が**「自分で行うべきこと」**がないか確認します。

• __array_ufunc__

  プロトコル:
  • NEP 133 で定義されたメソッド。
  • Dask や CuPy などのライブラリが主導権を握るためのフック。
• 動作例:

  class Diverted:
      def __array_ufunc__(self, ufunc, method, *inputs, **kwargs):
          return f"intercepted {ufunc.__name__}.{method}"
  
  >>> np.add(np.arange(3), Diverted())
  'intercepted add.__call__'
  

通常の配列（ndarray）の場合、このチェックを通過して計算権限を維持します。

4. ループの選択：プロモーションとディスパッチ

2 つの

float64

dispatching.cpp

ディスパッチステップ

1. オーバーライド適用:

   signature

   の固定があるか確認。
2. キャッシュチェック: 同じ組み合わせで呼ばれた場合はキャッシュヒット。
3. マルチディスパッチ: 最適な「ループ」を探す（サブクラスを持つ DType）。
4. プロモーター: マッチしない場合、型を変換して再試行するヘルパー。
5. 登録確認:

   ArrayMethod

   が発見され、ufunc ループに登録される。

クラシックなループへのアダプター

新しいインターフェース (

PyArrayMethodObject

PyUFuncGenericFunction

static int
generic_wrapped_legacy_loop(PyArrayMethod_Context *context, ...) {
    legacy_array_method_auxdata *ldata = (legacy_array_method_auxdata *)auxdata;
    
    // クラシックなループインターフェースへの呼び出し
    ldata->loop(data, dimensions, strides, user_data); 
    
    return 0;
}

5. 反復するか？：高速パスと NpyIter

実際に計算を行う前に、「単純なループ（trivial loop）」を使うべきかを判断します。

高速パス (

try_trivial_single_output_loop

)

条件を満たせば、イテレーター (

NpyIter

• 条件:
  • 形状が一致し、ブロードキャスト不要。
  • キャスト不要。
  • オペランドがすべて 1-D または連続（contiguous）。
• 動作: 単一の C 関数呼び出しで完了。

一般パス (

execute_ufunc_loop

+

NpyIter

)

高速パスの条件に合わなければ、一般的な配列イテレーターを使用します。

npy_uint32 iter_flags = ufunc->iter_flags |
             NPY_ITER_EXTERNAL_LOOP | // イテレーター管理
             NPY_ITER_REFS_OK | 
             NPY_ITER_BUFFERED |      // ブロードキャスト対応など
             ...;
             
// メインループ
do {
    res = strided_loop(context, dataptr, countptr, strides, auxdata);
} while (res == 0 && iternext(iter)); // 次のスライスへ進む

• GIL の解放: ループ中は GIL を解放し、マルチスレッド/プロセス環境で並行性を確保。
• 浮動小数点フラグ: 計算前後にクリアされ、オーバーフローなどを検出。

ストライドとメモリの関係

ストライド（stride）がある場合でも、1-D であれば高速パスが可能です。ただし、不連続な N-D ビューや

out

ベンチマーク例:

import numpy as np
import timeit

n = 10_000_000
a, b, out = (np.random.rand(2 * n) for _ in range(3))
contig_a, contig_b, contig_out = (x[:n].copy() for x in (a, b, out))
strided_a, strided_b, strided_out = (x[::2] for x in (a, b, out))

time_contig = timeit.timeit("np.add(contig_a, contig_b)", globals=globals(), number=10)
time_strided = timeit.timeit("np.add(strided_a, strided_b)", globals=globals(), number=10)

• 結果: 連続配列が約 2.6 倍速い。
• 理由: ストライド版はメモリ帯域幅を 2 倍消費するためだが、根本的な差は

  DOUBLE_add

  の内部実装にある。

6. ループそれ自身：SIMD と特殊化

DOUBLE_add

loops_arithm_fp.dispatch.c.src

コード生成の仕組み

/**begin repeat
 * Float types: #type = npy_double#
 * Arithmetic: # kind = add#
 */
NPY_NO_EXPORT void NPY_CPU_DISPATCH_CURFX(DOUBLE_add)
(char **args, ...) {
    // 内部実装は CPU 依存で変化する
}

• conv_template.py

  がテンプレートを展開し、複数の C 関数（

  FLOAT_add

  ,

  DOUBLE_add

  など）を生成。

特殊化の階層

1. 還元ケース: 結果への書き込みが必要ない場合。
2. SIMD 加速版:
   • 条件: 両方の入力が連続、または片方がスカラー。
   • 実装:

     npyv_add_f64

     （CPU 固有の SIMD 命令：AVX, NEON など）。
3. フォールバック (Scalar):
   • シグナル化されない場合の単純なループ。

ストライド呼び出しの性能低下理由:

• SIMD 特殊化は「オペランドが連続」であることを条件としている。
• ストライド配列はこれを満たさないので、**スカラーブランチ（単純な for ループ）**にフォールバックする。
• これにより、SIMD の恩恵を受けず、メモリ帯域幅も効率的でないため低速になる。

7. ループはどこから来るか：動的生成と CPU チェック

生成された関数 (

DOUBLE_add

ufunc->functions

ディスパッチテーブルの構成

• 生成元:

  generate_umath.py

  (Python) →

  __umath_generated.c

  (C)。
• 構造例:

  'add': {
      'dd->d': {
          dispatch=[('loops_arithm_fp', 'fdFD')], # 浮動小数点/複素数系
      },
  }
  

• ビルド時: マクロで CPU 機能をチェックし、最適なループを選択。

  NPY_CPU_DISPATCH_CALL(func = add);
  // X86: AVX512 -> AVX2 -> Baseline
  // ARM : NEON (条件付き)
  

• 確認方法:

  >>> from numpy.lib import introspect
  >>> introspect.opt_func_info("add", "float64")
  # CPU 依存の結果を表示（例: 'baseline(NEON ...)'）
  

8. 結論：地質学としての NumPy

np.add(a, b)

1. オブジェクト解析:

   ufunc

   オブジェクトの確認とオーバーライドチェック。
2. 型解決: データ型 (

   float64

   ) を特定し、キャッシュされたループを探索。
3. 戦略選択: 配列の性質（連続性など）から高速パスかイテレーター生成かを選択。
4. 実行: GIL 解放を行い、CPU の能力に応じた SIMD ループを実行。

技術的インサイト

• レガシーの層:

  legacy_array_method.c

  など古いインターフェースがラップされながらも、今も主力として稼働している。
• 進化の痕跡: 型システムや SIMD 対応など、数十年前のコードの上に新しいレイヤーが積み重ねられ、かつ内部でも継続的に改良されている（例：ハッシュテーブル化されたループ選択）。
• 動機: この探求は、NumPy を「スレッドフリー Python」への移行などの大きな変化を予感させる。

おわりに: 単純な加算操作が、これほど多くの設計と最適化を内包していることは驚異的です。より深い理解のために、ブレークポイントを置いたソースコード追跡をお試しください。