2012/12/03

TeXを電卓として使おう(あるいは、TeXでベキ乗根)

プログラミング言語 TeX にはインタプリタがついてきます(ツッコミ禁止)。そこで他の多くの言語のインタプリタのように、TeX インタプリタを「電卓」として使ってみましょう。(なお、この記事は TeX & LaTeX Advent Calendar のために書いたものです。明日は @neruko3114 さん!)

まずは TeX インタプリタの起動です。日本語の組版をするなら ptex を実行するところですが、電卓に日本語はいらないので、今回は世界標準らしい pdftex を使います。 TeX インタプリタの標準プロンプトは * (アスタリスク)ですが、起動時だけは興奮気味に ** と2つ重ねて出力されるので、お約束として一言目は毎回 \relax と入力して TeX インタプリタを落ち着けてください。

$ pdftex
This is pdfTeXk, Version 3.141592-1.40.3 (Web2C 7.5.6)
 %&-line parsing enabled.
** \relax
*

さっそく 1+2 と入力してリターンを押したいところですが、TeX で単に数字や演算子を書いてしまうと、通常の文脈ではドキュメントへと印字する文字列とみなされます。いかにもドキュメント処理システムっぽいですね。数字を整数値として扱うには、整数用の \count レジスタというハコを用意してそこに値を格納し、そのレジスタ上で演算を行います。整数の加算からやってみましょう。

* \newcount\a
* \a=1
* \advance\a by 2
* \showthe\a
> 3.
<*> \showthe\a

?      % ここで単にリターン

*

1行めの \newcount\a で新しい \count レジスタ \a を用意し、2行めで\a に整数値の1を代入、 \advance\a by 2 という命令で \a 「2」の値に2を加算し、 \showthe というコマンドで \a の中身を出力しています。結果は赤字の「3」です。 \showthe で結果を表示したあとはプロンプトが ? に変わりますが、気にせずリターンを押せば元の * に戻ります。簡単ですね。なお \advance コマンドの引く数にマイナス記号をつければ(\advance\a by -2)、減算ができます。減算専用のプリミティブは用意されていません。

続いて、整数の乗算と除算です。それぞれ \multiply\divide という命令が用意されています。先程のセッションに続けて以下のように入力してみましょう(したがって、いま \a には 3 が入ってます)。

* \multiply\a by 3
* \showthe\a
> 9.
<*> \showthe\a

?      % ここで単にリターン

* \divide\a by 2
* \showthe\a
> 4.

最後の行は「9÷2」の結果なので「4.5」になるはずですが、 \count は整数用のハコなので端数は切り捨てられます。キャストとかされません。というか、そもそも実数そのものを扱うハコがありません。

では、結果が実数になるような複雑な計算は実現できないのでしょうか。試しにベキ乗根のつもりで \sqrt\a とかしてみましょう。

* \sqrt\a
! Missing $ inserted.

$がない、つまり数式組版モードじゃないと言われました! 当たり前と言えば当たり前なのですが、TeX というのは本来はドキュメントシステムなので、 sqrt はもちろん、 + や * 、 さらには数字にいたるまで、ふつうの言語なら数値やその演算を提供するであろう要素はほとんど組版の機能に利用されてしまっています。

しかし TeX は、「かっこいい出力ができる汎用プログラミング言語」でもあるはずです。ベキ乗根の近似値くらい計算できないはずがない。そこで、ベキ乗根の近似値をかっこよく組版で使うための TeX マクロ \sqrtd を作ってみましょう。目標は、こんな入力から、

You can get \sqrtd{2}, \sqrtd{2.5}, \sqrtd{3} and \sqrtd{3.14} 
by Newton's method with initial guess $1.0$.

こんな出力を得ることです。

newtonmethod-output

(2のベキ乗根があまり精度よくないのは目をつぶってください><)

TeXでニュートン法を実装しよう!

ベキ乗根の近似値を得る方法として一番簡単なのはニュートン法でしょう。まずはニュートン法を使って値 x のベキ乗根を求める手順のおさらいです。

  1. ベキ乗根になりそうな値 y を適当に推測します。
  2. その推測値 y のベキ乗 y2 と、元の値 x との差を調べます。
    • この差が十分に小さかったら、その推測値 y はだいたいベキ乗根に近いということで、答は y です。
    • この差が十分に小さくなかったら、推測値 y と x/y の真ん中の値をあらためて推測値 y とし、1.に戻ります。

