はじめに
こんにちは.みなさんAtCoderやってますか?(私は最近あまり取り組めておりません.)
AtCoderでは三角関数を用いた問題が時折出題されます.
このような問題はC++やPythonを用いている人を桁や精度の関係で苦しめますが,><>
のようなそもそも三角関数が実装されていない言語を用いる人に致命傷を与えます.
そこで,マクローリン展開を用いて$sin$と$cos$を計算するプログラムを実装し,三角関数を用いた問題を解いてみることにしました.
この記事にはソースコードと,実際に問題を解いた結果を示します.ソースコードはまだまだ改善の余地があると感じているので,有識者の方がいらっしゃればアドバイスいただけると幸いです.
sin,cosの実装
$sin$と$cos$の実装は以下の通りです.
v
\ 1[ f 3141592653589793
v +,a < v {+*2}:<
>l3=?!^:{:}$2[>:}$:@{aaa**:@*@*$)?v :}$:@{0v
^ {-*2}:<^v!?(*-1<
>-*$},@$@::}**:}+$r@@v1::$-1$~$~~{] <
^1::+2$*-10$v?=0}-1:{<
v ]~$~~~<
>
v
\ 1[ f 3141592653589793
v +,a < v {+*2}:<
>l3=?!^:{:}$2[>:}$:@{aaa**:@*@*$)?v :}$:@{0v
^ {-*2}:<^v!?(*-1<
>-*$},@$@::}**:}+$r@@v0@:1$-1$~$~~{]<
^1::+2$*-10$v?=0}-1:{<
v ]~$~~~<
>
使用例
$sin(0)$と$cos(0)$の値を計算するプログラムを示します.(なお$sin(0)=0$,$cos(0)=1$です)
0 0 v
v <
\ 1[ f 3141592653589793
v +,a < v {+*2}:<
>l3=?!^:{:}$2[>:}$:@{aaa**:@*@*$)?v :}$:@{0v
^ {-*2}:<^v!?(*-1<
>-*$},@$@::}**:}+$r@@v1::$-1$~$~~{] <
^1::+2$*-10$v?=0}-1:{<
v ]~$~~~<
> v
v $ <
\ 1[ f 3141592653589793
v +,a < v {+*2}:<
>l3=?!^:{:}$2[>:}$:@{aaa**:@*@*$)?v :}$:@{0v
^ {-*2}:<^v!?(*-1<
>-*$},@$@::}**:}+$r@@v0@:1$-1$~$~~{]<
^1::+2$*-10$v?=0}-1:{<
v ]~$~~~<
> $ n " "o n ;
0 1
プログラムは非常にシンプルで0と0をスタックに追加し,sinを計算→swap→cosを計算→再びswap→出力という処理をしています.
使用方法
このプログラムはスタックの一番上の値(一番最後に追加した値)をスタックから取り除き,その値から$sin$,$cos$を計算するものです.
Forth風に書くならば(n -- sin_n)
,(n -- cos_n)
です.
また,以下のような特徴があるためソースコードのどこに書いても動かすことができます.スニペット的に使用できるということです.
- レジスタを使用しない
- スタックの一番上の値以外を変更しない
-
p
命令を使用しない - ジャンプ命令
.
を使用しない
命令ポインタの動き方
左上から始まって,左下へ出ていきます.ループをする箇所は赤色にしています.
先頭一行目は右端から左端へジャンプしますが,これ以外で端から端までのジャンプは行いません.
カスタマイズ方法
2行目に書かれている数字は,マクローリン展開で第何次の項まで計算するかを表す数と,円周率です.上記のプログラムだと,15次の項まで計算し,円周率は小数点以下15桁まで記述されています.(計算精度を考えると少々冗長かもしれません.)
2行目を変更することで,計算回数や精度を変更し,計算負荷を小さくすることが可能です.
例えば以下の通りです.下の例では7次の項まで計算し,円周率を3.14として計算するように設定しています.
v
\ 1[ 7 314
( 以下略 )
sinとcosの違い
ちなみに実は内容はほぼ同じで,6行目が少し違うのみです.(r@@v
の後が1::
か0@:a
の違いです)
$sin$: >-*$},@$@::}**:}+$r@@v1::$-1$~$~~{] <
$cos$: >-*$},@$@::}**:}+$r@@v0@:1$-1$~$~~{]<
実践
折角三角関数を実装したのでAtCoderで試し切りをしてみることにしました.
ABC259B Counterclockwise Rotation
ABC259のB問題は$a,b,d$を受け取ったうえで,
$acos\theta−bsin\theta, asin\theta+bcos\theta$
を計算する問題です.(ただし$\theta=d*\pi/180$)
提出結果
ACしたコードはこちらです.
> 10 > i :0(?v :48*=?v :a=?v :"-"=?v 68*-$a*+ 50.
>~~~ 02. >~*00. >~*00. >~ $ 01- * $ 50.
\ 1[ 3141592653589793
v +,a <
\l2=?!^ * 63a**, ] :
\ 1[ f 3141592653589793
v +,a < v {+*2}:<
>l3=?!^:{:}$2[>:}$:@{aaa**:@*@*$)?v :}$:@{0v
^ {-*2}:<^v!?(*-1<
>-*$},@$@::}**:}+$r@@v0@:1$-1$~$~~{]<
^1::+2$*-10$v?=0}-1:{<
v ]~$~~~<
\ $
\ 1[ f 3141592653589793
v +,a < v {+*2}:<
>l3=?!^:{:}$2[>:}$:@{aaa**:@*@*$)?v :}$:@{0v
^ {-*2}:<^v!?(*-1<
>-*$},@$@::}**:}+$r@@v1::$-1$~$~~{] <
^1::+2$*-10$v?=0}-1:{<
v ]~$~~~<
\ :}$:@{ {{ :}$:@{ }} }@{@
\ 4[ $@ * 01-* @* + n 48*o ]
> {*@*+n ;
コード長: 850 Byte
実行時間: 9 ms
メモリ: 3124 KB
解説
詳細に解説すると長くなってしまうので,簡単に.
1~2行目は入力を受けとる処理です.負の値を入力される可能性があるため少しりょりを工夫しています.
3~5行目は$d$を$\theta$に変換する(度数法→弧度法)処理をしています.またこの後$sin\theta$,$cos\theta$を計算するため値を複製しています.
6~12行目は前述した$cos$のコード.13行目でswap$
して,14~20行目で$sin$のコードで計算しています.
最後に21~23行目で$acos\theta−bsin\theta$, $asin\theta+bcos\theta$をそれぞれ計算して,出力しています.
仕組み
ここから先は興味のある人向けです.
やっていることは非常にシンプルで単に$sin$と$cos$をマクローリン展開したものを計算しているだけです.
私はあまり数学に明るくないですが,マクローリン展開とは
ある関数の値を、x=0を中心としたテイラー展開で表現する手法 (by AI)
らしいです.じゃあテイラー展開はなにかというと,
関数のある一点での導関数の値から計算される項の無限和として関数を表したものである。 (wikipediaより引用)
ということらしいです.詳しい仕組みはともかく,三角関数を単なる四則演算のみで計算できるようになるということです.
$sinx$と$cosx$のマクローリン展開はそれぞれ,
$sinx = x - \dfrac{x^{3}}{3!} + \dfrac{x^{5}}{5!} - \dfrac{x^{7}}{7!} + \dfrac{x^{9}}{9!} ...$
$cosx = 1 - \dfrac{x^{2}}{2!} + \dfrac{x^{4}}{4!} - \dfrac{x^{6}}{6!} + \dfrac{x^{8}}{8!} ...$
のようになっています.これをある程度高次の項まで計算することで,それなりの精度で(少なくともAtCoderでACできる程度に)計算できます.
マクローリン展開の性質上,$x$の絶対値が0に近いほど精度が良くなるので,$x$を$-\pi \leq x \leq \pi$の範囲に丸めてから計算を行います
おわりに
マクローリン展開によって結構精度よく計算することができたのでそれなりに満足しました.
また今回は$sin$,$cos$のみの実装でしたが,実際にAtCoderの問題を確認した結果atan2
も実装しないと解けない問題があると分かりました.
次回はatan2
の実装を記事にしたい...と考えております.