日本ソフトウェア科学会第 35 回大会 (2018 年度) 講演論文集

対称な関数呼び出しの融合佐藤 重幸

科学技術計算や機械学習において，対称な関数呼び出しが用いられることが良くある．例えば，物体 x から y への物理作用を計算する f px, yq という関数を考える． 作用・反作用の原理から，作用 f px, yq を計算するときには，反作用 f py, xq も計算して利用することが，物理モデルから自然である． この f px, yq と f py, xq は殆ど同じ計算であるため，両方を呼び出す場合には，両方をまとめて計算する効率的な融合関数によって代用できる． 本研究では，作用を定義する関数から，SMT ソルバとデータフロー解析に基づいて，相互作用を効率的に計算する融合関数を導出する方法を提案する．

コストだけで済むことになる．

1 序論

実際のプログラムでは，このような作用の計算式

時間発展系シミュレーション [1, 2] や様々な機械学
習アルゴリズム [6, 9] は，N 体問題，すなわち N 点

は，カーネルと呼ばれる関数（手続き）として実装される．このとき，f⃗12 を計算して p2 に積算するカー間の相互作用を計算する問題として定式化される．こネル gravity は，その重力作用の対称性に由来して，の N 点間の相互作用は，全ての 2 点間の相互作用の次のように対称な関数呼び出しが生じることになる．

和として定義される．そしてこの相互作用は，対称性があることがよくある．具体例として，物体 p1 から
p2 に作用する重力 f⃗12 を考える．

f⃗12 “

Gm1 m2¨ p⃗x1 ´ ⃗x2 q||⃗x1 ´ ⃗x2 ||3

ここで，G は重力定数であり，m1 と m2 はそれぞれ

p1 と p2 の質量，⃗x1 と ⃗x2 はそれぞれ p1 と p2 の位置ベクトルを意味する．作用・反作用の原理から，p2
から p1 への重力 f⃗21 は，´f⃗12 である．これは，上式における記号の添字 1, 2 を反転させることで代数的に確認できる．したがって，f⃗12 と f⃗21 をまとめて計算するならば，f⃗12 の計算コストと，ベクトル反転の˚ Fusion of Symmetric Function Calls.

This is an unrefereed paper. The copyright belongs
to the author. This is ongoing work. If you would
like to cite this paper, please contact the author forthe latest version.
Shigeyuki Sato, 高 知 工 科 大 学, Kochi University ofTechnology.

gravity ( p1 , p2 )gravity ( p2 , p1 )# p2 に f⃗12 を積算# p1 に f⃗21 を積算

この様な対称な関数呼び出しの対があるとき， p1と p2 の作用を，上述のように効率良くまとめて計算するカーネル（これを相互作用カーネルと呼ぶ）を用いることが望ましい．典型的な N 体ソルバ（例えば ExaFMM†1 ）では，この相互作用カーネルが元のカーネルと共にハードコードされて利用される．しかし，元のカーネルと，定義上重複部分の多い相互作用
カーネルをハードコードすることは，記述性や保守性の観点で望ましくない．
では，現代の最適化コンパイラ（例えば GCC やClang/LLVM）が対称な関数呼び出しの対があった
ときに，コード最適化によって最終的に相互作用カー
ネルと同等のコードを吐くかというと，そんなことはない．2 つの gravity の関数呼び出しをインライン†1 https://github.com/exafmm/exafmm

展開してまとめて最適化を適用しても，伝統的な低水準の演算に関する冗長性除去では，前述のように融合した効率良い関数は出てこない．

この例からわかるように，相互作用カーネルを構成するには，元のカーネル中に現れる計算式を対称（Sym） カーネル引数を反転させたときにそこで本研究では，対称性のある作用を与えるカーネルから，効率的な相互作用カーネルを導出する高水
準最適化技法を提案する．提案手法は，古典的な冗長性除去手法とは異なり，SMT ソルバを用いて計算式の対称性と歪対称性を判定し，必要に応じて反転した値が等しい．
歪対称（Skew） カーネル引数を反転させたときに，符号が反転する．
複製（Dup） 上 2 つ以外．カーネル引数を反転させた式を複製する必要がある．

