はじめに

三角関数が必要な以下の問題を Befunge-93 で解きました。

Befunge-93 とは？という方は前回の記事もご参照ください。

問題概要

点 $(a, b)$ を原点を中心として反時計回りに $d$ 度回転させた点 $(a', b')$ の座標を求めよ。
$-1000 \le a, b \le 1000$
$1 \le d \le 360$

通常の解法

問題自体は座標の回転の公式に当てはめるだけで解けます。

\displaylines{
a' = a\cos(d°) - b\sin(d°) \\
b' = a\sin(d°) + b\cos(d°)
}

しかし、Befunge には $\sin$ $\cos$ なんて高等な関数もなければ、小数の演算も用意されていません。整数の演算のみでこの三角関数の計算を実現する必要があります。

CORDIC アルゴリズム

整数の演算のみで三角関数の有理数近似を求めるアルゴリズムとして、CORDIC アルゴリズムがあります。

ここではざっくりした説明を書いていきます。

点 $(x, y)$ を $\pm\theta$ 回転させる操作を考えると、

\displaylines{
x \gets x\cos\theta \mp y\sin\theta = \cos\theta(x \mp y\tan\theta) \\
y \gets y\cos\theta \pm x\sin\theta = \cos\theta(y \pm x\tan\theta)
}

と変形できます。（ただし $\cos\theta \not = 0$ ）

この式から、$\pm\theta_0,\pm\theta_1,\dots,\pm\theta_N$ と順に回転させる操作は、まず各 $\theta_i$ について

(x, y) \gets (x \mp y\tan\theta_i, y \pm x \tan \theta_i)

と更新し、$\cos\theta_i$ 倍する部分を最後にまとめて

(x\times\prod_{i=0}^N\cos\theta_i, y\times\prod_{i=0}^N\cos\theta_i)

と処理すれば計算できることがわかります。

$\prod\cos\theta_i$ は有理数近似を前計算で求めておくとして、あとは「符号を適切に選べば回転量の総和 $\sum \pm \theta_i$ を（ある程度の範囲内で）好きな値 $\theta$ に収束させられる」「$\times \tan \theta_i$ の計算が楽」という $\theta_i$ の列があればうまく計算できそうです。

実際、$\theta_i=\arctan(2^{-i})$ とすることで $\tan\theta_i = 2^{-i}$ が計算しやすく、 $-\pi/2 \lt \theta \lt \pi/2$ の任意の $\theta$ で計算できるようになります。¹

以上のことから、以下のアルゴリズムで $(x, y)$ を $\theta$ 回転した座標を求めることができます。

前計算として、$\theta_i := \arctan(2^{-i})$ と $K := \prod_{i=0}^N\cos\theta_i = \prod_{i=0}^N 1/\sqrt{1+2^{-2i}} $ の値を求めておく
$0 \le \theta \lt \pi/2$ となるように、先に $\pi/2$ 単位の回転を処理しておく
各 $i$ について、
- $\theta \ge 0$ のとき、$(x, y) \gets (x - y/2^i, y + x/2^i)$、$\theta \gets \theta - \theta_i$ と更新する
- $\theta \lt 0$ のとき、$(x, y) \gets (x + y/2^i, y - x/2^i)$、$\theta \gets \theta + \theta_i$ と更新する
各成分を $K$ 倍したものが答えになる

最初の座標を $(1, 0)$ とすれば三角関数の値 $(\cos \theta, \sin \theta)$ が得られます。

実装(Python)

後で Befunge に移植することを意識した実装になっています。

出力も含めて全て整数演算で構成しました
値も基本的には32bit整数、大きくても64bit整数で扱える範囲に収まるようにしつつ、AC できる精度もちゃんと確保しています
入力に合わせて $\theta_i = \arctan(2^{-i})$ を $\frac{180}{\pi}$ 倍して度数法で記録しています

cordic.py

X = 26

# 分母
# SCALE = 1 << X
# ITER = X
SCALE = 67108864
ITER = 26

# 角度の分母
# DEGSCALE = 1 << 24
DEGSCALE = 16777216

# atan(2^-i) / pi * 180 の分子部分
# import math
# atan_table = [int(round(math.atan(2**-i) / math.pi * 180 * DEGSCALE)) for i in range(ITER)]
atan_table = [754974720, 445687602, 235489088, 119537938, 60000934, 30029717, 15018523, 7509720, 3754917, 1877466, 938734, 469367, 234684, 117342, 58671, 29335, 14668, 7334, 3667, 1833, 917, 458, 229, 115, 57, 29]

