はじめに

数学のプロならSageMathなんだろうが、
工学者は他と連携しやすい素のPythonの方が使いやすいだろうということで
sympyでやってみる。
場合によりscipy, numpyも使う。
といってもまずはPython上で数がどのように表されるのかを知らないといけない。

import sympy as sym
sym.init_printing() # Jupyter上でリッチな数式の表示をするためのメソッド

数論

Pythonでの数値演算と数学における実数との差

C言語における__int64やdoubleは64bitである。
すなわち高々$2^{64}$通りの数しか表せない。
(実際もこのオーダーのはずで$2^{63}$程度というビットの無駄遣いの実装はないだろう)
__int64は実数に均等に分布し、範囲内では整数環である。
Pythonの整数は最初から多倍長整数なのでメモリが許す限り数学的な整数と遜色ない。
一方のdouble含め、並の言語における浮動小数は
1の付近では密度が濃いが$2^{64}$の付近ではスカスカになる。(相対的には同じ密度)

x = 1234567890
print(x * x * x * x * x) # 整数はオーバーフローせずに自動拡張される

2867971860299718107233761438093672048294900000

f = 0
for i in range(9):
    f += 0.1
print(f)

0.8999999999999999 # 典型的な浮動小数による丸め誤差

実数である0.1を浮動小数に変換するときに最も近い浮動小数に丸められるが、
それは実数の0.1よりわずかに小さい。
それが9回足すという操作によって顕在化した。
この種の誤差を回避するには符号、整数部、仮数部を持ち、
整数部と仮数部を多倍長整数にした浮動小数クラスを実装し、
手計算するように仮数部を桁合わせしてから整数による演算を行えばよい。
Pythonにはその手のクラスとしてdecimal型があるがsympyを使うのでFloat型を使う。

s = sym.Float(0)
for i in range(9):
    s += sym.Float(0.1)
print(s)

0.900000000000000

s = sym.Float(1.0, 32) # 精度の指定
t = sym.Float(0.1)
s *= 10 ** 30
s += t
print(s)

1000000000000000000000000000000.1 # 桁落ちはギリギリ起こらない

print(sym.Float(0.3, 30))
print(sym.Float("0.300000000000000000000000000000", 30))

0.299999999999999988897769753748　# そもそもPythonの浮動小数を引数で渡すと精度が足りない
0.300000000000000000000000000000 # stringで渡す

有理数

とはいえ浮動小数の思想では精度以上での桁落ちは避けられず、
分配法則も成り立たないので厳密性に欠ける。
そこでRational型を使う。

a = sym.Rational(1,3)
b = sym.Rational(2,5)
c = sym.Rational(4,15)
print(a + b)
print(a + b + c)

11/15
1

Rational型は整数/整数でPythonの整数型は無限精度なので
これは数学的な有理数と同じように扱える。
つまり加減乗除や逆元などが想像通りの体として考えてよいということだ。
ただし大規模な演算をやろうとするとすぐに分母、分子が数万桁になって苦しくなる。

代数的数と特殊な数

冪根

a = 3 ** sym.Rational(1,4)
b = 6 ** sym.Rational(1,4)
a * b

$$\sqrt[4]{2}\sqrt{3}$$

べき乗はPowも使えるが**のオーバーロードが楽。
浮動小数ではなく有理数を渡せば
$x^{\frac{a}{b}}=(\sqrt[b]{x})^a$に従ってある程度根号のまま計算してくれる。

円周率

harf_pi = sym.S.Pi / 2
harf_pi * harf_pi * sym.Rational(3, 4) + harf_pi

$$\frac{\pi}{2}+\frac{3\pi^2}{16}$$

sym.sin(sym.S.Pi * 2 / 3)

$$\frac{\sqrt{3}}{2}$$

$\pi$も有理数係数ならば拡大体のように扱える。
三角関数に入れれば根号で返ってくる。

自然対数の底

e = sym.S.Exp1
e ** 4 + 1

$$1+e^4$$

sym.log(sym.S.Exp1 ** sym.Rational(10, 7))

$$\frac{10}{7}$$

$e$は$\log$に入れると有理数の世界に戻ってくる。
近似じゃなくて厳密なのがよい。

虚数単位i

i = sym.S.ImaginaryUnit
j = sym.sqrt(-1) # こう書いても同じ
print(i*j)
print(sym.S.Exp1 ** (sym.S.ImaginaryUnit * sym.S.Pi)) # オイラーの等式

-1
-1

$e^{i\pi}=-1$を検算。四則演算は略。

実数の壁

どこまでいっても所詮は有理数係数の拡大体で可算無限個の数しか表現できないので、
数学的な実数との間には越えられない壁があることは留意しよう。