多倍精度浮動小数点数や10進小数や乱数や準数値演算に関していろいろ自作した。
たしかTAOCPでDonald Knuth先生が「この分野はつまらないと思われがちだけど、実際に学んでみると面白い」みたいなことを言ってたけど、ほんとその通りだった。

実装

C言語→Python→Schemeの順に低レイヤから実装。

C言語
まずC言語で10進小数とその基本関数（四則演算や指数対数や一様乱数）を実装（srcディレクトリ）。
10進小数なのでメルセンヌツイスタとかは使えず、線形合同法で一様乱数を生成。

Python
上記のC言語での実装をPythonから使えるように、CFFI実装（bigfloat.py）。
既存のScheme処理系を改造して、10進小数を扱うSchemeを実装した（softscheme.py）。
ついでに、SchemeからPythonの関数を呼んだり、実行中にPythonモードとSchemeモードを切替えれるようにするなど改造。

Scheme
一様乱数を素にして、正規分布、ガンマ分布、T分布に従う乱数を実装した（libディレクトリ）。
ノンパラメトリックな検定が欲しかったのでKS検定(Kolmogorov-Smirnov検定)も実装した。

ソースコードと開発環境
ソースコードはここ。

開発環境

   CPU: Intel(R) Core(TM) i7-4770HQ CPU @ 2.20GHz 64 bits
    OS: Ubuntu 14.04 LTS
   gcc: 4.8.4 (Ubuntu 4.8.4-2ubuntu1~14.04.3)
python: 3.6.0
  cffi: 1.9.1

実用目的の実装ではないです。当時手を動かしながら学ぶのが目的だったので変な実装になっているけど、そのうち完成度高いやつを作りたいな。
以下にデモをいくつか紹介。

デモ①有効桁

任意桁の精度で計算できる（ここでは50桁に設定してある）。

demo/signigicant_digits.scm

(define f 
    (lambda (a b)
        (if (= a (+ a b))
            0
            (+ 1 (f a (/ b 10))))))
(print (int (f 1 1)))

bash

$python softscheme.py demo/significant_digits.scm
50

比較対象としてpythonで同じことをすると、普通にやると16桁になる。

demo/signigicant_digits.py

a=1
b=1
counter = 0
while True:
	counter = counter + 1
	b = b/10
	if a == a+b:
		break
print(counter)

bash

$python demo/significant_digits.py
16

デモ②丸め誤差

テント写像を使った数列

