CPythonのソースコードを読んでみた

estieでソフトウェアエンジニアをしている安東です。普段の業務で関わっているのはPython製のシステムが中心ですが、過去データを分析するのにちょっとRustを使ってみたりもしています。

こうやって普段からお世話になっているPythonですが、中身がわからないまま使い続けることに対してやや不安を感じることがあります。たしかに書き方だけ知っていれば大抵の場面でなんとかなってしまうのでしょうが、それだけではカバーしきれないところがどんなプロジェクトでも突然やってきます。カバーできる範囲を増やすためにも、突然の出来事の予兆を事前に嗅ぎつけるためにも、プログラムが動いている感覚が地に足ついた形でほしいのです。

ということで、今回はPythonの処理系でおそらく一番メジャーなCPythonのソースコードを読んでみようと思います。ただ、ソースコード全体を読むには時間も記事のスペースも足りないので、今回は視点を少し変えて、「初めて読むソースコードから目的の実装本体を見つけるにはどうするか」実況のような形でお届けします。

レギュレーション

筆者はPythonはよく使っていますが、CPythonを実装を読むのは初めてです。
言語処理系に関する一般的な知識は少しあります。
「下記のコードの2つのメソッドalloc_list とalloc_bytes の速度差の原因らしいところを見つける」ことを目標に読んでいきます。

def alloc_list(size):
    return [0] * size

def alloc_bytes(size):
    return b"A" * size

alloc_list とalloc_bytes の中身は非常に単純で、長さsizeで同じ値が並んだ領域を返すだけになっています。

それでは捜索スタート！

第1ステージ：再現性を確認せよ！

事象が手元で再現しなかったり、そもそもバイナリが動いていなかったりビルドが通っていなかったりするような状況で調査を進めるのはとてもつらいところがあります。そこで、まずは執筆時点でのCPythonの最新のソースコード (Python 3.12 α 0) を落としてきて、動作チェックと事象の再現確認を行います。

CPythonのビルドは、大筋Githubのreadmeと公式ガイドに従うことでM1 Macでもこのようにすんなりできました。ここでは公式ガイドと異なりmakeを4並列で動かすようにしていますが、目立った問題はなく、2分ほどでmakeが完了します。C言語のプロジェクトなのに一発でビルドが通るとは感動です。

$ git clone https://github.com/python/cpython
$ cd cpython
$ ./configure --enable-optimizations
$ brew install pkg-config openssl@1.1 xz gdbm tcl-tk
$ PKG_CONFIG_PATH="$(brew --prefix tcl-tk)/lib/pkgconfig" \
  ./configure --with-openssl=$(brew --prefix openssl@1.1) --enable-optimizations
$ make -s -j4
$ ./python.exe
>>> print("Hello")
Hello

python.exeも立ち上がり、ちゃんと動作していそうです。100回同じ操作を繰り返してパフォーマンス計測するデコレータをサクッと書いて、レギュレーションで挙げた2つのメソッドalloc_list とalloc_bytes の性能差がどれほどか、要素数を128M個（約1億個）にして実験してみます。

パフォーマンス計測用のデコレータの実装（bench.py）

# bench.py
import time
import math
from functools import wraps

def bench(num_try=100):
    def bench_impl(func):
        @wraps(func)
        def wrapper(*args, **kwargs):
            samples = []
            for i in range(num_try):
                print(f"{i+1}/{num_try}",end="\r")
                start = time.time_ns()
                func(*args, **kwargs)
                end = time.time_ns()
                samples.append(end - start)
            min_ns = min(samples)
            max_ns = max(samples)
            avg_ns = sum(samples) / num_try
            std_ns = math.sqrt(sum(item**2 for item in samples) / num_try - avg_ns**2)
            min_ms = min_ns / 1000000
            max_ms = max_ns / 1000000
            avg_ms = avg_ns / 1000000
            std_ms = std_ns / 1000000
            print(f"min/avg/max/stddev = {min_ms:.3f}/{avg_ms:.3f}/{max_ms:.3f}/{std_ms:.3f} ms")
        return wrapper
    return bench_impl

# list.py
from bench import bench

@bench()
def alloc_list(size):
    return [0] * size

if __name__ == "__main__":
    alloc_list(128 * 1024**2)

# bytes.py
from bench import bench

@bench()
def alloc_bytes(size):
    return b'A' * size

if __name__ == "__main__":
    alloc_bytes(128 * 1024**2)

$ ./python.exe ~/work/list.py
min/avg/max/stddev = 474.193/494.015/528.137/10.762 ms

$ ./python.exe ~/work/bytes.py 
min/avg/max/stddev = 14.424/18.388/26.094/1.688 ms