# CORDIC スケール補正係数 K = ∏_{i=0}^{N} 1/sqrt(1+2^{-2i})
# K = 1.0
# for i in range(ITER):
#     K /= math.sqrt(1 + 2**(-2*i))
# K_int = int(round(K * SCALE))  # の、分子部分
K_int = 40752055

def cordic_rotate(start_x, start_y, angle_deg):
    """CORDIC で (start_x, start_y) を angle_deg 回転"""
    # 回転量 90° 毎に調整
    mode = angle_deg // 90 % 4
    if mode == 0:
        x, y = start_x, start_y
    elif mode == 1:
        x, y = -start_y, start_x
    elif mode == 2:
        x, y = -start_x, -start_y
    elif mode == 3:
        x, y = start_y, -start_x
    angle_deg %= 90

    # 角度 deg の分子部分
    theta = angle_deg * DEGSCALE

    x *= SCALE
    y *= SCALE

    for i in range(ITER):
        x_shift = x >> i
        y_shift = y >> i
        if theta >= 0:
            # x[i+1] <- x[i] - y[i] / 2^i
            # y[i+1] <- y[i] + x[i] / 2^i
            x, y, theta = x - y_shift, y + x_shift, theta - atan_table[i]
        else:
            # x[i+1] <- x[i] + y[i] / 2^i
            # y[i+1] <- y[i] - x[i] / 2^i
            x, y, theta = x + y_shift, y - x_shift, theta + atan_table[i]
        # x, y は 32bit に収まらないが、64bitには収まるハズ
        assert abs(x) < 2 ** 62 and abs(y) < 2 ** 62

    # スケーリング補正
    assert abs(x * K_int) <= 2 ** 63 and abs(y * K_int) <= 2 ** 63
    x = (x * K_int) // SCALE
    y = (y * K_int) // SCALE
    return x, y

