Befunge-93で常用対数を求めた話

初投稿です。タイトル通り常用対数を求めていきます。

常用対数のおさらい

常用対数について(高校数学のおさらい)

$10^a =b$を満たす$a$のことを$b$の常用対数といい、$b$を真数といいます。
この時、$a=\log_{10}b$と表記します。
また、常用対数にはいくつかの公式があります。

• $\log_{10}10=1$
• $\log_{10}ab=\log_{10}a+\log_{10}b$
• $\log_{10}a^k=k\log_{10}a$

この公式を今回のアルゴリズムに利用していきます。

Befungeで対数を求めるアルゴリズム

$\log_{10}a=\frac{1}{2}\log_{10}a^2,\log_{10}10a=1+\log_{10}a$
この2つの公式を用いて対数を求めます。
なお、$10^m\approx a^n\Rightarrow\log_{10}a\approx\frac{m}{n}$は精度と速度の面から採用を見送られました。

具体的に書くと、

• 初期値として対数の近似値$=0$、真数$=a$とする
• もし$a<10$なら、$a$を$2$乗して、係数に$\frac{1}{2}$をかける
• もし$a>=10$なら、$\log_{10}a$の係数を対数の近似値に加算し、$a$を$10$で割る

これをそのまま実装できればいいのですが、言語仕様上複数の問題点があります。
以下、一つ一つ問題点とその対策について書いていきます。

1.Befungeが小数に対応していないのに10で割らないといけない問題

仕様上割り算は整数の範囲でしか行えないため、小数を記述することができません。
これについては、真数を$\frac{a^b}{10^k}$の形で記述することで、分母に$10$をかけるという動作で解決できます。

2.桁が大きすぎてオーバーフローしてしまう問題

$1$回の操作で$2$乗するので、$n$回の操作では元の数が$2^n$倍になることは容易に想像できます。
するとスタックの値がInfinityとなり、正しく計算できません。
これについては、ある一定以下の桁数を四捨五入することによって、解消することができます。
当然処理する桁を増やすほど精度は上がります。が、$20$桁もあればそれなりの精度は出せます。

3.係数の分母が指数関数的に増加する問題

分母が$2$のべき乗になるので、計算をすればするほど分母が大きくなってしまいます。
これについては、分子が常に奇数になることが証明済みなので計算回数を少なめに押さえればいいです。
分母が整数に収まる範囲として、$2^{64}$くらいまでの計算ができるので、かなりの精度は確保できます。

実装方法

有効桁数を指定する

この値より真数が大きければ四捨五入するという、ある種のしきい値になります。

"."145*    v
#v_1-\25**\>:!
 @

まずスタックに出力用の小数点と$1$と$20$を積み、スタックの一番上をコピーしてnot演算をします。
すると行目左の_で右に行くようになります(スタックの最後=指数が$0$なら左に行きます)。
このように、条件を満たしたら右に、そうでなければ左に、ループを抜けたら下の行にカーソルが移動するようにプログラムを構築すると、フローチャートと似たカーソルの動き方となりコーディングがしやすいです(あくまで個人の感想です)。
右に行くとスタックを入れ替え$10$をかけ入れ$1$を引きます。しきい値が$10$倍され指数が$1$減ります。

真数を指定する

v>6\p55+
<  v12345
**+>\:9`#v_\25
         @

スタックの$1$番上に$0$があるので、$6$を積んで入れ替えることで値を有効活用できます。
vの右側には対数を求めたい真数を書きます。この場合、$\log_{10}12345$を求めることになります¹。
55+で$10$を積んでいるので、$9$より$2$番目のスタックが大きいかを比較することで、
それが入力の一部なのか否かを判定し、入力の一部ならスタックの1番上に$10$をかけ2数を足します。
2行目の実行順序は左方向なので、最上位の桁が一番上に積まれます。
すると$123=3+12\times10$のように一番上のスタックに$10$をかけ2数を足すと
もとの自然数が再帰的に求まります。
なおこの構造上、12345のところに算用数字以外を入力すると不都合な動作が生じ得ます²。

真数の条件判定と加工

vp001p01$<>"RORRE",,,,,@
>10g:0`!#^_:25*%#v_25*/10p
00g25**00p       >00g10g`#v_
 v`1g01p05*56p040p030p020 <
v_$:,2+,@
@

ここでは真数の大きさを判定しています。
左上の部分で判定に用いた$10$を消去し得られた真数を$[1,0]$に、判定用に$1$を$[0,0]$に格納します。
真数は$[1,0]$の場所に格納されているので取り出し、$0$以下ならERRORと出力して停止します。
そうでなければ$2$行目で右に行くので、$10$の倍数か判定します。
$10$の倍数なら精度確保のため$10$で割ります³。
$3$行目では$[0,0]$の中身が初めて真数を超えるまでを$10$倍しています。
こうすることで、この数と真数をそれぞれ$2$乗して桁数を判定できるようになるため便利です。
こうして$10$の倍数でない自然数と桁数が求まったら各種の値を$4$行目で記録します。
$[2,0]$に係数の指数、$[3,0]$に近似値の分子、$[4,0]$に近似値の分母、$[5,0]$で加算回数を指定しました。
$5$行目で真数が$1$かどうか判定し$1$なら対数は整数なので、.0と出力して停止します。
.のアスキーコードは$46$,0のアスキーコードは$48$なのでコピーして$2$を足すことで0を生成しています。