alloc_list のほうが約27倍遅いという結果になりました。 alloc_list とalloc_bytes どちらも単純で似た中身のコードですが、かなり大きな性能差です。これだけ性能差があるなら、性能差があるということを頭の片隅に置いておくだけでも、いつか役に立ちそうな雰囲気があります。

ということで、これで調査に使うソースコードと動作確認用のバイナリが揃いました。さらに調査対象の「alloc_list とalloc_bytes とで動作速度が違う」も再現確認できました。第1ステージクリアです！

第2ステージ：取っ掛かりを探せ！

ここからはPythonに関する知識その他諸々を動員して、Visual Studio Codeでソースコードを追いかけていく起点になる場所を探します。C言語レベルだとlistに対する* とbytesに対する* とで走る処理が異なるはずなので、それぞれの処理の起点になる場所が目標です。

Pythonの* 演算子に対応するC言語の実装が簡単に見つかればよいが、記号1文字であまりに検索性が低い。

ところで、Pythonでモジュールをimportすると*.pycというキャッシュファイルができることを知っている。このファイルの目的からして、python.exeの中ではソースコードはなにかしらの中間表現に変換して持っているはず。

ということは、中間表現の中身と、* 演算子と中間表現の対応がわかれば中間表現のほうから追っていけるかも。

と考えてGoogle検索すると、標準で入っているdisモジュールを使うとPythonのメソッドをバイトコードとして読めることがわかりました。中の人くらいしか使わなさそうなモジュールにまでしっかりとしたドキュメント（しかもほぼ日本語！）があるのはPythonのとても嬉しいところに感じます。なお、*.pycの中身とバイトコードの関係は調べても結局よくわかりませんでした。

このdisモジュールを使ってalloc_list とalloc_bytes の中身を出力してみるとこうなりました。disモジュールのリファレンスの冒頭にあるとおり、この出力結果はPythonのバージョンによって異なることがあるので注意が必要です。

>>> import dis
>>> def alloc_list(size):
        return [0] * size
>>> def alloc_bytes(size):
        return b'A' * size
>>> dis.dis(alloc_list)
  1           0 RESUME                   0

  2           2 LOAD_CONST               1 (0)
              4 BUILD_LIST               1
              6 LOAD_FAST                0 (size)
              8 BINARY_OP                5 (*)
             12 RETURN_VALUE
>>> dis.dis(alloc_bytes)
  1           0 RESUME                   0

  2           2 LOAD_CONST               1 (b'A')
              4 LOAD_FAST                0 (size)
              6 BINARY_OP                5 (*)
             10 RETURN_VALUE

各命令の意味の細かいところはわからないものの、どこがボトルネックになっていそうか一応見てみます。

LOAD_CONST , BUILD_LIST , LOAD_FAST はいずれもここではごく小さいデータしか扱っていないと思われ、実行時間は無視していいはずです。RETURN_VALUE も参照を返しているだけのはずなので無視してよさそうです。

残るはBINARY_OP ですが、これが* 演算子に相当する命令（もう少し言うと、オペランドによって演算を切り替える二項演算）で、時間のかかっているところと思われます。さらにこの語、小文字が含まれていないので検索で引っかかるノイズが少なくてとてもよさそうです。

ということでここでようやくC言語のソースコードをVisual Studio Codeで開き、BINARY_OP で検索します。検索ノイズはまだまだあるのですが、ファイル名だけ見るとceval.cが怪しく見えます。このファイルの中でBINARY_OP の引っかかる適当な箇所の周辺を調べると、TARGET(hoge) が case hoge: （後略） に展開されることがわかりました。バイトコードの命令を解釈して実行する部分は巨大なswitch-case文になっていると予想してTARGET(BINARY_OP) で調べると、5592行目にありました。このcaseの中身を見るといかにもオペランドによって演算を切り替えていそうな処理（元ソースの5599行目・下記引用の11行目）があるので、この部分がBINARY_OP 命令の実装そのもののようです。

// CPythonでBINARY_OP命令を処理している箇所（コメントは筆者によるもの）
// Python/ceval.c の5592行目以降
        TARGET(BINARY_OP) {     // case BINARY_OP: ...
            // 引数の準備とassert
            PREDICTED(BINARY_OP);
            PyObject *rhs = POP();
            PyObject *lhs = TOP();
            assert(0 <= oparg);
            assert((unsigned)oparg < Py_ARRAY_LENGTH(binary_ops));
            assert(binary_ops[oparg]);
            // BINARY_OPの本体
            PyObject *res = binary_ops[oparg](lhs, rhs);
            // 戻り値のセットと次の命令へのジャンプ
            Py_DECREF(lhs);
            Py_DECREF(rhs);
            SET_TOP(res);
            if (res == NULL) {
                goto error;
            }
            JUMPBY(INLINE_CACHE_ENTRIES_BINARY_OP);
            DISPATCH();
        }

