最初にC++でどのように数を表現したらいいかを考えてみました。ANSI C++について勉強しているうちに、ブルバキ風の抽象的な表現がある程度可能であることに気が付きました。これについては次の文書にまとめてみました。

「C++による代数的構造の表現」

　実際に整数型を実装してみたのですが、前節で考えたようにはいかないようです。数学上の抽象概念をC++の抽象クラスとして記述するというのは難しいということがわかってきました。難しいというよりもC++がそのようにできていないといった方がいいかもしれません。例えばある数クラスが順序環であるときに、順序環という代数構造に対応する抽象クラスがあって、そのクラスを継承することによって、順序環であることを表現したいところです。そして例えば順序環に対して定義できる関数があったときに、その関数の定義を、順序環という抽象クラスを継承するという性質だけを使って記述するのが理想的です。ところが現在のC++ではこういうことを直接的に記述することは残念ながらできないのです。

　最初の理想は実現できないことがわかりましたが、それでもなおC++にはCやFORTRANにはない魅力的な機能がいっぱいつまっています。いろいろ実験しながらC++による数値計算の可能性をさぐっていきたいと思います。

　ソースを以下に置いておきます。

• s0_base.h
• s0_math.h
• s0_math.cpp
• main.cpp

Visual C++ 6.0 + Windows 98という環境で試しましたが、その他の環境でも少しの修正で動かすことができると思います。

　実験の結果ですが、上のソースをVisual C++で最適化無しで、インライン展開もしないで実行すると数クラスでやった方が組込型の約２０倍程度時間がかかるとことがわかりました。次回は数クラスをチューンナップして、どれだけ組込型に迫れるか実験してみたいと思います。

　前節のソースをVisual C++で実行速度優先最適化及びインライン展開を行なってコンパイルリンクすると約２０倍程度であった実行時間の差が約１．５倍程度に短縮されました。だいたい予想通りの効果があることがわかりました。なお、この実験はかなり遅いマシンでやらないとはっきりした差がでないと思います。さらに組込型は32ビット整数演算で、数クラスはPentium内臓FPUの64ビット整数演算で行なっていますので、必ずしも条件が対等ではないことをお断りしておきます。次回は数クラスの機能を拡張して、より複雑な計算をしてみたいと思います。

　a,bを整数とし、b≠0としますと、次のような整数q,rが一意的に存在します。このことによって、q≡a/b, r≡a%bという２項演算を定義することができます。このような除法定理は、整数だけでなく例えば多項式環においても成り立ちます。一般に整域Rの0でない各元に対して、その大きさと呼ばれる整数d(a)≧0が定義され、0でない任意の元a,b∈Rについて、次のようなq,r∈Rが存在するとき、Rをユークリッド整域といいます。さらにq,r∈Rが一意的に定まるという条件のときに/と%が定義できるのですが、逆に/と%が定義できるユークリッド整域を代数構造としてどのように特徴付ければいいかは承知しておりません。

　ユークリッド整域において/と%が定義できるとき、任意の２元a,bについて、ユークリッドの互除法によりその最大公約数gcd(a,b)を計算することができます。gcd()をテンプレート関数として書いてみました。ソースを以下に置いておきます。

• s0_base.h
• s0_math.h
• s0_math.cpp
• main.cpp

gcd()においてその引数クラスが/と%が定義されているユークリッド整域という性質をもっていることを検査したいところですが、実行効率という観点から型検査はしないことにします。

　前節で最小公倍数を求める関数：gcd()を作りましたが、最小公倍数を求める関数をlcd(a,b)とすれば、という性質からlcd()はgcd()を使って簡単に作ることができます。

　/と%が定義できるユークリッド整域を簡単のために、仮に「一意ユークリッド整域」と呼ぶことにします。一意ユークリッド整域Rにおいて、自明でない２元一次不定方程式：はcがgcd(a,b)の倍数であるときに限り解が存在し、その解は拡張ユークリッド互除法によって求めることができます。これをeuc()という関数にしてみました。

• s0_base.h
• s0_math.h
• s0_math.cpp
• main.cpp

　文献[C1]を読んでいて、整数クラスにおいて、を定義しておくことによって、基本型:intが必要とされる文脈でこの整数クラスが使われる場合に暗黙の型変換が行なわれるという便利な機能があることに気がつきました。そこでクラス定義を修正しました。またメンバ関数名を小文字にするという修正も施しています。

　これでエラトステネスのふるいを書き直してみます。

• s0_base.h
• s0_math.h
• s0_math.cpp
• main.cpp

　よく使われる整数べき乗をipow()という関数にしました。実数べき乗については標準関数のpow()を多重定義することにします。べき乗はRSA暗号における暗号化及び復号化で使われるためにその高速算法が現在でも盛んに研究されているようです。ここでは汎用性があって簡単な繰り返し２乗法を文献[N3]を参考に実装しています。

• s0_base.h
• s0_math.h
• s0_math.cpp
• main.cpp