def answer(a):
    # 出力も整数演算のみで
    if a < 0:
        print('-', end='')
        a = -a
    print(a // SCALE, end='.')
    a %= SCALE
    for _ in range(10):
        a *= 10
        print(a // SCALE, end='')
        a %= SCALE

start_x, start_y, angle_deg = map(int, input().split())
x, y = cordic_rotate(start_x, start_y, angle_deg)
answer(x)
print(' ', end='')
answer(y)

実装(Befunge)

コード内に ASCII 制御文字が含まれています。Qiita 上では制御文字を表示できないので、ここでは制御機能用記号に置換しています。正確なコードは AtCoder で提出したものを参照してください。²
きちんと実行できることが確認できたのは現状 AtCoder のジャッジのみです。（ジャッジに使われている Tim's Befunge Compiler をローカルで使っても、おそらく正常に動作すると思います。）
特に、Web上で動作する The BedroomLAN Befunge Interpreter³ や HTML5 Befunge93 Interpreter では正常に動作しませんでした。

"␐@@x`"****:2/"B␑Oxm"**+*+:2/"5␕;~Q"***++:2/"␏␇/OE"***++:2/"␈␛yG"**++:2/"}u␂"**v
v:/2:+1/2:/2:/2:/2:+1/2:/2:/2:/2:+1/2:/2:+1/2:+8/2:+"9"/2:+*"3␉"/2:++*"␜␡m"/2:+<
>2/:2/1+:2/:2/1+ "␚"01>p 01g:#v_ "␚"21p"@@@@␄"****:11p" "*1-01p v
v                     ^0p10:-1< v-1g30\*-10<                    >
>&&&:"Z"%"@@@@"****42p  "Z"/4%  >03p   03g#^_ 11g*:01g/22p01g%32p 11g*:01g/02p0v
vp311p300                                                                p21%g1<
>21g03g`#v_                 v
         >42g0`#v_"+>␈-K␉+!␉"v
                > "->␈+K␉-!␉">ppp 01g02g*12g+:01g22g*32g+13g/ - \01g22g*32g+\13v
^           p31*2g31p30+1g30 p24 -  g0g30g24p21%g10p20/g10:p23%g10p22/g10: + /g<
v                           >"_+Vt_"***+: 01g22g*32g+*11g/ \ 01g02g*12g+*11g/
>:0`!#v_       >:11g/:#v_"0",vv /*52\+"0"%*52:<  v  ,<
,     >"-",01-*^       >01-\  >             :#^_$>:0`|
"                      v,"."$<                       <
                       >11g%25*>:!#v_\25**:11g/"0"+,v
"                              ^-1\             %g11<
^p"␐␆@"                          $$<

各パーツで何をしてるか解説していきます。
以下、コード上の $i$ 行 $j$ 列を「座標 $(i, j)$」と表記します。座標 $(i, j)$ の値を取得・更新するとき、コードとしては j i g や j i p と行・列の順番を逆に指定する点に注意してください。⁴

定数を記録する

"␐@@x`"****:2/"B␑Oxm"**+*+:2/"5␕;~Q"***++:2/"␏␇/OE"***++:2/"␈␛yG"**++:2/"}u␂"**v
v:/2:+1/2:/2:/2:/2:+1/2:/2:/2:/2:+1/2:/2:+1/2:+8/2:+"9"/2:+*"3␉"/2:++*"␜␡m"/2:+<
>2/:2/1+:2/:2/1+ "␚"01>p 01g:#v_ "␚"21p"@@@@␄"****:11p" "*1-01p v
v                     ^0p10:-1<                                 >

atan_table をスタック上に作成する
　まず atan_table[0] = 754974720 を "␐@@x`"****、つまり 96 * 120 * 64 * 64 * 16 で作ります
　このように、"" の中の文字列が数値としてスタックに積まれることを利用して、大きな値を少ない文字数で構成することができます⁵⁶
　atan_table[1] 以降は atan_table[i] = atan_table[i-1] // 2 + n を繰り返す形で構成しています
atan_table[i] を座標 $(0, i)$ に記録する
ITER, SCALE, BASE = 2 ** 31 - 1 をそれぞれ座標 $(1, 2), (1, 1), (1, 0)$ に記録する
　BASE は 32bit の範囲を超える可能性がある値 x を記録するときに、代わりに x // BASE, x % BASE の 2 つに分けて記録するために使います

入力を受け取り、CORDIC ループ前の処理をする

v                               v-1g30\*-10<
>&&&:"Z"%"@@@@"****42p  "Z"/4%  >03p   03g#^_ 11g*:01g/22p01g%32p 11g*:01g/02p0v
                                                                         p21%g1<

入力 $a, b, d$ を受け取り、theta = d % 90 * 2**24 を座標 $(2, 4)$ に記録する
d // 90 回、座標を 90° 回転する
y = b * SCALE を座標 $(2, 2), (2, 3)$ に記録する
　座標 $(2, 2)$ には y // BASEを、座標 $(2, 3)$ には y % BASEを記録しています
　値の取得は 01g22g*32g+、値の更新は :01g/22p01g%32p と処理することになります
x = a * SCALE を座標 $(2, 0), (2, 1)$ に記録する

CORDIC ループ部分の処理を実行する

vp311p300
>21g03g`#v_                 v
         >42g0`#v_"+>␈-K␉+!␉"v
                > "->␈+K␉-!␉">ppp 01g02g*12g+:01g22g*32g+13g/ - \01g22g*32g+\13v
^           p31*2g31p30+1g30 p24 -  g0g30g24p21%g10p20/g10:p23%g10p22/g10: + /g<
v                           >"_+Vt_"***+: 01g22g*32g+*11g/ \ 01g02g*12g+*11g/

i = 0, shift = 1 をそれぞれ座標 $(3, 0), (3, 1)$ に記録する
　shift には 2 ** i が記録されます。Befunge にはビットシフトがないため、 x >> i の代わりに x // shift と計算します
i < ITER なら 3. へ、そうでないなら 6. へ
theta の正負によって、以降の計算の符号を切り替える
　" 文字 列番号 行番号 " p と実行することで、コード上の + - を直接書き換えています。例えば上段 " 直後の3文字 +>␈ は座標 $(8, 62)$ の文字を + に書き換えることを意味しています
$(x, y) \gets (x \mp y/2^i, y \pm x/2^i)$、$\theta \gets \theta \mp \theta_i$ の処理
i += 1, shift *= 2 として 2. へ
y * K_int // SCALE, x * K_int // SCALE をスタックに積む
　この値 / SCALE が出力すべき答えです

10 進小数に変換して出力する

>:0`!#v_       >:11g/:#v_"0",vv /*52\+"0"%*52:<  v  ,<
,     >"-",01-*^       >01-\  >             :#^_$>:0`|
"                      v,"."$<                       <
                       >11g%25*>:!#v_\25**:11g/"0"+,v
"                              ^-1\             %g11<
^p"␐␆@"                          $$<

値が負のときは - を出力し、符号を反転させる
整数部分を取り出し、0 なら 0 を出力して小数部分の出力へ、 0 でないなら -1 をスタックに乗せて整数部分の出力へ
[整数部分の出力] 整数部分を下の位からスタックに乗せていく
　while n != 0: stack.append(n % 10); n //= 10 という処理
-1 が出てくるまでスタックから出力

[小数部分の出力] 筆算の要領で 1 桁ずつ、10 桁分出力する
　Python でいう↓の処理

    a %= SCALE
    for _ in range(10):
        a *= 10
        print(a // SCALE)
        a %= SCALE

1 周目のみ、コードの終点 @ をコード上に設置して、答えの区切りとなるスペースを出力し、1. へ戻る

動作例

いくつかの角度 $\theta$ について、点 $(1, 0)$ を $\theta$ 回転させたときの出力 $(\cos\theta, \sin\theta)$ を見てみます。

$\theta$		$\cos\theta$	$\sin\theta$
$30°$	実数値	0.8660254038	0.5000000000
	出力	0.8660253882	0.5000000000
$45°$	実数値	0.7071067812	0.7071067812
	出力	0.7071067392	0.7071068286
$72°$	実数値	0.3090169944⁷	0.9510565163
	出力	0.3090169280	0.9510565400

6 ～ 7 桁くらいの精度で計算できていることがわかります。

おわりに

esolang には厳しそうな問題が無事に AC できて良かったです。書く前は 25 行という Befunge 特有の制限に収まるか不安でしたが、実際に書いてみると意外と行数に余裕を持って書くことができました。

p 命令を使ってコードを動的に書き換える機能を効果的に使えたのは良かった一方で、コードの後半部分がやや空白の目立つ Befunge らしくないコードになってしまいました。
もっと文字がぎっしり詰め込まれた、Befunge らしいコードを作っていきたいところです。

Befunge はまだ AtCoder のジャッジに正式実装されてませんが、最速で 10 月ごろには正式に実装されるようです。ぜひ A, B 問題を Befunge で解いてみてはいかがでしょうか？⁸

参考文献・関連記事

前回の記事

CORDIC アルゴリズムについて

Befunge で常用対数
今回この問題に挑戦するきっかけとなった記事です。

この $\sum \pm \theta_i$ が $-\pi/2 \le \theta \le \pi/2$ に収束させられる証明は https://mathlog.info/articles/2625 で書かれています。 ↩
AtCoder 上でも水平タブを使ってる 8, 9 行目はレイアウトが崩れています。v ^ 命令で上下に移動するときの移動先がわかりづらくなる…… ↩
The BedroomLAN Befunge Interpreter では、4 行目末尾の > の右で入力をスタックに積み、5 行目先頭の &&& を空白に変更し、8 ～ 10 行目を以下のように書き換えることで実行できました。
```
         >42g0`#v_"+>␈-K␉+␗␉"v
                > "->␈+K␉-␗␉">ppp 01g02g*12g+:01g22g*32g+13g/ - \01g22g*32g+\13v
^ p31*2g31p30+1g30 p24 -  g0g30g24p21%g10p20/g10:p23%g10p22/g10:      +      /g<
```
& 命令で無限ループするのは有名(?)ですが、p 命令の列番号が 70 までしか使えないみたいです。
また、途中で数値が大きくなる可能性があるため、AtCoder で実行した場合よりも誤差が大きくなるかもしれません。 ↩
よく逆にしてバグらせます…… ↩
1 ～ 127 のほぼ全ての値をスタックに積めますが、10(␊改行), 13(␍復帰) は改行されてしまうため、34("ダブルクォーテーション) は文字列モードが終了してしまうためそれぞれ使用不可です。
また、1 ～ 31 と 127 は制御文字なので入力するのが難しいです。今回は Python の pyperclip を利用して pyperclip.copy(chr(16)) などとしてクリップボードから入力しました。 ↩
具体的な構成は、「演算に +* だけを使うもののうち最短のもの」という条件で求めました。$10^9$ 近い値なので全探索は厳しいですが、半分全列挙を使うことで高速に求めることができます。（こんなややこしいことせずとも、127 進数として構成すればそれなりに短いものが得られる……と後から気づきました。） ↩
$\cos(72°) = \frac{\sqrt{5}-1}{4}$ ↩
レートが下がっても責任は負いません。 ↩