計算式を導入することで，相互作用カーネルを構成す

これら 3 つに分類することが役立つ．計算式が分類る．この手法を，SSA 形式上のデータフロー解析とされていれば，相互作用カーネルのコードは，元のして定式化した．本論文では，前述の重力カーネルにカーネルのコードに沿って素直に生成できる．した

対する予備実験の結果も示す．

がって，本論文で扱う最適化の技術的要点は，カーネ
ル中の計算式を分類する静的解析である．

2 動機付けする例
図 1 は，前節で挙げた重力作用カーネルの Python

3 入力言語と制約に よ る 具 体 例 で あ る ．尚 ，本 論 文 で は ，プ ロ グ

本研究において最適化の対象となるカーネルは，任ラ ム 記 述 に は Python（ の サ ブ セット ）を 用 い る ．

意の関数ではない．ある程度お行儀の良いプログラム

gravity(p1,p2) は，p1 から p2 への重力を計算し，であることを仮定する．本節では，想定する入力プロ

それにしたがって p2 の加速度を積算する関数である．
この gravity における局所変数を観察すると，次グラムの制約について説明する．
まず，入力言語の構文は，図 3†2 に示される Pythonのことが分かる．まず，diff は，⃗
x1 ´ ⃗x2 に対応し，のサブセットを扱う．意味論は，動的言語の特徴を備

p2 から p1 への反作用を計算する時には符号反転した値を用いることができる．r2 は，||⃗
x1 ´ ⃗x2 || に対応

えた Python ではなく，それを型付き言語の形で制約した RPython†3 に従う．したがって，し，p1 から p2 への作用を計算するときも，その反作‚ 型宣言は無いものの，型推論によって全ての変2用を計算する時も，その値は等しい．r3 も同様である．delta は，p1 から p2 への作用と p2 から p1 の作用を計算するときに，異なる値となる変数である．以上より，重力相互作用カーネル gravity_mutual が
図 2 のように実装できる．ここで，delta_r は，引数が反転した場合の delta の計算式の値を保持する．
この delta と delta_r の計算には共通部分がある
が，乗算*の結合性と可換性を利用しつつ伝統的なコンパイラ最適化を適用することで，次のように冗長性を除去できる．
tmp = G / r3 * diff [ i ]p2 . acc [ i ] = p2 . acc [ i ] + tmp * p1 . mass
p1 . acc [ i ] = p1 . acc [ i ] - tmp * p2 . mass

この 3 行目だけが，反作用の計算に追加で必要な部分である．これは，元の gravity 関数を適用する場合に比べて，大幅に効率化されている．

数に型が付く．

‚ 変数への再代入は許すが，初期化値と同じ型の値に制限される．

‚ 未初期化変数を参照しうると制御フローから判断されるプログラムは，静的に弾く．

‚ カーネル引数の型（それが持つ全ての属性の型）は具体的に与えられている．
という 4 点が仮定される．
さらに，解析を簡単にするために，別名に関する制約を加える．

‚ カーネルの引数 p1 と p2 ，及びその属性は別名†2 本論では図 3 で用いられた記号は，集合としての解釈を持つ．つまり，Exp は式を意味する抽象構文木の集合を意味する．同様に，Var は，変数記号の集合ではなく，変数出現の集合である．
†3 https://rpython.readthedocs.io/en/latest/def gravity ( p1 , p2 ):diff = [ p1 . pos [ i ] - p2 . pos [ i ] for i in range ( D )]r2 = sum ( diff [ i ]* diff [ i ] for i in range ( D ))r3 = r2 * sqrt ( r2 )for i in range ( D ):delta = G * p1 . mass / r3 * diff [ i ]p2 . acc [ i ] = p2 . acc [ i ] + delta
図1

重力作用カーネル．引数 p1 から p2 への作用を計算して，それを p2 に累積させる関数．

