2011年1月11日火曜日

数式処理システム Maxima を使ってみた一例

少し前に「独学Linux」さんで数式処理システムの Maxima が取り上げられていて、 twitter のフォロアーさんの一人である @labunix さんがツイートしたものに思わず反応してしまいました。


私の場合、 Maxima はたまに使う程度であるのですが、うまくキマるとなかなか楽しいのでそんな例(のようなもの?)を紹介したいと思います。

今回紹介するのは Aizu Onlline Judge の中からで AOJ 0081 の問題です。問題の概要としては
平面上の異なる3点 P1, P2, Q の座標の組をファイルから読み込んで、それぞれの組について、点P1 点P2 を通る直線を対称軸として点 Q と線対称の位置にある点 R を出力して終了するプログラムを作成してください。

というものです( 詳細については AOJ 0081 の問題文を参照してください )。問題そのものはそんなに難しいものではなく、物理的に紙の上でであればコンパスがあれば中学生でも解けるものですが、具体的に座標を求めるとなると 点Qを通る直線P1P2上の直交点を求めてQの反対側に伸ばすということを考えるのが正攻法なのではないかと思います。

ですが、手順が増えるとバグを埋め込んでしまったり例外条件などを考慮しないといけないし、よく考えると、コンパスを使って点Rを特定するように線分 P1QとP1R、P2QとP2Rの長さが等しいということを使ってダイレクトに2つの変数を2つの条件式から特定するということができそうな気がします。とは言っても手で計算するのはなかなか大変そうだということで Maxima の登場です(あえて wxMaxima のスクリーンショットを載せましたが、実際には Terminal 上で maxima コマンド使いました )。


ここで、[x = xq, y = yq] というのは点Qそのものを示す自明な解ですので、そうではない他方のほうが点Rを示す座標ということになります。あとは、写経(苦笑)のごとく、モリモリ実装すれば完成です。
x = ((x2-x1)*y2+(x1-x2)*y1)*2*yq +
    (2*x1-xq)*y2*y2 + 2*(xq-x2-x1)*y1*y2 + (2*x2-xq)*y1*y1 +
    (x2*x2 - 2*x1*x2 + x1*x1)*xq;
x/= y2*y2 - 2*y1*y2 + y1*y1 + x2*x2 - 2*x1*x2 + x1*x1;

y = ( y2*y2 - 2*y1*y2 + y1*y1 - x2*x2 + 2*x1*x2 - x1*x1 ) * yq +
    ( (x2-x1)*xq - x1*x2 + x1*x1 ) * 2 * y2 +
    ( (x1-x2)*xq - x1*x2 + x2*x2 ) * 2 * y1;
y/= y2*y2 - 2*y1*y2 + y1*y1 + x2*x2 - 2*x1*x2 + x1*x1;
と Maxima の活用事例というより数式写経ねたを書いてきましたが、この問題は「パソコン甲子園 2005」本戦で出題された問題で、レギュレーションを読む限りではこの大会では Maxima を使用することは許されていないようです。実際には正攻法で点Qを通る直交点を計算して云々という実装をすることになるのでしょう。