binary_opsがどこで設定されているのかまだよくわかりませんが、それを追いかけていけば何かつかめそうです。

ひとまず、listに対する* とbytesに対する* 、それぞれの処理の起点になる場所が見つかったので、これで第2ステージクリアということにしておきます。もしかすると次の最終ステージで行き詰まるかもしれませんが、そのときはこのステージでおこなった推測のうち確度の低いところを見直してやり直せばいいです。

最終ステージ：本丸に向かえ！

ここからはVisual Studio Codeの力を使ってコードを追いかけていきます。

ひとまず上で引用した5592行目からのコードは、5599行目以外は名前からしてお釣りのようなコードなのでひとまず忘れることにします。5599行目で呼ばれる実体が知りたいのでVisual Studio Codeの機能でbinary_ops の定義に移動すると、* 演算子を実行するときはPyNumber_Multiply 関数が呼ばれることがわかります。

今知りたいのは数に対する演算ではないので一瞬ドキッとするのですが、さらにabstract.cの1101行目にある定義に飛んで中身を見ると、数以外に対する実装も書かれていて安心できます。直上にあるsequence_repeat 関数の実装もあわせて見ると、PyNumber_Multiply 関数の中身はほとんどPySequenceMethod 構造体にあるsq_repeat の指す先を呼んでいるだけとわかります。

問題はsq_repeat の指す先が何かなのですが、まずこのメンバをFind all referencesで探しても値を設定していそうな箇所が見つかりません（注：執筆する際に改めて確認したところ、OTHER REFERENCES RESULTSにちゃんと載っていました。初見では気づきませんでした）。次に構造体の型名 PySequenceMethods で調べてもよくわかりません。そしてsq_repeat の型名ssizeargfunc で調べると、listobject.cやbytesobject.cで、list_as_sequence やbyte_as_sequence を初期化している箇所が見つかりました。定義元に移動すると、どちらの関数も求めていたlistやbytesに対する* 演算子の実装本体に見えます。

// 同じ要素を繰り返したlistを生成する処理のCPython内部での実装（コメントは筆者によるもの）
// Objects/listobject.c の551行目以降
static PyObject *
list_repeat(PyListObject *a, Py_ssize_t n)
{
    // 空のlistは何回繰り返しても空のまま
    const Py_ssize_t input_size = Py_SIZE(a);
    if (input_size == 0 || n <= 0)
        return PyList_New(0);
    assert(n > 0);

    // メモリ確保
    if (input_size > PY_SSIZE_T_MAX / n)
        return PyErr_NoMemory();
    Py_ssize_t output_size = input_size * n;

    PyListObject *np = (PyListObject *) list_new_prealloc(output_size);
    if (np == NULL)
        return NULL;

    // 同じ要素を繰り返したlistを生成する処理の本体
    // 特に input_size == 1 の場合が最初に速度を比較した対象
    PyObject **dest = np->ob_item;
    if (input_size == 1) {  // (a-1) コピー元のlistの要素数が1つだけの場合
        PyObject *elem = a->ob_item[0];
        _Py_RefcntAdd(elem, n);
        PyObject **dest_end = dest + output_size;
        while (dest < dest_end) {
            *dest++ = elem;
        }
    }
    else {  // (a-2) コピー元のlistの要素が複数の場合
        PyObject **src = a->ob_item;
        PyObject **src_end = src + input_size;
        while (src < src_end) {
            _Py_RefcntAdd(*src, n);
            *dest++ = *src++;
        }

        _Py_memory_repeat((char *)np->ob_item, sizeof(PyObject *)*output_size,
                                        sizeof(PyObject *)*input_size);
    }

    Py_SET_SIZE(np, output_size);
    return (PyObject *) np;
}

// Include/internal/pycore_list.h の60行目以降
static inline void
_Py_memory_repeat(char* dest, Py_ssize_t len_dest, Py_ssize_t len_src)
{
    assert(len_src > 0);
    Py_ssize_t copied = len_src;
    while (copied < len_dest) {
        Py_ssize_t bytes_to_copy = Py_MIN(copied, len_dest - copied);
        memcpy(dest + copied, dest, bytes_to_copy);
        copied += bytes_to_copy;
    }
}

