Cython(8) 素数配列を作り、素数リストを返す


2022年 09月 15日

Cythonのオンラインマニュアルで、Python内で配列を使う方法が紹介されている。
しかし、Python自体には配列はなく、配列を使いたいとなったらNumPyの配列を使うのが普通と思うが、違う方法のプログラムが載っていた。
配列だけでなく型も宣言していて、その方法がCythonの型宣言機能もPythonの中で利用するというものである。
そのため、最初に import cython がある。

# [Python] cythonの型宣言と配列を利用した素数の計算

import cython

def primes(p_num:cython.int):
    MAXSIZE = 10000
    p_array: cython.int[MAXSIZE]
    p_len: cython.int = 0
    n: cython.int = 2
    p: cython.int

    p_array = [0] * MAXSIZE

    if p_num > MAXSIZE:
        p_num = MAXSIZE
    
    while p_len < p_num:
        for p in p_array[:p_len]:
            if n % p == 0:
                break
        else:
            p_array[p_len] = n
            p_len += 1
        n += 1
    return list(p_array[:p_num])

型宣言部分だけハイライトしておいたが、整数型を宣言するのに単にintではなく、Pythonのファイルで、C言語の整数型宣言になっている訳ではない。

In [1]: from primes_cyarray import primes

In [2]: timeit(primes(10000)) 
3.12 s ± 32.6 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

なんだか遅くなった。
Pythonの中で無理やりCythonの型宣言を行ったためゴチャゴチャしたプログラムになっており、実行時間が長くなるのも当然な気がする。

Pythonの中にCythonの型宣言をするのではなく、Cythonに書き直してみた。

# [Cython] 配列を利用した素数の計算

from cython cimport boundscheck, wraparound

@boundscheck(False)
@wraparound(False)

def primes(int p_num):
    MAXSIZE = 10000
    cdef int[10000] p_array
    cdef int    p_len = 0
    cdef int    n = 2, p
    cdef int    idx
    
    while p_len < p_num:
        for p in p_array[:p_len]:
            if n % p == 0:
                break
        else:
            p_array[p_len] = n
            p_len += 1
        n += 1
    return [p_array[idx] for idx in range(p_len)]

型宣言を変数毎に宣言するのではなく、cdef: を用いて、次行以降にcdef無しの型宣言を並べてみたが、この書き方のほうがスッキリする。

    cdef:
        int p_array[10000]
        int p_len = 0
        int n = 2
        int p, pr
In [1]: from primes_cyarray import primes 

In [2]: timeit(primes(10000)) 
193 ms ± 636 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

C++のvectorを使ってみる

配列のときと同じ様に、Pythonプログラムから直接C++のvectorを使うプログラムは書けるのだが、Cythonのマニュアルに書いているように、そのような呼び出しは、Pythonの標準ではサポートしていない。実際、Cythonのマニュアルの様に書くと、実行時にエラーが発生するので、これ以上の説明は止める。

Cython版では、C++のvectorを使うため、CではなくC++のコードに変換されるよう指定する。

また、vectorを使うために、vectorをcimportしておく。import ではない。

素数のvectorの名前はp_vectorとし、要素の型はintとし、vectorの初期リザーブサイズを指定して、希望数の素数が入るだけのサイズにしておく。

# [Cython] C++のvectorを利用した素数の計算

# distutils: language=c++

from libcpp.vector cimport vector
from cython cimport boundscheck, wraparound

@boundscheck(False)
@wraparound(False)

def primes(int p_num):
    cdef vector[int]    p_vector
    p_vector.reserve(p_num)
    cdef int    p, n = 2

    while p_vector.size() < p_num:
        for p in p_vector:
            if n % p == 0:
                break
        else:
            p_vector.push_back(n)
        n += 1
        
    return p_vector
In [1]: from primes_vector import primes 

In [2]: primes(10)  
Out[2]: [2, 3, 5, 7, 11, 13, 17, 19, 23, 29]

In [3]: primes(10000)[-10:] 
Out[3]:
[104677,
 104681,
 104683,
 104693,
 104701,
 104707,
 104711,
 104717,
 104723,
 104729]

In [4]: timeit(primes(10000)) 
193 ms ± 2.12 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

前回(Cython(7))で185 msだったのが、193msになったので、多少遅くなったが、僅かである。
どうも、リストでもvectorでも、速度にあまり差がでないようである。
実は、Pythonのリストは、いわゆるリストではなく、動的配列(dynamic array)として作られる。だから、当然vectorとの差がほとんど出ない訳である。
ということは、リストなら速いはずのことが、とても遅くなってしまうことがあるということだが、そのあたりの話は別の機会にしよう。

関数の最後のreturn p_vectorにより素数のvectorを戻り値にしているが、関数はdef primesにより宣言されているのでPythonから呼べるようになっている。そのため、vectorの戻り値は、Pythonのリストに自動的に変換され、Python側からすると、素数のリストを返すPythonの関数に見える。

なお、p_vector.reserve(p_num)により、指定個数の素数を入れられるvectorサイズを予約したのだが、この行をコメントアウトしても、実行時間に変化は見られなかった。vectorのサイズの管理は、プログラマが気にしなくても適当に上手に管理してくれているようだ。