def grav ity_mut ual ( p1 , p2 ):diff = [ p1 . pos [ i ] - p2 . pos [ i ] for i in range ( D )]r2 = sum ( diff [ i ]* diff [ i ] for i in range ( D ))r3 = r2 * sqrt ( r2 )for i in range ( D ):delta
= G * p1 . mass / r3 * diff [ i ]delta_r = G * p2 . mass / r3 * - diff [ i ]p2 . acc [ i ] = p2 . acc [ i ] + delta
p1 . acc [ i ] = p1 . acc [ i ] + delta_r
図2

重力相互作用カーネル．引数 p1 から p2 への作用とその反作用を計算して，それぞれを p2 と p1 に累積させる関数．

FunctionDef

::“

Block

::“

Stmt`

Stmt

::“

for Var in range (Exp) : Block

AssignStmt|

AssignStmt::“

Target = Exp|

Target[Exp] = Exp|

Var = [Exp for Var in range (Exp) ]|Var = sum (Exp for Var in range (Exp) )|Var = ϕpVar , Var q

Target::“Var | pp1 |p2 qp.Attr q`

Exp

::“

Target | Const | Exp.Attr | Exp[Exp]|図3

def Id (p1 ,p2 ) : Block

Exp BinOp Exp | UnOp Exp | Exp if Exp else Exp

入力言語の構文．Python サブセット．タイプライタ体は終端記号を意味し，` は 1 回以上の繰り返しを意味する．

ϕ 関数の代入文は実際の入力には含まれず，SSA 変換で中間的に導入され，SSA 逆変換で最終的に除去されるものである．

を持たない．

‚ 配列変数は別名を持たない．

が成立しないからである．
実装上の都合から，式 Exp については，既存のここで，別名を持たないというのは，入力においてSMT ソルバがそのまま解釈できるような範囲の表現

別名があったとしても，別名解析によって全て除去できるということを含意している．特に前者の制約

能力に制限されている．したがって，再帰関数は呼び出せないし，sum は式 Exp には含まれない．Z3†4 なは，本論文で扱う最適化において事実上必須である．

どの SMT ソルバと同様に，マクロとして解釈できる

gravity_mutual の例で示したように，p1 と p2 を反関数は許すが，簡単のために，その呼出しは前処理転させた計算式を効率的に実装できなければ，最適化†4 https://github.com/Z3Prover/z3

によって全て展開されているものとする．実用上のある 0 に対応するものであり，次小節での一般化のた配慮から，組み込み関数（sqrt などの数学関数）もめに導入される．

許可するが，SMT 式上では uninterpreted functions
として解釈される．

変数解釈 E :: Var Ñ L の下で反転した式への変換

RE は，図 4 のように構文に沿って素直に定義でき最後に，相互作用カーネルを構成するための前提とる．そして，RE を用いれば，変数解釈 E の下で式して，作用と反作用を独立に計算して積算できるといExp を L 上で分類する関数 check E は，図 5 のようう性質が必要である．これを形式的に言い換えると，に即座に得られる．この check E に基づいて，代入文

p1 の属性 a が更新されるとき，p2 の属性 a の値に依の右辺の評価と E の更新を繰り返す抽象解釈を行え存してはならなず，p1 と p2 を逆にしても同様である

ば，図 1 の例は，2 節で述べたように最適化できる．

ということである．これは言語レベルではなくプログ

RE において，Dup である変数について，対応す

ラムレベルの制限であり，本論文では事前解析によっる fresh な変数を用いて SMT ソルバで検証することてこれが保証されていると仮定する．その事前解析

は，安全な近似である．例えば加算が結合的且つ可は，field-sensitive taint analysis に相当するものに換的な数値 x と y について，x ` y は対称であるのなる．前述の別名に関する制約から，その実装は難し

で，s = x + y の代入文を解釈すると，s は Sym といものではない．

なる．しかし，t = x と s = t + y に分解し，それ以上の入力に関する制限は，プログラム一般を記述

