ツリー法、FMMのどちらにしても、直接計算よりも速く計算できるというとこ ろに値打ちがある。アルゴリズムの理論的な研究としては、計算量のオーダー 評価とか理論的なオペレーション数を数えるとかして少なければよい方法であ るということになるが、実際に使う側から見れば、理論的に良いかどうかより も
ツリー法については、多重極展開が 
個あってそれを評価するというの
はもうどうしようもないので、あまり改良の余地はない。これに対し、 FMM
では特に多重極展開を局所展開に変換するところでさまざまな方法が提案され
ている。いくつか見ていくことにしよう。
多重極展開から局所展開への変換は、結局展開係数のベクトルに変換行列
を掛ける操作になる。 Greengard と Rokhlin [11]
はこの変換を FFT を使って高速化できることを示唆した。 Elliott と Board
はこの方法を実現し、詳細な性能評価を行なった[10]。
例えば p=16までとったときには変換自体は3倍ほど高速化されている。ただ
し、上に述べた直接計算とのトレードオフがあるので、計算全体の高速化は
倍程度ということになる。また、通常は p=16 というようなところま
でとることはまずあり得ないので、実際の効果はさらに小さい。
FFT を使う方法は、理屈はともかく実装は非常に大変で、原論文を読むと苦労
がうかがえるのと同時にあまり自分でやりたいという気は起きなくなる。これ
に対し、「コロンブスの卵」的なアイディアで変換の計算量を 
から
に減らす方法も提案されている。この方法では、変換を一
度に行なうのではなく、まずz軸を移動の方向に持ってくるような回転を行
なう。すると、 
方向の次数 m が違うものはカップリングがなくな
るので、変換行列は変換前と変換後で m が同じ 
個以外はゼロになり、
計算量が減るわけである。変換前と変換後に回転する必要があるが、この回転
の計算量自体も 
である。したがって、変換は 
で済み、最適化す
れば 
となる。これは FFT よりも漸近的には計算量が大きい
が、実用的な p の値では同じ程度の効果があり、実装はずっと簡単である。
この他にも同じようなアイディアの方法がいくつも提案されているが、効能も 大差ないようなのでここでは省略する。
FMM を実装する上での大きな障害の一つは、定式化、特に球面調和関数のシフ トや多重極展開から局所展開への変換が複雑であり、プログラムの開発や高速 化が難しいということであった。
Anderson [3]は実装をずっと単純にする画期的な方法を提案 した。この方法を使うとプログラムが簡単になるだけでなく実行速度の観点か らもメリットがあるという結果が得られている[9]。以 下、この方法を簡単に紹介しよう。
基本的なアイディアは、球面調和関数の展開係数をデータとして使う代わりに 球面でのポテンシャルの値そのものを使おうというものである。半径 a の球の 中の粒子が外に及ぼすポテンシャルは、球の表面でのポテンシャルの値がわかっ ていればラプ ラス方程式の境界値問題を解けば与えられる。解はポア ソンの公式を使って以下のように書ける
ここで 
は n次ルジャンドル関数である。この積分を、球面上で
の適当な標本点を使った数値積分に置き換える。2次元であれば、円を等分割
すればいい。、3次元では、数値積分が球面調和関数の直交性を満たしている
必要がある。そのような置き方は良く研究されていて、文献
[23,19]を見れば数値的にそのような点が求めら
れているので、それを使えばよい。なお、技術的な注意としては、エリア
シングによる誤差を防ぐために数値積分で表現出来ている次数に応じて中の無
限和を適当なところで打ち切る必要がある。
子セルの外の多重極展開をその一つ上のセルでのそれに変換するには、上の式
に従って親セルを含む球面上でのポテンシャルを求めればいい。また、局所展
開も、
での展開が rでの展開に変わるだけで上と同じ式である。さら
に、多重極展開を局所展開に変換するには、上の式をつかって局所展開す
る球面上でポテンシャルを評価するればよい。つまり、上の式だけで、 FMM
に必要な数学がすべて表現されていることになる。これに対し、オリジナルの
FMM のための変換式は全部書くと 1 ページいっぱいに式が並ぶようなもので
あった。
Anderson の方法では別に計算量のオーダーが減るわけではなく、変換の計算
量は 
のままであるので理論的には前節で述べた方法に比べて良い点は
ない。しかし、式が簡単になっていることでプログラムの最適化やチューニン
グがずっと容易になるという現実的なメリットはけっして小さくない。
ここで、著者らが最近提案したP
M
法 (pseudo-particle multipole
method)[22,20] にも簡単に触れさせていただき
たい。基本的な考え方は Anderson の方法とほとんど同じものであるが、違
いは球面上でのポテンシャル分布を使うのではなく球面上の質量分布によって
多重極展開を表現するというものである。定式化は Anderson の方法と同様
に行なえ、
という形で実粒子の分布から仮想粒子の質量が表現される。ここで
,
はそれぞれ実粒子、仮想粒子の質量であり、 
は展開の中心から見て実粒子と仮想粒子がなす角度、
は実粒子の原点か
らの距離、 r は仮想粒子を置く球の半径、Kは仮想粒子の数である。
子セルの展開を親セルの展開にシフトするには、子セルの仮想粒子を実粒子で
あると思えば同じ式が使える。
局所展開については、Anderson の方法を使ってもよいし、上と同様に仮想 粒子を置くという形で表現することもできる。現在のところ、我々はツリー法 との組み合わせでしか実装していないので、こちらについては省略する。
普通のプログラマブルな計算機の上では、P
M
法には Anderson の方法
同様にプログラムが簡単であるということ以上になにかメリットがあるわけで
はない。しかし、専用ハードウェアとの組み合わせでは大きなメリットがある。
我々は GRAPE[26,24] という重力相互作用だけ
を高速に計算する専用計算機を開発しており、これと組み合わせると
P
M
法はセルからの重力も専用計算機で計算出来て画期的な高速化がで
きるのである。
ツリー法や FMM が提案されたのは 1980 年代中頃であり、大規模計算にはベ クトル計算機や場合によっては並列計算機が使われるようになったころである。 特に実装が簡単なツリー法は提案後すぐにベクトル化、並列化が精力的に研究 された。
まず、ベクトル化について簡単にまとめておこう。並列化についても、共有メ モリ型の機械では同じアイディアで行なえる。以下、アルゴリズムのなかでもっとも 計算量が多いツリーを再帰的にたどって粒子への重力を計算するところに話を 絞る。全体については文献[25]や原論文 [21,15,5]を参照されたい。
ベクトル化のためには、アルゴリズムの並列性を生かす必要がある。ツリーを たどるところでは2種の並列性がある。一つは別の粒子への重力は独立に計算 出来るという並列性であり、もう一つはある粒子への重力のうち、ツリーのな かで上下関係にないセルからの力は独立に計算できるというものである。
前者の並列性を使うには、粒子全部を一度に見るような形にプログラムを書き 換えればよい。これには、まず再帰によって表現されたアルゴリズムをスタッ クやポインタを使う明示的な繰り返しループに変形し、さらに一番外側の粒子 を順番に処理するというループを一番内側に持ってくる。これはかなり繁雑な 作業になるが、ツリー法の場合はもともとのアルゴリズムが単純であるので注 意すれば出来る。後者の並列性を使うには、やはり再帰を明示的な横型探索 (Breadth First Search) のループに書き換え、ツリーの1階層づつをベクトル 化されたループで処理していく。前者の方法のほうがベクトル長は長く、一般 には高い性能が出るようである。
分散メモリの並列計算機の場合には、アルゴリズムの並列性だけでは十分では なく、どのように空間分割をして粒子と領域をプロセッサに割り付けるかが問 題となる。逆にいえば、それが決まれば並列化そのものはそれほど複雑ではな い。メッセージパッシングの機械では Warren と Salmon による ORB（再帰的 直交分割）や Morton ordering による領域分割が良く知られている [27,28]。物理的に共有メモリである機械 や、分散共有メモリの場合は、Morton ordering によって粒子をソートしてお くことで自動的に空間分割が実現でき、かなり良い性能が得られている [18]。また、最近は HPF などのデータ並列型の言語を使っ た並列化も試みられている[16]。
前節で触れたが、専用計算機との組み合わせについてもここでまとめておこう。
筆者らのグループでは、GRAPE という名前で粒子間相互作用を計算することだ けに特化した計算機を開発し、実際に天文学の研究に使っている [24,26,17]。FMMやツリー法のよう な速い方法があるのに、そもそもこういった計算機を作ることに意味があるの だろうか。
もちろん、我々は意味があるから作っているつもりであるが、その意味は実は 2つあると考えている。一つは、粒子毎にタイムスケールが大きく違うような 問題では、ツリー法や FMM を使うのは非常に難しい、特に並列計算機で使え るような実現法はまだ存在しないということである。直接計算の場合には、多 重時間刻みに適した高精度な積分公式が知られていて、これは並列化や専用計 算機への実装も可能である。しかし、ツリー法や FMM の場合、特に並列化で は多数の粒子への重力を並列に計算できるということを利用しており、多重時 間刻みにするのは難しい。また、このような問題では要求される計算精度も高 いのが普通であり、直接計算に対する優位性がそもそも小さいのである。
もうひとつの理由は、ツリー法も FMM も、実は直接計算の部分を高速化する だけで非常に速くなるということである。例えばツリー法の場合、それほど計 算精度がいらない場合には、セルからの力はその重心からの力（モノポール） と思って良い。これはモノポールといっても実際には1次の項まで入っている し、セルが立方体なので高次の項の寄与は想像するより小さくなる傾向がある からである。この時には、粒子からの重力が計算出来ればよい。また、実際に 高次の項まで入れる場合でも、先に述べた疑似粒子多重極法を使えば粒子で表 現できる。
もっとも、順番にツリーをたどっていくような計算自体をハードウェアで実現 するのは、不可能ではないかもしれないが面倒である。我々は、Barnes によ るベクトル化の方法[5]を応用して、汎用計算機の方でツリー をたどるがその回数を減らすような工夫をして使っている。
専用計算機の効果であるが、以下に実例をあげてみよう。昨年度の Gordon-Bell 賞のコスト・パフォーマンス部門で、DEC Alpha のクラスター70 台を並列に使うという計算が、賞をとったわけではないが最終選考まで残った。 これは 533MHz の Alpha 70 台で並列ツリー法のプログラムを走らせ、実効 6.8 Gflops を達成したというものである。これに対し、我々は今年度の Gordon Bell 賞に昨年完成した GRAPE-5 を使った同じツリー法の計算で応募 した。我々の計算は、 DEC (Compaq) Alpha のホスト一台にGRAPE-5 ボード2 枚をつけた計算機で、昨年の Alpha 70 台とほぼ同じ実効速度を得た。（理論 ピーク性能も同程度である）コストパフォーマンスでは数倍良くなっているし、 占有スペースや消費電力ではほぼ2桁良い。これは一例であるが、問題にあま りよらないでこの程度の性能を得ることができ、専用計算機のメリットは非常 に大きい。