いきなり TeX はきついので、とりあえず Scheme で書いてみます(詳しくはSICPの1.1.7節を参照)。

(define (sqrt-iter guess x)
  (if (good-enough? guess x)
      guess
      (sqrt-iter (improve guess x)
                 x)))

(define (good-enough? guess x)
  (< (abs (- (square guess) x)) 0.001))

(define (improve guess x)
  (/ (+ guess (/ x guess))) 2))

このアルゴリズムを TeX で実装したいわけですが、 TeX にはそのものずばりな実数値がありません。ただ、TeX本来の機能である自動組版を制御するために、値が実数になる「寸法(dimension)」というデータ型は用意されています。このデータ型上にも、整数値のときと同じ演算命令(\advance\multiply\divide)が制限付きながら用意されているので、そこに上記のアルゴリズムの x や y を格納すれば原理的にはベキ乗根が求められるだろう、という算段です。

ここで厄介なのが、寸法上の演算の制限です。寸法の値を自由に掛け算したり割り算したりすることはできません。寸法に対して可能な乗算と除算は、具体的には以下のとおり。

  1. 寸法は、整数倍したり、整数で割ったりすることはできる。
  2. 寸法の前に、何か数字っぽいものを書くと、その数字の値が寸法の値に掛けられる。

上記の 1. は、寸法の値に実数を掛けたり割ったりすることは基本的にできない、つまり、実数どうしの乗算と除算が寸法型のままではできないということです。

上記の 2. の手段を使えば、寸法に実数を掛けることはできます。ここでちょっと実際の TeX コードを見てみましょう。寸法データ型を用意して「3.14 ポイント」を格納し、その寸法を 1.5倍するには、例えばこう書きます。

\newdimen\x      % 新しい寸法データ\xを用意
\x=3.14 pt       % 寸法\xに単位付きの実数を格納
\x=1.5\x         % 寸法\xを1.5倍

寸法の単位はいくつかありますが、ここではポイント(pt)を使っています。2行めでポイントの単位「pt」を省略して 3.14 という実数値を寸法に代入しようとしてはいけません(TeX に苦情を言われます)。3行めの結果、 \x には「4.71 ポイント」に相当する値が入ります。

ここで、実数に実数を掛けようとして次のように書いたらどうなるか考えてみましょう。

\newdimen\x      % 新しい寸法データ\xを用意
\newdimen\y      % 新しい寸法データ\yを用意
\x=3.14 pt       % 寸法\xに単位付きの実数を格納
\y=1.5 pt        % 寸法\yに単位付きの実数を格納
\x=\number\y\x   % 寸法\xを1.5倍?

最後の行に出てくる \number は、直後のデータ型の数字表現を返すコマンドだと考えてください。上記を TeX に渡すと「Dimension too large.」というエラーになり、 \x16383.99998pt という寸法になるはずです。これは、 TeX の寸法の上限(をポイント単位で表したもの)です。

何が起こったのでしょう。\number\y は 1.5 ではなく、「ポイントで表すには "too large" だ」と TeX に言われてしまう何か巨大な寸法であり、TeX はとりあえずの処置としてポイント表示で最大な数字を \x に入れたのです。ちなみにその巨大な寸法とは、 1.5 ポイントを sp という単位で表した数字に 3.14 ポイントを掛けた寸法であり、308674.56pt です。実をいうと、TeX は内部では sp という単位の整数倍で寸法を扱っています。だから、寸法で実数が扱えるというのも実はまやかしで、ポイントくらい大雑把な単位で見れば小数表示になる寸法を実数とみなして扱える、ということだったりします。

というわけで、寸法そのものは単位付きの値なので、その単位で表した場合の実数値を示す数字として使うには、何とかして単位の部分をひっぺがす必要があります。そのための仕組みがこちらの \strippt です。