ぞれの代入文を順に評価すると，t も s も Dup となするには制限が厳しいものだが，局所的な計算を記述

る．Dup ではない変数を Dup と判定しても，その変するカーネルにおいては，妥当な範囲である．実際，数が利用される部分に対応するコードの質が落ちる

図 1 の重力作用カーネルで示されるように，N 体問だけで，誤った結果は生まないので安全である．

題のカーネル程度ならば，素直に記述できる．

尚，この観察から，RE において，v 1 を導入する代
わりに，その定義式を可能な限り展開して SMT ソル

4 提案手法

バで検査すれば，check E の判定精度が向上しうるこ

2 節で述べた重力作用カーネルを用いて，提案するとがわかる．しかし，一般に，事態はそれほど簡単で最適化手法を述べる．

はないことを 5 節にて議論する．

4. 1 基本的なアイデア

4. 2 SSA 形式上の大域的解析

提案手法の基本的なアイデアは，2 節の説明の通り，

さて，check E があれば図 1 の例が扱えると単純に1) 計算式を Sym ，Skew ，Dup に分類し， 2) その値述べたが，これは gravity の例に，たまたま変数のを保持する変数にラベル付けし，3) 各代入文につい再定義がなかったからである．gravity と同等の計て作用・反作用の 2 通りの意味解釈に沿ってコード生

算を，累積的計算として実装した gravity2 （図 6）成するということにまとめられる．

については，check E を考えるだけでは，きちんと扱

今，入力の Exp は SMT ソルバで解釈可能なものに限定しているので，任意の e :: Exp について，そのえない．制御フローを捉えた大域的解析として定式化すべきである．

