C++で有理係数多項式

けっこう探した挙げ句にとうとう見つけた.日本語の情報は皆無だが,少なくともちょっとしたサンプルを書くにはなかなか使い勝手がよい.

トップページから辿れる ftp://ftpthep.physik.uni-mainz.de/pub/gnu/cln-1.1.9.tar.bz2 を落として展開し,configure.私はすでにaptでGMPを入れていたので ./configure --with-gmp としたが,GMPは必須ではない(GMPは手書きのasmコードを使っており,使うと速くなる).

configure後はもちろんコンパイル.私は make >& makelog.1 & とかするのが好みだ(途中で出力が気になったら tail makelog.1tail -f makelog.1 で覗けばよい).なおコンパイルには40分以上かかった(他のプロセスを動かしてたり,私のマシンがショボいのもあるが).ついついpstreeとかtopとかしちゃうよね(笑).終ったら make check して sudo make install.configure時に --prefix を指定しなかったので /usr/local 下にインストールされる(/usr/local/lib/libcln.so とか).これが嫌なら ./configure --prefix=$HOME/local とかにすること(学校で使う人とか).
インストールが終ったら適当なディレクトリを掘って(私の場合は ~/prog/test/cln/)そこにサンプルコードを色々置く.1ファイル1プログラムの単純なものをたくさん置くときにはGNU Makeのワイルドカード機能が便利だ:

targets = $(basename $(wildcard *.cpp))

all : $(targets)

.PHONY : all clean

CXXFLAGS = -Wall
LDFLAGS = -lcln -lgmp

% : %.cpp
        $(CXX) $(CXXFLAGS) $(CPPFLAGS) $< -o $@ $(LDFLAGS)

clean :
        -rm *~ $(targets)
        ctags -R

これを使えば,例えば foo.cpp には g++ -Wall foo.cpp -o foo -lcln -lgmp が実行される.詳しくは http://www.ecoop.net/coop/ から辿れるGNU Make 3.77マニュアル邦訳とかを参照.

なおCLNのインストール時に --prefix で特殊な場所を選んだ人は,LDFLAGS = `cln-config --cppflags --libs` とすること.外部プロセスの出力を取り込むための逆クォートに注意.cln-configはpkgconfigにならった補助プログラムで,PREFIXに合わせて適切に -I/usr/local/include -L/usr/local/lib -lcln -lgmp などのオプションを出力してくれる.なおインデントはもちろん8スペースじゃなくてタブね*1
で,肝心なのはどういうプログラムを書くかということだが,試しに有理係数多項式のGCDを求めるプログラムを書いてみたのでそっちを見てほしい.M. Bronstein, Symbolic Integration I ---Transcendental Functions, Springer の疑似コードをそのまま実装してみたものだ.オブジェクトの値渡しを多用したりしてるがソースがキモいとか言わないように.

コンパイルして実行するとこんな感じになる.

[nodakai@sprawl cln]$ ./chap1
 --- Section 1-2 ---
f = 3*x^3 + 1*x^2 + 1*x + 5
g = 5*x^2 + -3*x + 1
f = (3/5*x + 14/25) * g + 52/25*x + 111/25

f = 3*x^3 + 1*x^2 + 1*x + 5
g = 5*x^2 + -3*x + 1
25 * f = (15*x + 14) * g + 52*x + 111

 --- Section 1-3 ---
f = 1*x^4 + -2*x^3 + -6*x^2 + 12*x + 15
g = 1*x^3 + 1*x^2 + -4*x + -4
gcd(f, g) = 5*x + 5

f = 1*x^4 + -2*x^3 + -6*x^2 + 12*x + 15
g = 1*x^3 + 1*x^2 + -4*x + -4
gcd(f, g) = 5*x + 5 = (-1*x + 3)*f + (1*x^2 + -6*x + 10)*g

f = 1*x^4 + -2*x^3 + -6*x^2 + 12*x + 15
g = 1*x^3 + 1*x^2 + -4*x + -4
gcd(f, g) = 5*x + 5 = (-1*x + 3)*f + (1*x^2 + -6*x + 10)*g

もちろんきれいな出力はプログラマが自力で書式を定めて実装するしかないので期待しないで欲しい.それとシンボリックな構文木を扱う能力も当然ながらない(上の 25 * f とかいうのは手書きで出力したものだ).まぁ数学系のアルゴリズムを試すには十分便利だろう.

*1:GNU Make 3.80では「タブのつもりですか?」という有名なメッセージはもう修正されてしまっている.こないだこの話題で盛り上がっていつ修正されたかをCVSで調べたんだが,もう忘れた...