昔作った関数電卓の供養。実装した計算が正しいことも効率的なことも保証できないが、せっかく調べて考えた数式が詰まっているので残しておく。
使用した言語がかなり制限があったため苦労した。
-
計算できるのは加減乗除、剰余、切り落とし(int)のみ
- 三角関数、対数関数などの数学関数は無し
- ビット演算も無し
- 関数はまともに定義できない
- ジャンプして戻るだけ ＝ 仮引数もスコープも無し
- 配列を使うのも面倒
- evalで「名前＋数字」の変数を作ることで実現
数学関数
ほとんどの場合はテイラー展開に持ち込む。その際に元の数式通りではなく、級数の収束を速めるよう引数の範囲を制限している。
平方根
冪演算の特別な場合であり √a = a1/2 だが、よく使うので別に実装する。
これはニュートン法で求める。2次収束（1ステップで精度が倍増）なので基本的に速いが、 √0 のときは遅いので直接答えを返すようにする。
u_{n+1} = \frac{1}{2} \left( u_n + \frac{a}{u_n} \right)
対数関数
底の変換公式により、ひとつの底について計算できれば全ての底に対応できる。なので最も簡単な自然対数を実装する。
\log_a{x} = \frac{\log_e{x}}{\log_e{a}}
自然対数はテイラー展開を利用して求める。ただし真数 x が 1 から離れていると収束しない/収束が遅いので、真数を 1 に近づける処理ができるよう式変形しておく。
\log_e{x} = \log_e{(p^m x')} = m \log_e{p} - \sum_{n=1}^{\infty} \frac{(1-x')^n}{n} \\
\left( p > 1, \;\; m \in \mathbb{Z}, \;\; x' = \frac{x}{p^m}, \;\; |1-x'| \le \frac{p-1}{p+1} \right)
|1-x'| をなるべく小さい範囲に収めることと、コンピューター上での計算のしやすさから、 p = 2 を選んだ。
指数関数・冪関数
指数が整数のときはバイナリ法が使える。
a^m =
\begin{cases}
1 / a^{-m} & (m<0) \\
1 & (m=0) \\
\left( a^{\lfloor m/2 \rfloor} \right) ^2 & (m \text{ is even}) \\
\left( a^{\lfloor m/2 \rfloor} \right) ^2 a & (m \text{ is odd}) \\
\end{cases}
指数が小数のときは、底を自然対数の底 e に変換することでテイラー展開が使いやすくなる。ただし、指数が 0 から離れていると収束が遅いので、指数を整数部と小数部に分けて、小数部にのみテイラー展開を使う。
a^x = e^{x \log_e{a}} = e^{m+x'} = e^m \cdot \sum_{n=0}^{\infty} \frac{{x'}^n}{n!} \\
\left( m \in \mathbb{Z}, \;\; |x'| \le 0.5 \right)
なお、テイラー展開した級数の各項が前の項に何かを掛けた形になっていれば、以下のようにして高次の項から足していき計算誤差を減らせる。
\begin{align}
\sum_{n=0}^{\infty} \frac{{x}^n}{n!} &= 1 + \frac{x}{1!} + \frac{x^2}{2!} + \frac{x^3}{3!} + \cdots \\
&= 1 + \frac{x}{1} \left( 1 + \frac{x}{2} \left( 1 + \frac{x}{3} \left( 1 + \cdots \right) \right) \right) \\
\end{align}
三角関数
三角関数は周期的なので、1周分の角度（-π≦x≦π）に対して計算できればいい。さらに、sinとcosは90°毎に入れ替わるため、 -π/4 ≦ x ≦ π/4 だけ求められればいい。（さらに非負のみにもできるが、テイラー展開を見越すと不要と考えた）
\left( \cos{(x+m\pi/2)}, \sin{(x+m\pi/2)} \right) =
\begin{cases}
\left( +\cos{x}, +\sin{x} \right) & (m \equiv 0 \pmod 4) \\
\left( -\sin{x}, +\cos{x} \right) & (m \equiv 1) \\
\left( -\cos{x}, -\sin{x} \right) & (m \equiv 2) \\
\left( +\sin{x}, -\cos{x} \right) & (m \equiv 3) \\
\end{cases}
これで |x| ≦ π/4 < 0.8 と十分小さくなったので、テイラー展開ですぐ求められる。
\begin{align}
\cos{x} &= \sum_{n=0}^{\infty} (-1)^n \frac{x^{2n}}{(2n)!} \\
\sin{x} &= \sum_{n=0}^{\infty} (-1)^n \frac{x^{2n+1}}{(2n+1)!} \\
\end{align}
tan は sin から cos を割ればいい。
逆三角関数
こちらは atan2(y,x) を軸に実装する。座標 (x,y) を引数に与えるとx軸との角度が求まり（複素数の偏角と同じ）、他の逆三角関数はこれひとつで模倣できる。
\begin{align}
\operatorname{asin}{x} &= \operatorname{atan2}{\left( x, \sqrt{1-x^2} \right)} \\
\operatorname{acos}{x} &= \operatorname{atan2}{\left( \sqrt{1-x^2}, x \right)} \\
\operatorname{atan}{x} &= \operatorname{atan2}{\left( x, 1 \right)} \\
\end{align}
atan2(y,x) の計算は、三角関数と逆の要領で 0 ≦ |y| ≦ x の領域（x軸正方向±45°）だけ求められればいい。
\operatorname{atan2}{(y,x)} =
\begin{cases}
\operatorname{atan2}{(-y,-x)} - 2 \cdot \pi/2 & (x \le y \lt 0) \\
\operatorname{atan2}{(+x,-y)} - 1 \cdot \pi/2 & (y \lt -|x|) \\
\operatorname{atan}{(y/x)} & (|y| \le x) \\
\operatorname{atan2}{(-x,+y)} + 1 \cdot \pi/2 & (|x| \lt y) \\
\operatorname{atan2}{(-y,-x)} + 2 \cdot \pi/2 & (0 \lt y \le -x) \\
\end{cases}
これによって atan の無限級数が使える。普通のテイラー展開では |y/x| が1に近いと遅いので、オイラーの方法を使う。
\begin{align}
\operatorname{atan}{z}
&= \frac{z}{1+z^2} \sum_{n=0}^{\infty} \prod_{k=1}^{n} \frac{2kz^2}{(2k+1)(1+z^2)} \\
&= \frac{z}{1+z^2} \left( 1 + \frac{2}{3} \frac{z^2}{1+z^2} \left( 1 + \frac{4}{5} \frac{z^2}{1+z^2} \left( 1 + \cdots \right) \right) \right)
\end{align}
双曲線関数
指数関数で求める。
\begin{align}
\sinh{x} &= \frac{e^{x} - e^{-x}}{2} \\
\cosh{x} &= \frac{e^{x} + e^{-x}}{2} \\
\tanh{x} &= \frac{e^{x} - e^{-x}}{e^{x} + e^{-x}} \\
\end{align}
逆双曲線関数
対数関数で求める。
\begin{align}
\operatorname{asinh}{x} &= \log_e{\left( x + \sqrt{x^2+1} \right)} \\
\operatorname{acosh}{x} &= \log_e{\left( x + \sqrt{x^2-1} \right)} \\
\operatorname{atanh}{x} &= \frac{1}{2} \log_e{\left( \frac{1+x}{1-x} \right)} \\
\end{align}
上の asinh の計算は x が負の方向に大きいと桁落ちするので、 asinh(x) = -asinh(-x) と符号を分離して計算したほうが良さそう。 acosh も x が1に近いと桁落ちが危ない。
階乗・ガンマ関数
階乗 n! は、 1 ～ n の整数を全て掛けた数。この計算を n が小数（や複素数）でもいいように拡張したものがガンマ関数 $\Gamma(z)$ 。関数電卓として小数に対応したかったが、他と比べて段違いに難しかった…。
普通の階乗計算と同じように、 n = 0 の近くに達するまでは再帰的に縮小することができる（負の整数はゼロ除算が発生するため不可）。 n が整数であれば基本ケース 0! = 1 によって簡単に終われる。
n! = \Gamma(n+1) =
\begin{cases}
n \Gamma(n) & (n \gt \frac{1}{2}) \\
\frac{1}{n+1} \Gamma(n+2) & (n \le -\frac{1}{2}) \\
\end{cases}
あとは基本ケース |n| ≦ 1/2 を計算できればいい。ガンマ関数の自然対数をテイラー展開した式が比較的シンプルだったので使うことにした。 γ をオイラーの定数、 ζ(k) をゼータ関数として、
\log_e{\Gamma(1+z)} = -\gamma z + \sum_{k=2}^{\infty} \frac{\zeta(k)}{k} (-z)^k
ゼータ関数も無限和なので計算する必要がある。以下の式で s が小さいときは収束が遅いので、 s ≦ 4 は値をハードコーディングすることにした。
\zeta(s) = \sum_{n=1}^{\infty} \frac{1}{n^s}
いま見返すと、ゼータ関数から1を引いておけば各項の係数が急速に0に近づくため収束を改善できそう。1引いた分を寄せ集めると自然対数になり（下式左辺）、後で指数関数を作用させるので直接計算する必要が無い。
\log_e{\Gamma(1+z)} + \log_e{(1+z)} = (\gamma-1) (-z) + \sum_{k=2}^{\infty} \frac{\zeta(k)-1}{k} (-z)^k
おまけ
関数電卓には学校で使いそうな二項演算をいくつか入れていた。四則演算での実装は難しくない。
順列・組合せ
組合せの計算で計算途中の値が大きくなり過ぎないように少し気を使ったが、あまり工夫はしていない。
- 順列 ${}_{n}P_{r} = n \cdot (n-1) \cdot \dots \cdot (n-r+1)$
- 組合せ ${}_{n}C_{r} = {}_{n}P_{r} / {}_{r}P_{r}$
- r を小さくする： ${}_{n}C_{r} = {}_{n}C_{n-r}$
- 割り算を織り交ぜる： ${}_{n}P_{r} / {}_{r}P_{r} = \frac{n-r+1}{1} \cdot \frac{n-r+2}{2} \cdot \dots \cdot \frac{n}{r}$
- 重複順列 ${}_{n}\Pi_{r} = n^r$
- 重複組合せ ${}_{n}H_{r} = {}_{n+r-1}C_{r}$
最大公約数・最小公倍数
最大公約数はユークリッドの互除法で計算する。それが求まれば最小公倍数も簡単に計算できる。
\gcd{(a,b)} =
\begin{cases}
|a| & (b=0) \\
\gcd{(b,a \bmod b)} & (\text{otherwise}) \\
\end{cases}
\\
\operatorname{lcm}{(a,b)} =
\begin{cases}
0 & (a=b=0) \\
\frac{|ab|}{\gcd{(a,b)}} & (\text{otherwise}) \\
\end{cases}
直角三角形の斜辺
図形の計算でよく出てくるし、ベクトルや複素数の大きさを求めるのにも使う。プログラミング言語によっては hypot という名前で入っている。
平方根は実装済みなので三平方の定理で普通に計算する。厳密に実装するなら2乗したときにオーバーフロー/アンダーフローしないよう対策が必要だが、面倒なのでサボった。
c^2 = a^2 + b^2 \\
c = \sqrt{a^2 + b^2}
逆算用の関数も欲しかったので入れておいた。
a = \sqrt{\left| c^2 - b^2 \right|}
逆数和の逆数
正式な名前が分からない演算。仕事算、並列抵抗や直列コンデンサの合成、レンズの公式、二体問題での換算質量、調和平均などに使う。
\frac{1}{f} = \frac{1}{a} + \frac{1}{b} \\
f = \frac{ab}{a+b}