// 同じ要素を繰り返したbytesを生成する処理のCPython内部での実装（1行コメントは筆者によるもの）
// Objects/bytesobject.c の1443行目以降
static PyObject *
bytes_repeat(PyBytesObject *a, Py_ssize_t n)
{
    // 以下しばらく PyBytesObject のメモリ領域の確保
    Py_ssize_t size;
    PyBytesObject *op;
    size_t nbytes;
    if (n < 0)
        n = 0;
    /* watch out for overflows:  the size can overflow int,
     * and the # of bytes needed can overflow size_t
     */
    if (n > 0 && Py_SIZE(a) > PY_SSIZE_T_MAX / n) {
        PyErr_SetString(PyExc_OverflowError,
            "repeated bytes are too long");
        return NULL;
    }
    size = Py_SIZE(a) * n;
    if (size == Py_SIZE(a) && PyBytes_CheckExact(a)) {
        Py_INCREF(a);
        return (PyObject *)a;
    }
    nbytes = (size_t)size;
    if (nbytes + PyBytesObject_SIZE <= nbytes) {
        PyErr_SetString(PyExc_OverflowError,
            "repeated bytes are too long");
        return NULL;
    }
    op = (PyBytesObject *)PyObject_Malloc(PyBytesObject_SIZE + nbytes);
    if (op == NULL) {
        return PyErr_NoMemory();
    }
    _PyObject_InitVar((PyVarObject*)op, &PyBytes_Type, size);
_Py_COMP_DIAG_PUSH
_Py_COMP_DIAG_IGNORE_DEPR_DECLS
    op->ob_shash = -1;
_Py_COMP_DIAG_POP
    op->ob_sval[size] = '\0';

    // 同じ要素を繰り返したbytesを生成する処理の本体の呼び出し
    _PyBytes_Repeat(op->ob_sval, size, a->ob_sval, Py_SIZE(a));

    return (PyObject *) op;
}

// Objects/bytesobject.c の3561行目以降
void
_PyBytes_Repeat(char* dest, Py_ssize_t len_dest,
    const char* src, Py_ssize_t len_src)
{
    // 空のbytesは何回繰り返しても空のまま
    if (len_dest == 0) {
        return;
    }
    if (len_src == 1) { // (b-1) コピー元のbytesの要素数が1の場合
        memset(dest, src[0], len_dest);
    }
    else {  // (b-2) コピー元のbytesの要素が複数の場合
        if (src != dest) {
            memcpy(dest, src, len_src);
        }
        Py_ssize_t copied = len_src;
        while (copied < len_dest) {
            Py_ssize_t bytes_to_copy = Py_MIN(copied, len_dest - copied);
            memcpy(dest + copied, dest, bytes_to_copy);
            copied += bytes_to_copy;
        }
    }
}

これでlistとbytesどちらも実装本体の場所が突き止められたので、第2ステージに戻ることなくめでたく最終ステージもクリアです！

Exステージ：速度の違いを解明せよ！

ここまででCPythonのソースコードを読んでみつつ実装箇所を探し出すという目的は一応達成できました。ただ、題材の「性能差はどこから来るのか？」には答えられていないので、まだ歯切れが悪い感じです。最終ステージと同様に追っていって見つけられたのは次の2点でした。

中で確保されるデータのサイズが違う
- list: 1要素あたり8byte (sizeof(PyObject*)) ※ 同じポインタが使い回されるので、細かいアロケーションが頻発することはない
- bytes: 1要素あたり1byte (sizeof(char))
コピーの仕方が違う
- list
  - (a-1) 元オブジェクトの要素数が1の場合：1要素ずつコピー
  - (a-2) 元オブジェクトの要素が複数の場合：コピーする長さを倍々にしながらmemcpy
- bytes
  - (b-1) 元オブジェクトの要素数1の場合：memsetだけ
  - (b-2) 元オブジェクトの要素が複数の場合：コピーする長さを倍々にしながらmemcpy

なんと元オブジェクトの要素数によってコピーの仕方がが変わるんですね。それぞれの実装に対応する場所は最終ステージで引用したコードにコメントで示してあります。一見するとGCのための参照カウントの管理コストのぶんだけlistのほうが不利に見えますが、参照カウントはデータのコピー前に一気に増やしているのでほとんど関係しないと思われます。

とすると、データサイズもコピーの方式も揃えて、[0] * (1024**2 * 128)と b"A" * (1024**2 * 128)の対決ではなく[0, 1] * (1024**2 * 64)とb"ABCDEFGH" * (1024**2 * 128)の対決にすればどちらも同じくらいの速度になる、ということになりそうです。ところがそう簡単にはいかなくて、依然として2.5倍ほどlistのほうが遅くあまり説明になっていません。

他に考えられそうなのは、処理本体ではないので読むのを飛ばしたメモリ確保の方法になります。ところが軽く追った感触では最終的に同じ関数を中で呼んでいそうに見えて結局よくわからずじまい、というところで時間切れです。

Exステージはクリアならず……