\begin{eqnarray}
x_{n+1}=\left\{ \begin{array}{l}
\mu x_n & (x_n < 1/2) \\
\mu (1-x_n) & (1/2 \leq x_n) \\
\end{array} \right.
\end{eqnarray}

について、μ=2として初期値0.1とすると、0.1→0.2→0.4→0.8→0.4→0.8→…となって「0.4→0.8→」の繰り返しになる。

demo/tent_map.scm

(define tent_func
    (lambda (x)
        (if (< x 0.5)
            (* x 2)
            (- 2 (* x 2)))))
(define extend
    (lambda (L)
        (append L (list (tent_func (nth -1 L))))))
(define f
    (lambda (L n)
        (if (> (length L) n)
            L
            (f (extend L) n))))
(plot (f (list 0.1) 100))

上記だけだと当たり前だけど、下記のように普通に（2進小数で）これをやろうとすると、想定外のカオス的振る舞いが発生する。
原因は、丸め誤差と関数の初期値敏感性が組み合わさったため。
参考文献は「カオスとフラクタル」（山口昌哉）など。

demo/tent_map.py

import matplotlib.pyplot as plt
def tent(x):
	if 2*x < 1:
		return 2*x
	else:
		return 2-2*x
v = 0.1
data = []
for i in range(100):
	data.append(v)
	v = tent(v)
plt.plot(data)
plt.show()

デモ③素数

有限群のラグランジュの定理によると、下記が成り立つ。

pを2でも5でもない素数とすると、分母がpの分数について、
p-1 = （循環節の長さ） * （循環節の種類）.

参考文献は「素数はめぐる」（西来路文朗、清水健一）など。
数論好きな人向けの循環小数を眺めるデモ。確かに成り立ってる。

bash

$python softscheme.py demo/prime_reciprocal.scm

prime:3
length_of_repetend:1
num_of_subsets:2
3-1=1*2
+0.33333333333333333333333333333333333333333333333333
+0.66666666666666666666666666666666666666666666666667
prime:7
length_of_repetend:6
num_of_subsets:1
7-1=6*1
+0.14285714285714285714285714285714285714285714285714
+0.42857142857142857142857142857142857142857142857143
+0.28571428571428571428571428571428571428571428571429
+0.85714285714285714285714285714285714285714285714286
+0.57142857142857142857142857142857142857142857142857
+0.71428571428571428571428571428571428571428571428571
prime:11
length_of_repetend:2
num_of_subsets:5
11-1=2*5
+0.090909090909090909090909090909090909090909090909091
+0.90909090909090909090909090909090909090909090909091
+0.18181818181818181818181818181818181818181818181818
+0.81818181818181818181818181818181818181818181818182
+0.27272727272727272727272727272727272727272727272727
+0.72727272727272727272727272727272727272727272727273
+0.36363636363636363636363636363636363636363636363636
+0.63636363636363636363636363636363636363636363636364
+0.45454545454545454545454545454545454545454545454545
+0.54545454545454545454545454545454545454545454545455
prime:13
length_of_repetend:6
num_of_subsets:2
13-1=6*2
+0.076923076923076923076923076923076923076923076923077
+0.76923076923076923076923076923076923076923076923077
+0.69230769230769230769230769230769230769230769230769
+0.92307692307692307692307692307692307692307692307692
+0.23076923076923076923076923076923076923076923076923
+0.30769230769230769230769230769230769230769230769231
+0.15384615384615384615384615384615384615384615384615
+0.53846153846153846153846153846153846153846153846154
+0.38461538461538461538461538461538461538461538461538
+0.84615384615384615384615384615384615384615384615385
+0.46153846153846153846153846153846153846153846153846 
+0.61538461538461538461538461538461538461538461538462
prime:17
length_of_repetend:16
num_of_subsets:1
17-1=16*1
+0.058823529411764705882352941176470588235294117647059
+0.58823529411764705882352941176470588235294117647059
+0.88235294117647058823529411764705882352941176470588
+0.82352941176470588235294117647058823529411764705882
+0.23529411764705882352941176470588235294117647058824
+0.35294117647058823529411764705882352941176470588235
+0.52941176470588235294117647058823529411764705882353
+0.29411764705882352941176470588235294117647058823529
+0.94117647058823529411764705882352941176470588235294
+0.41176470588235294117647058823529411764705882352941
+0.11764705882352941176470588235294117647058823529412
+0.17647058823529411764705882352941176470588235294118
+0.76470588235294117647058823529411764705882352941176
+0.64705882352941176470588235294117647058823529411765
+0.47058823529411764705882352941176470588235294117647
+0.70588235294117647058823529411764705882352941176471

デモ④乱数

一様乱数は線形合同数列を用いて実装した。


X_{n+1} = a X_n + c 　(mod 　m)

ここで有効桁をeとして、m=10^e．
c = 0 の場合は最長周期を得られないが、定数を上手く選べば十分長い周期になる。
今回の実装では、a = 117 or 171 (mod 200)のとき、周期は5*10^(e-2)となることを利用した。
参考文献は「The Art of Computer Programming Volume 2」(Donald Knuth)など。

上記の一様乱数を用いて、正規分布乱数、ガンマ分布乱数、T分布乱数を生成した。
参考文献は「確率分布乱数生成法」(四辻哲章)など。

lib/rand_norm.scm

(define rand_norm
    (lambda ()
        (begin 
            (define rand_s
                (lambda ()
                     (- (* 2 (rand)) 1)))
            (define u (rand_s))
            (define v (rand_s))
            (define w 
                (+ (* u u) (* v v)))
            (if (< w 1)
                (* u (sqrt (/ (* -2 (ln w)) w)))
                (rand_norm)))))

lib/rand_gamma.scm

(define rand_gamma
	(lambda (alpha beta)
		(begin
			(define d (- alpha (/ 1 3)))
			(define c (/ 1 (sqrt (* 9 d))))
			(define gen_z 
				(lambda ()
					(begin 
						(define z (rand_norm))
						(define v (+ 1 (* c z)))
						(if (> v 0)
							z
							(gen_z)))))
			(define z (gen_z))
			(define v (+ 1 (* c z)))
			(define w (* v (* v v)))
			(define y (* d w))
			(define u (rand))
			(if (+ 	(<  u 
					  (- 1 (* 0.0331 (* (* z z) (* z z)))))
					(<= (ln u)
					  (+ (/ (* z z) 2) (+ (* d (ln w)) (- d y)))))
				(* beta y)
				(rand_gamma alpha beta)))))

lib/rand_student.scm

(define rand_student
	(lambda (r)
		(begin
			(define rdiv2 (/ r 2))
			(define d     (sqrt rdiv2))
			(/ (* d (rand_norm))
			   (sqrt (rand_gamma rdiv2 1))))))

Scipy.statsの乱数と比較。各サンプル数2000．

正規分布：

ガンマ分布（形状母数shape parameterは2、尺度母数scale parameterは1）：

T分布（自由度3）：

感想

ここでは紹介しなかったけど、テイラー展開とかニュートン法とかも役に立ってたんすね。
計算機の数学体系は、純数学的な論理体系とは異なっていて独特だと感じた。
てか人類の英知すごいわ。

参考資料

Donald Knuth 「The Art of Computer Programming Volume 1」
Donald Knuth 「The Art of Computer Programming Volume 2」
Peter Norvig 「(An ((Even Better) Lisp) Interpreter (in Python))」
William H. Press 「Numerical Recipes in C」
四辻哲章「確率分布乱数生成法」
山口昌哉「カオスとフラクタル」
山口昌哉「数値解析と非線形現象」
大石進一「精度保証付き数値計算」
西来路文朗、清水健一「素数はめぐる」