\newdimen\zero \zero=0pt
\begingroup
  \catcode`P=12
  \catcode`T=12
  \lowercase{
    \def\x{\def\rempt##1.##2PT{##1\ifnum##2>\zero.##2\fi}}}
  \expandafter\endgroup\x

\def\strippt{\expandafter\rempt\the}

何を言ってるのか分からないと思うが、私もこんなのゼロから書けるわけありません。これは LaTeX のコアで定義されているマクロをほぼそのまま拝借してきたものです。 \expandafter のとこだけ説明しておくと、ポイントを表す「pt」の各文字をちょっと特殊な文字とみなす閉じた環境で「なんとか.なんとかpt」というトークン列にマッチさせて \rempt というマクロを定義したいのだけど、それを外から使えるようにするために環境を閉じる \endgroup の展開を遅らせているのが最初の \expandafter で、その\remptを「なんとか.なんとかpt」という形の引数に対して使うべく \the で先に展開するために配置されているのが 2つめの \expandafter です。

以上の下準備でようやく、ニュートン法を実装するのに必要なベキ乗のためのマクロ \square が定義できるようになりました。

\newdimen\tempdimen \tempdimen=\zero

\def\square#1{\tempdimen=#1%
  #1=\strippt \tempdimen \tempdimen
  \tempdimen=\zero}

この \square を使ってニュートン法のアルゴリズム全体を書き下してみます。

\newdimen\prevguess
\newdimen\val
\newdimen\guess

\def\sqrtiter{%
  \prevguess=\guess
  \square\prevguess
  \advance\prevguess by -\val  % 推測値のベキ乗と元の値との差
  \ifdim\zero>\prevguess       % 差の絶対値
    \prevguess=-\prevguess\fi
  \ifdim\prevguess<0.001pt     % 差が十分小さかったら、
    \let\next=\relax           % おしまい。
  \else                        % 十分小さくなかったら、
    \improve\guess\val         % 推測値を改善して、
    \let\next=\sqrtiter        % 再帰。
  \fi\next}

上記でまだ未定義なのは、ニュートン法のアルゴリズムで推測値を改善させる部分を担うマクロ \improve です。次はこれを実装します。\improve には「実数÷実数」の演算がでてきますが、寸法の値を使った除算には直接の手段はなさそうなので、やや面倒な工夫が必要です。全体はこんな定義になりました。

\newcount\tempcountA
\newcount\tempcountB
\newdimen\tempdimenA
\newdimen\tempdimenB

\def\improve#1#2{%
  \tempcountA=#2%                   % いったん割られる数を整数にして、
  \multiply\tempcountA by 1000%     % 十分な倍率で大きくしておく。
  \tempcountB=#1%                   % 割る数も整数にする。
  \tempdimenB=#1%
  \divide\tempcountA by \tempcountB%% ここで割り算。
  \tempdimenA=\tempcountA pt%       % 結果をポイント単位で寸法に戻す。
  \divide\tempdimenA by 1000%       % 倍率を戻す。
  \advance\tempdimenA by \tempdimenB% 
  \divide\tempdimenA by 2%
  #1=\tempdimenA}

\newcount は、整数を格納できるハコを作ります。そこに「寸法」を代入すると、先に説明したsp単位の整数に変換されて格納されます。整数どうしは \divide というプリミティブで割り算できるので、これだけでもよさそうに思えますが、 \divide は剰余を切り捨てて結果を整数に丸めてしまうので、これだけだと結果の誤差が大きくなりすぎます。そこで、割られる数だけ十分な倍率を掛けておいて、割り算して結果を寸法に変換し直したあとで倍率を戻しました。

\improve の実装に必要な残りの演算は、単純な足し算と整数2による除算なので、特にこれ以上の工夫はいりません。これで寸法 \val\guess にそれぞれベキ乗根を求めたい数と初期推測値を入れておいて \sqrtiter を起動すると、ベキ乗根の値にポイントの単位がついたものが \guess に代入されます。あとは、こんなコマンドを定義すれば目標達成です。

\def\sqrtd#1{\val=#1 pt%
  \guess=1.0pt%
  \sqrtiter%
  $\sqrt{\number#1} \approx \strippt \guess$%
}

You can get \sqrtd{2}, \sqrtd{2.5}, \sqrtd{3} and \sqrtd{3.14} 
by Newton's method with initial guess $1.0$!

今回書いた TeX スクリプトの全体はこちら。なお、 LaTeX で実数を本気で使いたいなら、こんな面倒なことをしなくても、fpパッケージという便利な浮動小数点数演算のための仕組みがあります。そちらを使いましょう。

No comments :