中のカーネル引数を反転させた式 er :: Exp を構成で本研究では，SSA 形式上のデータフロー解析としきれば，e “ er も，e “ ´er も検証できる．e :: Expて，この大域的解析を定式化する．まず，次の有界半の解釈を定めるには，それを評価する環境，つまり変束 pL, [q を考える．

数の解釈が必要である．変数の値は，対称性・歪対称性を判定できる程度の情報であれば充分なので，

L “ tZero, Sym, Skew , Dupuと抽象化する．ここで，Zero は，対称かつ歪対称で‚ 任意の x P L について x [ Dup “ Dup （Dupが最小元であり [ の吸収元）‚ 任意の x P L について x [ Zero “ x（Zero が最大元であり [ の単位元）

RE :: Exp Ñ Exp
RE pp1 q “ p2RE pp2 q “ p1$’’-v pEpvq “ Skew q’’&
RE pvq “ v 1 pEpvq “ Dupq
’
’
’
’
%v
potherwiseq$’’-v[e] pEpvq “ Skew q’’&
RE pv[e]q “ v 1 [e] pEpvq “ Dupq
’
’
’
’
%v[e]
potherwiseq
RE pcq “ cRE pe.aq “ RE pe1 q.a
RE pe1 [e2 ]q “ RE pe1 q[RE pe2 q]
RE pe1 f e2 q “ RE pe1 q f RE pe2 q
RE p a eq “ a RE peq
RE pe1 if e2 else e3 q “ RE pe1 q if RE pe2 q else RE pe3 q
図4

変数解釈 E :: Var Ñ L の元で反転した式への変換 RE ．ここで，v ，e，f，a は，それぞれ Var ，Exp ，BinOp ，UnOp 上を動くメタ変数であり，v 1 は，反転した場合の v に対応する fresh な変数である．

‚ Sym [ Skew “ Dupcheck E :: Exp Ñ L

これは，L による抽象化の意味から自然な半束であcheck E p0q “ Zero
$
’
’
Sym
’
’
&
check E peq “ Skew
’
’
’
’
%Dup
図5

る．これをデータフロー値の定義域とする．

pRE peq “ eqpRE peq ‰ e ^ RE peq “ ´eqので，代入文の解釈 F :: AssignStmt Ñ L があるとpotherwiseq

都合がよい．変数解釈 E :: Var Ñ L の下で代入文を変数解釈 E :: Var Ñ L の元で Exp の対称性・歪対称性を検査する check E ．ここで，0 は，数値 0 と等価な定数を意味するメタ変数．

def gravity2 ( p1 , p2 ):diff = [0.0 for i in range ( D )]r2 = 0.0
for i in range ( D ):diff [ i ] = p1 . pos [ i ] - p2 . pos [ i ]r2 = r2 + diff [ i ]* diff [ i ]r3 = r2 * sqrt ( r2 )for i in range ( D ):delta = G * p . mass / r3 * diff [ i ]q . acc [ i ] = q . acc [ i ] + delta
図6

コード生成時には各代入文について扱いを定めたい

累積的に記述された重力作用カーネル．図 1 のプログラムと同等．

解釈する関数 eval E を図 7 に与える．そして，eval E
によって定められるデータフロー値に基づいた大域的解析は，ワークリストアルゴリズムとして，図 8 に示される F で実装できる．

F において注意すべきは，Use における配列変数の扱いである．v を配列変数（の出現）としたとき，

Usepvq は，v の要素を更新する全ての AssignStmtからなる集合を返すとする．これは，全ての配列の値が，制御フロー非依存で合流するという安全な近似である．この近似によって，配列変数については，SSA
形式の対象から除外する．
入力言語の制約から，配列変数についても SSA 変
換し，制御フローに沿って値を合流させる解析も可能であるが，あまり賢明な選択ではない．まず，現実的eval E :: AssignStmt Ñ Leval E px = eq “ check E peq
eval E px = [e for _ in range(_)]q “ check E peq
eval E px = sum(e for _ in range(_))q “ check E peq
eval E px1 = ϕpx2 , x3 qq “ Epx2 q [ Epx3 q$’&check E pe3 q pe1 P Var ^ check E pe2 q Ě Symq
eval E pe1 [e2 ] = e3 q “’%Duppotherwiseq
eval E p_._ = _q “ Dup.
図7

変数解釈 E :: Var Ñ L の元での AssignStmt のデータフロー値を与える eval E ．ここで，_ は don’t care 項を意味する．

F :: PpAssignStmtq Ñ pAssignStmt Ñ Lq
FpSq “ F 1 pS, ts ÞÑ Zero | s P Suq
F 1 pH, F q “ F
F 1 pW Y tsu, F q “ F 1 pW Y W∆ , F 1 q,$’
&H
(F 1 psq “ F psq)
where W∆ “’
%Usepsq (otherwise)
F 1 “ pF zts ÞÑ F psquq Y ts ÞÑ F psq [ eval F ˝Def psqu.
図8

SSA 形式上のデータフロー解析アルゴリズム F ．入力は AssignStmt の全体集合．ここで，Usepsq は，

s P AssignStmt で定義される変数を利用する AssignStmt の集合，Def pxq は，x P Var を定義する AssignStmt を意味する．SSA 形式を仮定しているために，Def pxq は一意に定まる．

なカーネルコードの記述において，精度が向上する例

式に制限してしまっても問題ないとも考えている．

があるとは想定しがたい．そして，例え精度が向上したとしても，配列変数のコピー代入はコストが大きい
ために，性能利得を最大化するコード生成（SSA 逆変換）が非自明になる．したがって，配列変数の値に

4. 3 SSA 形式の利点
ところで，提案する最適化は，SSA 形式との相性がとても良い．

ついてはフロー非依存に解析するのが妥当である．

SSA 形式は一般に，データフロー解析の実装を簡最後に，F 及び evalE の定式化に関して，range

単にするとか効率的にするとか言われる．確かに Fの引数を無視している点には留意が必要である．これにおいて，変数解釈 E を，F ˝ Def という（実装レは，range の引数式が必ず Sym になると仮定していベルでは大域的な）写像の合成として表現できることることを意味する．F の定式化の簡潔さためにそれは，定式化の見通しが良い．しかし，E をデータフを仮定したが，range の引数を考慮に入れる実装にす

ロー値として定式化すれば，同等の解析結果を通常形

る拡張することは難しくない．Use で for 文も返すよ

式の上で得ることは特段難しくない．

うにした上で，range の引数が Sym でないときには，本最適化において，SSA 形式が明白に有利に働く

その for 文全体を Dup だと判定すればよい．それ故のは，静的解析ではなく，コード生成においてであに，3 節で range の引数について制限を特別加えなる．通常形式である場合，変数解釈がフローに依存かった．しかし実際上は，自明に Sym であるようなするため，変数 x の値が，ある点では Sym で，別のところでは Dup になることがある．このとき，初め
1バで検証した方が良いのだろうか．とりわけ歪対称性て Dup になった点で，その定義式の複製と x の導入については，そんなことはない．単純な例としては，がなされるべきであり，さらに x1 の値は後に x と合x ´ y は歪対称であるが，d ´ px ´ yq は対称でも歪対流させなければならない．これはまさに，x について称でもない．x ´ y の部分式に関してのみ歪対称性を

SSA 変換と逆変換を求めていることと等しい．別の利用した変換を適用しないと，Dup と判断される計

言い方をすると，各変数の定義式の複製点は，単一代算式が増大し，生成されるコードの品質が低下する．

入性から自明に定まり，その複製した値を伝播させる

この様に，特に歪対称性に基づく変換は機械的にため変数導入も，ϕ 関数の代入文によって自然に達成

扱いにくく，既存の最適化コンパイラが利用しないのされる．したがって，本最適化において，SSA 形式

は，設計選択として妥当かつ合理的である．ではなを経由することは，定式化として見通しが良く，実装

ぜ，それほどの扱いにくさを内包する最適化が，前節上合理的である．

で述べた単純な解析アルゴリズムで実装でき，人手で変換することも容易なのだろうか．

5 最適化の難しさ

そのキーは，プログラマの正気度である．何らかの前節で述べた静的解析アルゴリズムは，SMT ソル

数理モデル式をプログラムとして記述するとき，正気バの恩恵もあって，とても単純なものである．実際，なプログラマなら，意味のある区切りで変数に式値

gravity の定義から，gravity_mutual の定義を人手を代入する．したがって，式の性質を判定する粒度とで作るのは，ほとんど機械的な作業である，では，対して代入文区切りを採用する F は，このプログラマ称性・歪対称性を活用する最適化は簡単なものなのでの正気さ，言い換えれば，イディオム性に依存していあろうか．

ることになる．また，人手で容易に変換できるのも，

少なくとも一般には，提案する最適化は簡単では暗黙にそのイディオム性にしたがって代入文を解釈しない．これは 6 節でも述べるように，最新の最適化ているからである．つまり，計算の解釈にちょうど良コンパイラ（GCC や Clang/LLVM）でも，gravity

い区切りが残されているという，プログラムのお行儀のような単純な例ですら対象の冗長性が除去されなの良さが，最適化の実装を単純化させている．

ことがその証左である．

まとめると，提案する最適化は，網羅的且つ完全な既存のコンパイラ最適化が効かない主な原因は，歪実装が難しい．しかし，入力プログラムにおけるイ対称性の扱いにくさにある．歪対称性を利用する時

ディオム性を頼るならば，単純に実装できる．イディには，x ´ y “ ´py ´ xq のように符号反転を外に括オム性に基づくヒューリスティクスが重要であり，そり出すことになるが，これ自体は演算量を増やす変れなしに実装するのは合理的とは言い難い．最適化器

換である．その露出した y ´ x が共通部分式にならにとっての頼みの綱は，プログラマの正気度である．

ない限り，その変換は適用すべきではない．その上，仮に y ´ x が共通部分式として除去されて ´t に置換

6 予備実験

されたとしても，それ単独では依然として演算量が減提案手法をプロトタイプ実装し，生成されたコード

らない．gravity の例のように，その ´t を利用するの予備評価を行った．

計算式が対称・歪対称になって初めて目立った性能利3 節で述べた入力言語は，Python3.6 の言語内 DSL得が生じる．つまり，歪対称を利用した変換自体が，として実装し，最適化器と C++14 へのコード生成

巨視的な解析を要する高水準最適化なのである．計算式を細分化して微視的に冗長性を発見・除去する伝統

器を Python ライブラリとして実装した．SMT ソル
バは，Z3†5 を採用し，その Python バインディング

的な冗長性除去の枠組みには，全く適してない．

である Z3Py を使って最適化器に組み込んだ．

巨視的な解析が必要だというなら，4 節で述べたよ
うに，可能な限り式をインライン展開して SMT ソル†5 https://github.com/Z3Prover/z3

2 節で挙げた重力作用カーネルに提案する最適化を
4

るように一般化されている．Ginsbach ら [5] は，イ

適用し，アセンブリサイズの変化と N “ 10 のとき

ディオム記述言語を設計し，それで記述されたイディの N 体問題の直接解法の実行時間の変化を調べた．

オムを SMT 制約式として翻訳して，LLVM IR 上で生成された C++14 コードのコンパイルには，Clang

イディオムにマッチさせる枠組みを開発した．これに6.0.4 と G++ 7.3.0 を用いた．最適化レベルどちら

よって，レガシーコードの中にある行列積コードを，も最高（Clang は-Ofast，G++が-O3）とし，両方BLAS ライブラリの呼出しに置換するという最適化に-ffast-math を指定した．これは，浮動小数点数を実装した．

の加算と乗算について，結合性と可換性を仮定した最本研究では，対称なカーネル呼出しを融合させた，適化を用いるフラグである．アセンブリサイズを比較

相互作用カーネルを構成した．確かに関数を融合させする時には，直接解法に関わらないコードは全て除去てはいるが，これはプログラム変換の文脈では，関数しつつ extern 大域変数を利用して，最適化によって融合というよりも，tupling [3] に相当する変換であ

直接解法の計算部分が除去されないようにした．

る．この tupling 変換は，プログラム一般に適用する

TM
R Core
実験に用いた計算機は，Intel⃝
i7 CPU

戦略を作るのが困難であることが知られている．Hu

M620 (2.67 GHz) と 8 GB DDR3-1333 を備えたラッら [7] は，再帰データ構造に沿った行儀の良い再帰関

プトップ PC であり，Debian x86_64 (Linux kernel

数に関して，系統的に tupling を適用する手法を開発

4.17.0) を実行していた．実行時間は，1 回の warm-up

した．適用戦略を立てにくい変換を，よくあるイディ実行の後に，15 回測定し，その中央値を採用した．

オムを捉えて効果的に適用するという点は，本研究と表 1 に実験結果をまとめる．最適化前は，Clang

共通している．

でも G++でも実質的にほぼ同じアセンブリが生成さ
れ，実行時間もほぼ一致していた．最適化後は，どち8 結論

らもコードサイズが 14%程度削減され，実行時間は

本研究では，対称な引数で呼び出される関数についClang のとき 43%，G++は 31%削減された．この結て，対称な呼出し対を効率的に 1 つに融合させる最適果より，提案する最適化技法は，最先端の最適化コン化手法を提案した．提案する最適化は，既存の最適化

パイラでもカバーできていない冗長性除去を実現し，コンパイラでは捕捉できていなかった内在する冗長性

実行性能の大幅な向上をもたらすことが確認できた．

を，計算式の歪対称性に基づいて特定し，除去する．
現状の最適化器はプロトタイプだが，簡単な例に

7 関連研究ついては動作している．複素数・数学関数・ラムダ式

我々の知る限り，本論文で提案する最適化それ自体（マクロ）についての完全に実装すれば，ExaFMM

は，コンパイラ最適化の文脈で扱われた事はない．しの複雑なカーネルにも，そのまま対応できるようになかし，その根底ある，等価性に基づく式変換 [8]，代る．特に 3 角関数の歪対称性については，単なる符数的性質に基づく高水準最適化 [4]，SMT ソルバを

号反転を導入するのではなく，π の周期性を使う必要利用したイディオム認識 [5] というアイデアは，コンがあり，SMT 式への翻訳には一工夫を要する．そのパイラ最適化の文脈では珍しいものではない．

翻訳がハードコーディング可能であることは確認でき

Tate らが提案した equality saturation [8] は，グているが，拡張可能な規則の形で，系統的に実装でき

ラフに基づく中間表現上の最適化手法であり，ループ

ることが望ましい．そうすれば標準ライブラリ以外の越しの代数変換に基づく最適化を主眼として設計さ関数についても，容易に対応できるようになる．

れた．彼らが扱った最適化自体は，古典的で微視的な

5 節で述べたように，検証するのにちょうど良い式

冗長性除去である．Ding らは [4]，三角不等式に基づのまとまりを作り出すことは，ヒューリスティクスとく演算子強度軽減を実装した．通常の辺に関する三して重要な要素である．これもハードコーティング

角不等式に加えて，角度に関する三角不等式を用いするのではなく，拡張可能な規則の形で設計できるこ

表1

重力作用カーネルを提案手法で最適化したことによるコードサイズとコード品質の変化．

gravityClangG++gravity_mutual
Clang

G++

アセンブリサイズ [LoC]

165

168

142

145

実行時間 [ms]

1738

1738

997

1202

とが望ましい．将来的には，高水準のイディオム識別器 [5] の形で基盤を設計し，他の様々なイディオム・

API 固有の最適化と統合して扱えるようにすることが合理的だと考えている．
謝辞
本研究は，JST ACT-I の支援を受けたものである．
参 考 文 献[1] Appel, A. W.: An Eﬃcient Program for ManyBody Simulation, SIAM J. ScI. Stat. Comput.,Vol. 6, No. 1(1985), pp. 85–103.
[2] Barnes, J. and Hut, P.: A hierarchical OpN log N q
force-calculation algorithm, Nature, Vol. 324(1986),
pp. 446–449.
[3] Chin, W.-N.: Towards an Automated Tupling
Strategy, Proc. the 1993 ACM SIGPLAN Symposium on Partial Evaluation and Semantics-based
Program Manipulation (PEPM ’93), ACM, 1993,
pp. 119–132.
[4] Ding, Y., Ning, L., Guan, H., and Shen, X.: Generalizations of the Theory and Deployment of Triangular Inequality for Compiler-Based Strength Re-

duction, Proc. the 38th ACM SIGPLAN Conference
on Programming Language Design and Implementation (PLDI ’17), ACM, 2017, pp. 33–48.
[5] Ginsbach, P., Remmelg, T., Steuwer, M., Bodin,
B., Dubach, C., and O’Boyle, M. F. P.: Automatic
Matching of Legacy Code to Heterogeneous APIs:
An Idiomatic Approach., Proc. the 23rd International Conference on Architectural Support for Programming Languages and Operating Systems (ASPLOS ’18), ACM, 2018, pp. 139–153.
[6] Gray, A. G. and Moore, A. W.: ‘N-Body’ Problems in Statistical Learning, Proc. NIPS ’00, 2000,
pp. 521–527.
[7] Hu, Z., Iwasaki, H., Takeichi, M., and Takano,
A.: Tupling Calculation Eliminates Multiple Data
Traversals, Proc. the 2nd ACM SIGPLAN International Conference on Functional Programming(ICFP ’97), ACM, 1997, pp. 164–175.
[8] Tate, R., Stepp, M., Tatlock, Z., and Lerner,
S.: Equality Saturation: A New Approach to Optimization, Proc. the 36th Annual ACM SIGPLANSIGACT Symposium on Principles of Programming
Languages (POPL ’09), ACM, 2009, pp. 264–276.
[9] van der Maaten, L. and Hinton, G.: Visualizing Data using t-SNE, J. Mach. Learn. Res.,
Vol. 9(2008), pp. 2579–2605.