　そろそろ多倍長整数の実装に取り掛かりたいと思いますが、その評価用にLucasによるMersenne数の判定プログラムを書いておきます。現段階では演算のオーバーフローが検出できませんので、判定が本当に正しいかどうかを知るにはオーバーフローの検出機能を付け加える必要があります。

　また整数クラスが２進表現のときにbit毎の論理演算による高速算法を記述できるように組込型に対するbit毎の論理演算と同じ演算子関数を実装しています。

• s0_base.h
• s0_math.h
• s0_math.cpp
• main.cpp

　試行錯誤の末、多倍長整数クラスを次のようにするのがいいのではないかという結論になりました。まず基数Pを定めて、任意の正整数aをとP進表現します。このときのsをaの次数と呼ぶことにします。多倍長整数クラスはPを定数として定め、配列の長さを可変とするテンプレートクラスとして表現します。

UINT32 status;
  UINT32 degree;
  UINT32 v[n];

degreeに次数を、v[i]にa_iを格納します。statusは後にPentium用に手書きで最適化することも考えて、Pentiumのフラグレジスタと同じbit配列にしました。現在はすべてC++で記述しているので、使っているのはOF(オーバーフロー),ZF(ゼロ),SF(負)の３bitです。a=0のときにはdegree=0と定めることにします。

　簡単な方法(文献[N5],[N6])で演算を実装して、Mersenne数の判定をしてみました。

• s0_base.h
• s0_uint32.h
• main.cpp

とりあえずはうまく動いているようですが、他の方法では検証しておりませんので、不具合があるかもしれません。Pentium-133MHzでp=1279のときに２４秒もかかってしまいます。UBASICでは１秒もかかりませんので話になりません。

　次回は高速乗算法を少し研究してみたいと思います。

　２つの多倍長整数a,b:の積abを計算することについて考えます。まず基数Pを変数Xに置き換えることによって多項式a(X),b(X)を作ります。この積c(X)≡a(X)b(X)は高々2s-2次ですが、長さを2sに揃えるためにと書くことにします。このとき数列(c₀,…,c_2s-1)は、長さ2sの数列(a₀,…,a_s-1,0,…,0),　(b₀,…,b_s-1,0,…,0)の畳み込みになっていますが、次のことが成り立ちます。

　この証明については、例えば文献[S1]を参照してください。これを利用して多倍長整数の乗算を行うことにします。その手順を簡単にまとめておきます。

　これを実装することが目標ですが、まずはDFTの高速算法であるFFT(高速フーリエ変換)について実験してみようと思います。

　多項式の高速乗算法を応用して多倍長整数の高速乗算法を実装しようという予定だったのですが、T.W.ケルナーの本に書いてある「このやり方自体かなり手間が省けるとはいえ、エレガントとはいえない。」という言葉が気になって、ストラッセン・ショーンハーゲ法に挑戦する気になってきました。これについては、[N5],[N6]にも大雑把な説明があるのですが、やはりクヌースの本[S2]を熟読しないと独自に実装できるレベルまで理解することは難しいようです。

　クヌースの本を読んでわかったのですが、n桁の乗算を計算量O(nlog nlog log n)で行なうというストラッセン・ショーンハーゲ法は、理論上の計算量の上界を改善するためのものであって実装を目指したものではないということです。クヌース自身も、実際には別の方法がいいということを述べています。したがってクヌースの本においても具体的なストラッセン・ショーンハーゲ法の紹介はされていなくて、論文が示してあるだけでした。和田秀男さんの[N5],[N6]には、ストラッセン・ショーンハーゲ法の概要が説明されていて、実装するときはクヌースの本を参考にしなさいということが書いてあるのですが、どうもこの意味を勘違いしていたのかもしれません。

　ストラッセン・ショーンハーゲ法では、計算する数の大きさによって語長が変わっていきます。ですからそのまま作譜しようとすると演算が多倍長演算になってしまいます。また2ⁿ+1を法としているために効率の良い最適化ができません。このように不利な状況が多いので、実装しようかどうか迷っています。ただ、すべて整数演算のみで、フーリエ変換がシフト演算と加法のみでできるという点に非常に魅力を感じています。特に、整数の乗算のために浮動小数点演算を使うというのがどうも引っかかっています。整数の演算は整数演算命令だけでやってみたいということもあります。

　このようなわけで夏休みに、ストラッセン・ショーンハーゲ法を効率よく実装する方法を考えることにしました。

　ストラッセン・ショーンハーゲ法を実装するというのは、長手順の詰め将棋のような感じで、頭の中ではすべてを読みきったつもりでいるのですが、いざこれを文章で説明したり、作譜しようとするとかなりの情報量になってしまいそうです。そこで、論点を整理して見通しよくするために離散Fourier変換を抽象化することについては、別にひとつの文書としてまとめてみました。[F1]

　さて「ストラッセン・ショーンハーゲ法を効率よく実装する方法」なのですが、効率がいいかどうかは別にして、少なくとも実装する方法についてはなんとかなりそうです。たばになったメモをめくりながらこれをまとめている途中です。完成予定は９月中というところでしょうか。多倍長整数クラスをUBASIC並みの性能にするのは、何とか今世紀中に達成することを目標にしています。