はじめに
メタプログラミングとは
メタプログラミング (metaprogramming) とはプログラミング技法の一種で、ロジックを直接コーディングするのではなく、あるパターンをもったロジックを生成する高位ロジックによってプログラミングを行う方法
出展 : wikipedia
要するに、コードを生成するコードを書いて複雑なコードを書いちゃいましょう！ってこと（だと私は解釈しています。間違っていたら教えてください）。
概要
Pythonの記号計算用ライブラリであるSympyを使って数式を生成し、C言語にコンバートします。SympyはMathematicaのように数式処理に使われることが多いようですが、数式を様々な言語にコンバートする機能も兼ね備えています。私はあまり使ったことがありませんが、ヘッダファイル付きで関数を出力する機能もあります。
この記事では入門的に、n×n行列の掛け算のコードをSympyでメタプログラミングしていきます。
Sympyの記法や使い方は以下にまとめられているので、興味のある方は読んでみてください。
インストール
Sympyのインストールはpipでできます。
pip install sympy
ダウンロードの場合は以下から。
インポート
まずはSympyをインポートします。
from sympy import *
今回は行列計算を実装したいので、Matrixもインポートしておきます。
from sympy import Matrix
これでMatrixも宣言できるようになりました。
シンボルの作成
数式に用いる変数はsymbolとして事前に宣言しておく必要があります。例えば、xとyをsymbolとして宣言する場合は
x = symbols('x')
y = symbols('y')
とすればOKです。今回は行列計算のために、n×nのMatrixの要素としてsymbolを定義していきます。
n = 3
A = Matrix([[symbols('a['+str(j)+']['+str(i)+']') for i in range(n)] for j in range(n)])
B = Matrix([[symbols('b['+str(j)+']['+str(i)+']') for i in range(n)] for j in range(n)])
これで行列A,Bがそれぞれ以下のように宣言されました。ここではn=3としました。
⎡a[0][0] a[0][1] a[0][2]⎤
⎢ ⎥
⎢a[1][0] a[1][1] a[1][2]⎥
⎢ ⎥
⎣a[2][0] a[2][1] a[2][2]⎦
⎡b[0][0] b[0][1] b[0][2]⎤
⎢ ⎥
⎢b[1][0] b[1][1] b[1][2]⎥
⎢ ⎥
⎣b[2][0] b[2][1] b[2][2]⎦
行列の出力はpprint(A)と記述するとこのように視覚的にわかりやすい形になります。
行列の掛け算
先ほどの行列A,Bの積を行列Cとします。行列の掛け算は簡単に、以下でOKです。
C = A*B
今度はprint(C)で中身を見てみると
Matrix([[a[0][0]*b[0][0] + a[0][1]*b[1][0] + a[0][2]*b[2][0], a[0][0]*b[0][1] + a[0][1]*b[1][1] + a[0][2]*b[2][1], a[0][0]*b[0][2] + a[0][1]*b[1][2] + a[0][2]*b[2][2]], [a[1][0]*b[0][0] + a[1][1]*b[1][0] + a[1][2]*b[2][0], a[1][0]*b[0][1] + a[1][1]*b[1][1] + a[1][2]*b[2][1], a[1][0]*b[0][2] + a[1][1]*b[1][2] + a[1][2]*b[2][2]], [a[2][0]*b[0][0] + a[2][1]*b[1][0] + a[2][2]*b[2][0], a[2][0]*b[0][1] + a[2][1]*b[1][1] + a[2][2]*b[2][1], a[2][0]*b[0][2] + a[2][1]*b[1][2] + a[2][2]*b[2][2]]])
と、きちんと計算出来ていることがわかります。
C言語へのコンバート
C言語へのコンバートは
ccode(<数式>, <代入する変数>, standard='C99')
と記述します。今回は行列Cの全要素をコンバート・出力したいので、
for i in range(n):
for j in range(n):
idx = i*n+j
code = ccode(C[idx],assign_to=('c['+str(i)+']['+str(j)+']'), standard='C89')
print(code)
としました。idx = i*n+jは、i, jを一次元に変換する式です。上記コードを実行すると以下が出力されます。
c[0][0] = a[0][0]*b[0][0] + a[0][1]*b[1][0] + a[0][2]*b[2][0];
c[0][1] = a[0][0]*b[0][1] + a[0][1]*b[1][1] + a[0][2]*b[2][1];
c[0][2] = a[0][0]*b[0][2] + a[0][1]*b[1][2] + a[0][2]*b[2][2];
c[1][0] = a[1][0]*b[0][0] + a[1][1]*b[1][0] + a[1][2]*b[2][0];
c[1][1] = a[1][0]*b[0][1] + a[1][1]*b[1][1] + a[1][2]*b[2][1];
c[1][2] = a[1][0]*b[0][2] + a[1][1]*b[1][2] + a[1][2]*b[2][2];
c[2][0] = a[2][0]*b[0][0] + a[2][1]*b[1][0] + a[2][2]*b[2][0];
c[2][1] = a[2][0]*b[0][1] + a[2][1]*b[1][1] + a[2][2]*b[2][1];
c[2][2] = a[2][0]*b[0][2] + a[2][1]*b[1][2] + a[2][2]*b[2][2];
行列Cの各要素が計算され、C言語の形で出力されました。ちゃんと ; も付いていますね。あとはこれを自分のコードに張り付ければOKです。
ソースコード
今回作成した行列計算コードです。
from sympy import *
from sympy import Matrix
from sympy.utilities.codegen import codegen
n = 3
A = Matrix([[symbols('a['+str(j)+']['+str(i)+']') for i in range(n)] for j in range(n)])
B = Matrix([[symbols('b['+str(j)+']['+str(i)+']') for i in range(n)] for j in range(n)])
pprint(A)
pprint(B)
C = A*B
print(C)
for i in range(n):
for j in range(n):
idx = i*n+j
code = ccode(C[idx],assign_to=('c['+str(i)+']['+str(j)+']'), standard='C89')
print(code)
まとめ
Sympyを使って行列の掛け算を計算するコードをメタプログラミングしました。このくらいの計算ならメタプログラミングするまでもなく書けますが、「元の式をtで何階偏微分して～」みたいな複雑な数式を実装する時には非常に有用です。（行列掛け算もfor文がなくなった分少し実行速度が上がるメリットはある？）
ご参考になれば幸いです!