公式から近似値を求める

>>00g:10g`!#v_:*25*/00p10g:*10p20g1+2v
_^#        `g06g01_v#!`*:/*25g06g01p0<0/\g00p01/*25+**25`4%*25:/\g01<\*::*25
  v-g04g02 $<      >60g25*/:25*/                                    ^
  >:!#v_1-30g2*30p
p|    >30g1+30p20g40p00g25**00p50g1-:50
 @

$1$行目で真数を取り出し桁数参考値と比較し、真数のほうが小さければ右に行って
$\log_{10}a=\frac{1}{2}\log_{10}a^2$を使って式変形して下に降ります。
$2$行目では桁あふれの検知と四捨五入による打ち切りの処理を行っています。
概ね桁数はしきい値の桁数より$2$少ない⁴ので、それを考慮して$3$行目右側で効率的に計算しています。
$10$で割ったあまりが$4$より大きいならば$10$を加算してから$10$で割り、そうでなければそのまま$10$で割ります。
$a>10$判定用の数は$1$の位が$0$なのでそのまま$10$で割ります。
$1$行目で下に行った場合、余分な値をgでまとめて$で削除し、4行目で近似値を通分します。
さらに$5$行目に降り分数の加算を行い残り加算回数が$0$になるか判定します。
$0$なら下に、そうでなければ上に行き、また近似値計算を繰り返します。

分数を小数に変換し出力

v>$,30g25**30p150p40g
>:!#v_1-50g2*50p
30pv>50g2*50p30g2*1+
#@_>30g:50g/86*+,50g%25**30p30g!

空いた$[5,0]$(加算回数でここまでで$0$になっている)を利用して分母を復元します。
1行目で初期値を設定し残り乗算回数である$[4,0]$の値をスタックに積みます。
その値が$0$になるまで、分母に$2$をかけ続けます。
$3$行目で少々細工をしています。最下位ビットより更に下の桁をが$0$の場合と$1$の場合で
上下から評価しその平均値を取ることで、微小なりとも精度を上げています。
分子が$0$になるまで、$4$行目で筆算式に割っています。今回のアルゴリズムでは分母の素因数は$2$のみ
なので、分数は必ず有限小数になります。そのため、割り算の打ち切り条件をあまり$0$としました。

完成品

最終的に以下のプログラムが完成します。

"."145*    v
#v_1-\25**\>:!
v>6\p55+
<  v2
**+>\:9`#v_\25
vp001p01$<>"RORRE",,,,,@
>10g:0`! #^_:25*%#v_25*/10p
00g25**00p        >00g10g`#v_
 v`1g01p05**147p040p030p020<
v_$:,2+,@
>>00g:10g`!#v_:*25*/00p10g:*10p20g1+2v
_^#        `g06g01_v#!`*:/*25g06g01p0<0/\g00p01/*25+**25`4%*25:/\g01<\*::*25
  v-g04g02 $<      >60g25*/:25*/                                    ^
  >:!#v_1-30g2*30p
p|    >30g1+30p20g40p00g25**00p50g1-:50
v>$,30g25**30p150p40g
>:!#v_1-50g2*50p
30pv>50g2*50p30g2*1+
#@_>30g:50g/86*+,50g%25**30p30g!

$4$行目のvより右に整数を記述することでその対数を求められます。

出力例

こちらのサイトに貼り付け⁵て実行。

正確な値がおよそ$0.30102999566398119...$なのを考えると、やや正確⁶な値が出ています。

---

1. 正確には、求める部分は小数部分のみです。常用対数表のようなイメージです。 ↩

2. 不都合な動作の例として、小数点が正しく表示されないことなどが考えられます。 ↩

3. $\log_{10}10a=\log_{10}10+\log_{10}a=1+\log_{10}a$となり小数部分に影響を与えないため、
   2乗したときの誤差を小さくするため10で予め割っておきます。
   今回は小数部分のみ求めるため整数部分の考慮をする必要がありません。 ↩

4. ループ数にもよりますが計算を繰り返すとしきい値より1桁または2桁小さい状態になります。 ↩

5. Ctrl+Vだと書式が保たれたまま貼り付けられで実行不可能となりますが、
   Ctrl+Shift+V(プレーンテキストとして貼り付け)を行うことで対応できます。
   1行目に謎の改行が発生することもありますが、これはページをリロードすると解決する場合が多いです。 ↩

6. どれほどの精度を求めるかによりますが、2のべき乗の桁数を求める程度なら無視できる誤差です。 ↩