====== 【数値情報処理a】第1回 数の表現と誤差 ======
===== 数値情報処理a について=====
数値計算とは、人間が解析的に解くのが難しかったり、解くのに時間がかかりすぎる問題をコンピュータ（計算機）に解かせることである。
コンピュータの発達に伴い、数値計算の手法も大きく発達してきた。
AI でも数値計算の手法が至るところで使われている。
本講義では R を使って数値計算の基本的な手法を学ぶ((一般的なプログラミング言語（C 言語や Ruby など）では同じ手法が使える。))。
* 数の表現と誤差
* 代数方程式
* 数値積分
* 乱数とモンテカルロ法
* 数値微分
* 常微分方程式の数値解法
===== Moodle =====
Moodle（ムードル）とは Web ブラウザでアクセスするオンラインの学習管理システムである。
このクラスでは皆さんが小テストを受験したり、課題を提出するために利用する。
端末や RStudio のアカウントとは異なるので注意。
----
==== Moodle にログイン ====
Moodle へのログインにはユーザ名とパスワードが必要である。
ユーザ名と初期パスワードはメールで配布する。
山本の他の講義で Moodle を使っていた人はそのまま使えるのでメールは送らない。
次のリンクをクリックして Moodle にログインしてみよう。
[[KMS>|Moodle Server（非公式）]]
* パスワードは10文字(以上)に
* パスワードは大文字小文字・数字と記号を交ぜるように
* パスワードは自分は読めるが他人は読めないものに
* 他の場所のパスワードと使い回しをしないように
{{ :moodle:moodle_student.pdf |Moodle の使い方（学生編）}}
===== 学内 RStudio サーバの利用 =====
本講義では R を使ったプログラミングを行うため、
R を Web ブラウザ上で使える IDE(統合開発環境) として RStudio を利用する。
* [[r:top]]
学内に RStudio サーバを用意してある。ただし、学内ネットワークからのみアクセスできる。
* [[RS>|学内 RStudio サーバ1（学内限定）]]
* [[RS2>|学内 RStudio サーバ2（学内限定）]]
学外から利用したい場合は[[https://www.koeki-prj.org/roy/openvpn/|VPN 接続]]でアクセスする。
----
==== ログイン ====
与えられたユーザ名とパスワードで学内 RStudio サーバにログインする。
ユーザ名は大学の端末のアカウントと同じユーザ名にしてある。
----
==== プロジェクトの作成 ====
プロジェクト機能とは、特定の作業ディレクトリでスクリプトファイルなどを管理する機能である。
「数値情報処理a」用のプロジェクトを作成するには以下の手順で行う。
- 右上のプロジェクトボタンをクリックして ''**[New Project]**'' を選ぶ、もしくはメニューの ''**[File]**''%%→%%''**[New Project]**'' を選ぶと「Create Project」のウィンドウが開く。
- ''**[New Directory]**'' を選ぶ。
- ''**[New Project]**'' を選ぶ。
- 「Directory name:」にプロジェクト名「数値情報処理a」を入力する((RStudio Cloud では日本語は文字化けするのでプロジェクト名は英語にした方がよい。))。
- 「Create project as subdirectory of:」 が「~」になっていなければ ''**[Browse]**'' をクリックして ''**[Home]**'' ボタンをクリックした後 ''**[Choose]**'' ボタンをクリックする。
- ''**[Create Project]**'' をクリックする。
これでプロジェクトが作成できる。
他の講義でこの学内 RStudio サーバを使う場合は、ファイルが混ざってごちゃごちゃにならないように講義ごとにプロジェクトを作成し、プロジェクトを切り替えて使うと良い。
別のプロジェクトへの切り替えは右上のプロジェクトボタンから行う。
===== Rのプログラム例 =====
R で実行させるコマンドを記述するファイルを__**スクリプトファイル**__という。
プログラムでいうソースコードである。
スクリプトファイルの拡張子は .R になる。
----
==== スクリプトの作成 ====
左上の ''**[+]**'' ボタンをクリックして ''**[R Script]**'' を選ぶ。
もしくは、メニューの ''**[Files]**''%%→%%''**[New File]**''%%→%%''**[R Script]**'' を選ぶ。
左上のペインに 「Untitled1」 という空のスクリプトのタブが作成される。
スクリプトのタブでフロッピーディスクのボタンをクリック（ショートカットキーは C-s）でファイルとして保存できる。
保存先は RStudio サーバになる。
スクリプトのタブに次のコードを入力して factorial.R というファイル名で保存しよう。
# 階乗の計算
r <- 1
for (x in 2 : 10) {
r <- r * x
}
print(r)
----
==== スクリプトの実行 ====
左上のペインのスクリプトのタブにある ''**[Source]**'' ボタンをクリックすることで実行できる。
factorial.R を実行するとうまくいけば以下の結果が表示される。
[1] 3628800
これは $10!=1\times 2\times 3\times \cdots \times 9\times 10$ の結果である。
----
==== スクリプトの処理 ====
スクリプトは上から下に書いた順番に処理が行われる（順次処理）。
''for'' というコマンドは、ある条件に従って同じ処理を繰り返す（反復処理、ループ処理）。
factrial.R のスクリプトの具体的な処理は以下の通り。
* 1行目：コメントなので何もしない((''%%#%%'' 以降の文字はコメントとみなされる。))。
* 2行目：変数 ''r'' に 1 を代入する((''%%<-%%'' の代わりに ''%%=%%'' でもよい。))。
* 3行目：変数 ''x'' に 2, 3, 4, ..., 9, 10 を一つずつ代入して「''%%{ }%%''」の中の処理を繰り返す。
* 4行目：変数 ''r'' に ''r * x'' の結果を代入する。
* 6行目：変数 ''r'' の値を出力する。
3〜5行目と同じ処理を ''for'' を使わなければ以下のコードになる。
# 階乗の計算
r <- 1
r <- r * 2
r <- r * 3
r <- r * 4
r <- r * 5
r <- r * 6
r <- r * 7
r <- r * 8
r <- r * 9
r <- r * 10
print(r)
R では $n!$ の計算用に ''factorial(n)'' という関数が用意されているので、上記のような計算をわざわざ行わなくてよい。
> print(factorial(10))
[1] 3628800
===== 2進数と16進数 =====
コンピュータ内部では、データ（数値、文字、画像など）は全て__**2 進数**__の状態で保存している。
2 進数の一桁の情報量を__**ビット（bit）**__という。
さらに、8 ビットの情報量を__**バイト（byte）**__という。
人間にとっては普段使っている__**10 進数**__の方が 2 進数より扱いやすいが、2 進数と 10 進数では桁の対応がシンプルではない。
そこで2 進数四桁（4 ビット）を一桁で表すことができる__**16 進数**__がよく使われる。
10 進数で一つの桁を表すには $0$ から $9$ までの 10 文字を使う。
一方、16 進数で一つの桁を表すには 16 文字必要で、$0$ から $9$ までの 10 文字だけでは足りない。
したがって $\mathrm{a}$ から $\mathrm{f}$ までのアルファベット 6 文字を加えた 16 文字で一つの桁を表すことになっている。
^ 2進数 ^ 10進数 ^ 16進数 ^
| 0 | 0 | 0 |
| 1 | 1 | 1 |
| 10 | 2 | 2 |
| 11 | 3 | 3 |
| 100 | 4 | 4 |
| 101 | 5 | 5 |
| 110 | 6 | 6 |
| 111 | 7 | 7 |
| 1000 | 8 | 8 |
| 1001 | 9 | 9 |
| 1010 | 10 | a |
| 1011 | 11 | b |
| 1100 | 12 | c |
| 1101 | 13 | d |
| 1110 | 14 | e |
| 1111 | 15 | f |
| 10000 | 16 | 10 |
| … | … | … |
例えば 10 進数の整数で $51966$ は、2 進数と 16 進数では以下のように表される。
\[
51966=(1100\ 1010\ 1111\ 1110)_{2}=(\mathrm{cafe})_{16}
\]
**かっこの下付きの数字は、かっこの中の数字の進数を表す。**
{{ nary.png?nolink |2 進数と 16 進数の例}}
===== N進数の整数 =====
==== N進数の整数から10進数への変換 ====
\(m\)桁の N 進数の整数 $(a_{m}a_{m-1}\cdots a_{3}a_{2}a_{1})_{N}$ から 10 進数への変換は以下で計算する。
\begin{eqnarray}
&&\sum_{i=1}^{m}a_{i}\times N^{i-1}\\
&=&a_{m}\times N^{m-1}+a_{m-1}\times N^{m-2}+\cdots+a_{3}\times N^{2}+a_{2}\times N^{1}+a_{1}\times N^{0}
\end{eqnarray}
=== 変換例1 ===
16 進数の整数 $(\mathrm{cafe})_{16}$ を 10 進数に変換する。
まず、各桁を 10 進数に変換すると以下のようになる。
\begin{eqnarray}
a_{4}&=&(\mathrm{c})_{16}=12\\
a_{3}&=&(\mathrm{a})_{16}=10\\
a_{2}&=&(\mathrm{f})_{16}=15\\
a_{1}&=&(\mathrm{e})_{16}=14
\end{eqnarray}
これより
\begin{eqnarray}
&&a_{4}\times 16^{3}+a_{3}\times 16^{2}+a_{2}\times 16^{1}+a_{1}\times 16^{0}\\
&=&12\times 16^{3}+10\times 16^{2}+15\times 16^{1}+14\times 16^{0}\\
&=&12\times 4096+10\times 256+15\times 16+14\\
&=&51966
\end{eqnarray}
となる。
=== 変換例2 ===
2 進数の整数 $(1100\ 1010\ 1111\ 1110)_{2}$ を 10 進数に変換する。
\begin{eqnarray}
&&a_{16}\times 2^{15}+a_{15}\times 2^{14}+a_{14}\times 2^{13}+\cdots+a_{3}\times 2^{2}+a_{2}\times 2^{1}+a_{1}\times 2^{0}\\
&=&1\times 2^{15}+1\times 2^{14}+0\times 2^{13}+\cdots+1\times 2^{2}+1\times 2^{1}+0\times 2^{0}\\
&=&32768+16384+2048+512+128+64+32+16+8+4+2\\
&=&51966
\end{eqnarray}
となる。
----
==== 10進数の整数からN進数への変換 ====
10 進数の整数を N 進数に変換するには $N$ で割った余りを使う。
$N$ で割って商と余りを求め、その商をさらに $N$ で割って商と余りを求めることを繰り返す。
商が $0$ になると終了する。
この原理を簡単に説明するために整数 $x$ が次の形に表せたとしよう。
\[
x = a_{3}\times N^{2}+a_{2}\times N^{1}+a_{1}\times N^{0}
\]
ここで $a_{3},a_{2},a_{1}$ はそれぞれ $0$ から $N-1$ の間の整数である。
$N$ で割って商と余りを求め、その商をさらに $N$ で割って商と余りを求めることを繰り返した結果は以下である。
| ^ 商 | ^ 余り ^
| $(a_{3}\times N^{2}+a_{2}\times N^{1}+a_{1}\times N^{0})\div N=$ | $a_{3}\times N^{1}+a_{2}\times N^{0}$ | $\cdots$ | $a_{1}$ |
| $(a_{3}\times N^{1}+a_{2}\times N^{0})\div N=$ | $a_{3}\times N^{0}$ | $\cdots$ | $a_{2}$ |
| $(a_{3}\times N^{0})\div N=$ | $0$ | $\cdots$ | $a_{3}$ |
最後は商が $0$ になったので終了した。
余りを**下から上に**並べたものが整数 $x$ の N 進数である。
\[
x=(a_{3}a_{2}a_{1})_{N}
\]
\\
=== 変換例 ===
10進数の $6$ を 2 進数に変換する。
| ^ 商 | ^ 余り ^
| $6\div 2=$ | $3$ | $\cdots$ | $0$ |
| $3\div 2=$ | $1$ | $\cdots$ | $1$ |
| $1\div 2=$ | $0$ | $\cdots$ | $1$ |
余りを下から上に並べると $6$ の 2 進数が得られる。
\[
6=(110)_{2}
\]
確認のために、この 2 進数を 10 進数に変換してみると、確かに 6 になることが分かる。
\begin{eqnarray}
&&1\times 2^{2}+1\times 2^{1}+0\times 2^{0}\\
&=&4+2+0\\
&=&6
\end{eqnarray}
===== RにおけるN進数の整数 =====
==== N進数の整数から10進数への変換 ====
N 進数の整数から10進数への変換は上述の式で計算すればよい。
例えば、16 進数の整数 $(\mathrm{cafe})_{16}$ から 10 進数への変換は次のようになる。
\[
12\times 16^{3}+10\times 16^{2}+15\times 16^{1}+14\times 16^{0}
\]
> 12 * 16 ^ 3 + 10 * 16 ^ 2 + 15 * 16 ^ 1 + 14 * 16 ^ 0
[1] 51966
$16^{n}$ は R では ''16 ^ n'' で計算できる。
より簡単な方法として ''strtoi(文字列, N)'' で変換する方法がある。
> strtoi("cafe", 16)
[1] 51966
> strtoi("1100101011111110", 2)
[1] 51966
ここで N 進数の数は文字列として入力する必要がある。
R では「''%%"%%''」（ダブルクォーテーションマーク）で囲ったものは文字列として扱われる。
----
==== 10進数の整数から16進数への変換 ====
10 進数から 16 進数への変換は ''sprintf()'' か ''as.hexmode()'' を使う。
> sprintf("%x", 51966)
[1] "cafe"
> as.hexmode(51966)
[1] "cafe"
10 進数から 2 進数への変換は ''intToBits()'' を使う方法がある（省略）。
===== コンピューターにおける整数の表現 =====
整数は符号部（$\pm$：プラスかマイナス）と正の整数で表現される。
\[
整数=\pm(正の整数)
\]
整数は無限に存在するので、有限な容量しか持たないコンピュータでは全ての整数を表現することはできない。
そこで、整数にある決まったデータ長を割り当てて、十分大きな有限の範囲の整数が扱えるようにしてある。
R の "integer"（整数型）では一つの整数に**32ビット（4バイト）のデータ長**があてがわれてメモリに格納される。
各ビットは以下のように割り当てられる。
* 符号部：1 ビット
* 正の整数：31 ビット
各ビットの格納のイメージは下図である((実際の格納の順番は R が動いているサーバによって違うかもしれない。))。
{{ integer.png?nolink |整数のイメージ}}
ここで $a_{i}$（$i=1,2,\cdots,31$）はそれぞれ 0 か 1 のどちらかである（2進数）。
各ビットから正の整数は以下の式で求められる（2進数から10進数への変換）。
\[
正の整数=\sum_{i=1}^{31}a_{i}\times 2^{i-1}=a_{1}\times 2^{0}+a_{2}\times 2^{1}+\cdots+a_{31}\times 2^{30}
\]
=== 正の整数の最大値 ===
正の整数の最大値は $a_{i}$ がすべて $1$ の場合なので
\[
正の整数の最大値=\sum_{i=1}^{31}1\times 2^{i-1}=\frac{1-2^{31}}{1-2}=2^{31}-1=2147483647
\]
になる。
ここで、等比数列（初項 $1$、公比 $2$）の和の公式を使った。
=== 正の整数の最小値 ===
正の整数の最小値は $a_{i}$ がすべて $0$ の場合なので
\[
正の整数の最小値=\sum_{i=1}^{31}0\times 2^{i-1}=0
\]
になる。
=== 整数の範囲 ===
正の整数の最大値より、R の "integer"（整数型）で表現可能な整数の範囲は
\[
-2147483647\le 32ビット整数\le 2147483647
\]
になる。
最大値が $2147483647$ であることを確かめてみよう。
> class(2147483647L)
[1] "integer"
> class(2147483648L)
[1] "numeric"
警告メッセージ:
非整数値 2147483648L に L がついています; 実数値を使用します
R では数字は一般的に "numeric"（実数型）として扱われるため、
数字の最後に ''L'' を付けて "integer"（整数型）の値であることを明示している。
$2147483647$ を超えると強制的に "numeric" になる。
===== N進数の小数 =====
==== N進数の小数から10進数への変換 ====
\(m\)桁の N 進数の小数 $(0.b_{1}b_{2}b_{3}\cdots b_{m-1}b_{m})_{N}$ から 10 進数への変換は以下で計算する。
\begin{eqnarray}
&&\sum_{i=1}^{m}b_{i}\times N^{-i}\\
&=&b_{1}\times N^{-1}+b_{2}\times N^{-2}+b_{3}\times N^{-3}+\cdots+b_{m-1}\times N^{-m+1}+b_{m}\times N^{-m}
\end{eqnarray}
=== 変換例 ===
2 進数の小数 $(0.101)_{2}$ から 10 進数に変換する。
\begin{eqnarray}
&&1\times 2^{-1}+0\times 2^{-2}+1\times 2^{-3}\\
&=&\frac{1}{2}+\frac{1}{2^{3}}\\
&=&0.5+0.125\\
&=&0.625
\end{eqnarray}
となる。
----
==== 10進数の小数から2進数への変換 ====
10 進数で $1$ より小さい小数を 2 進数に変換するには $2$ をかけた整数部を使う。
$2$ をかけて整数部と小数部に分け、その小数部にさらに $2$ をかけて整数部と小数部に分けることを繰り返す。
小数部が $0.0$ になると終了する。
この原理を簡単に説明するために小数 $y$ （$0\le y < 1$）が次の形に表せたとしよう。
\[
y = b_{1}\times 2^{-1}+b_{2}\times 2^{-2}+b_{3}\times 2^{-3}
\]
ここで、$b_{1},b_{2},b_{3}$ はそれぞれ $0$ か $1$ のどちらかである。
$2$ をかけて整数部と小数部に分け、その小数部にさらに $2$ をかけて整数部と小数部に分けることを繰り返した結果は以下である。
| ^ 整数部 | ^ 小数部 ^
| $(b_{1}\times 2^{-1}+b_{2}\times 2^{-2}+b_{3}\times 2^{-3})\times 2=$ | $b_{1}$ | $+$ | $b_{2}\times 2^{-1}+b_{3}\times 2^{-2}$ |
| $(b_{2}\times 2^{-1}+b_{3}\times 2^{-2})\times 2=$ | $b_{2}$ | $+$ | $b_{3}\times 2^{-1}$ |
| $(b_{3}\times 2^{-1})\times 2=$ | $b_{3}$ | $+$ | $0.0$ |
最後は小数部が $0.0$ になったため終了した。
この整数部を**上から下に**並べたものが小数 $y$ の 2 進数である。
\[
y = (0.b_{1}b_{2}b_{3})_{2}
\]
\\
=== 変換例1 ===
10 進数の小数 $0.625$ を 2 進数に変換する。
| ^ 整数部 | ^ 小数部 ^
| $0.625\times 2=1.25=$ | $1$ | $+$ | $0.25$ |
| $0.25\times 2=0.5=$ | $0$ | $+$ | $0.5$ |
| $0.5\times 2=1=$ | $1$ | $+$ | $0.0$ |
整数部を上から下に並べると $0.625$ の 2 進数が得られる。
\[
0.625=(0.101)_{2}
\]
この 2 進数を 10 進数に変換してみると、確かに $0.625$ になることが分かる。
\[
1\times 2^{-1}+0\times 2^{-2}+1\times 2^{-3}=0.625
\]
=== 変換例2 ===
10 進数の小数 $0.3$ を 2 進数に変換する。
| ^ 整数部 | ^ 小数部 ^
| $0.3\times 2=0.6=$ | $0$ | $+$ | $0.6$ |
| $0.6\times 2=1.2=$ | $1$ | $+$ | $0.2$ |
| $0.2\times 2=0.4=$ | $0$ | $+$ | $0.4$ |
| $0.4\times 2=0.8=$ | $0$ | $+$ | $0.8$ |
| $0.8\times 2=1.6=$ | $1$ | $+$ | $0.6$ |
| $0.6\times 2=1.2=$ | $1$ | $+$ | $0.2$ |
| $0.2\times 2=0.4=$ | $0$ | $+$ | $0.4$ |
| $0.4\times 2=0.8=$ | $0$ | $+$ | $0.8$ |
| $0.8\times 2=1.6=$ | $1$ | $+$ | $0.6$ |
| $\vdots$ | $\vdots$ | $+$ | $\vdots$ |
すぐに分かるが、この例では手順をいくら繰り返しても小数部が $0.0$ になることはないため、終わりがない。
整数部を上から下に並べると $0.3$ の 2 進数が得られる。
途中から $1001$ が無限に繰り返し現れる。
\[
0.3=(0.0\color{blue}{1001}1001\color{blue}{1001}1001\color{blue}{1001}1001\color{blue}{1001}1001\color{blue}{1001}1001\color{blue}{1001}1001\color{blue}{1001}\cdots)_{2}
\]
10 進数では $\frac{1}{3}=0.33333\cdots$ のように無限に同じ数を繰り返す小数（__**循環小数**__）があるので、2 進数でそのような小数があってもおかしくはない。
----
==== N進数の桁の上げ下げ ====
N 進数の数の桁を 1 つ上げるには $N$ をかける。
\[
(a_{4}a_{3}a_{2}a_{1})_{N}\times N=(a_{4}a_{3}a_{2}a_{1}0)_{N}
\]
N 進数の数の桁を 2 つ上げるには $N$ を 2 回かける（$N^{2}$ をかける）。
\[
(a_{4}a_{3}a_{2}a_{1})_{N}\times N^{2}=(a_{4}a_{3}a_{2}a_{1}00)_{N}
\]
\begin{eqnarray}
(a_{4}a_{3}a_{2}a_{1})_{N}\times N
&=&(a_{4}\times N^{3}+a_{3}\times N^{2}+a_{2}\times N^{1}+a_{1}\times N^{0})\times N\\
&=&a_{4}\times N^{4}+a_{3}\times N^{3}+a_{2}\times N^{2}+a_{1}\times N^{1}+0\times N^{0}\\
&=&(a_{4}a_{3}a_{2}a_{1}0)_{N}
\end{eqnarray}
\\
逆に N 進数の数の桁を 1 つ下げるには $N$ で割る。
\[
(a_{4}a_{3}a_{2}a_{1})_{N}\div N=(a_{4}a_{3}a_{2}.a_{1})_{N}
\]
N 進数の数の桁を 2 つ下げるには $N$ で 2 回割る（$N^{2}$ で割る）。
\[
(a_{4}a_{3}a_{2}a_{1})_{N}\div N^{2}=(a_{4}a_{3}.a_{2}a_{1})_{N}
\]
\begin{eqnarray}
(a_{4}a_{3}a_{2}a_{1})_{N}\div N
&=&(a_{4}\times N^{3}+a_{3}\times N^{2}+a_{2}\times N^{1}+a_{1}\times N^{0})\div N\\
&=&a_{4}\times N^{2}+a_{3}\times N^{1}+a_{2}\times N^{0}+a_{1}\times N^{-1}\\
&=&(a_{4}a_{3}a_{2}.a_{1})_{N}
\end{eqnarray}
\\
=== 10進数の例 ===
^ 10進数 ^
| $1234\times 10=12340$ |
| $1234\times 10^{2}=123400$ |
| $1234\div 10=123.4$ |
| $1234\div 10^{2}=12.34$ |
=== 2進数の例 ===
^ 2進数 ^ 10進数との対応 ^
| $(1001)_{2}\times 2=(10010)_{2}$ | $9\times 2=18$ |
| $(1001)_{2}\times 2^{2}=(100100)_{2}$ | $9\times 2^{2}=36$ |
| $(1001)_{2}\div 2=(100.1)_{2}$ | $9\div 2=4.5$ |
| $(1001)_{2}\div 2^{2}=(10.01)_{2}$ | $9\div 2^{2}=2.25$ |
===== コンピュータにおける実数の表現 =====
実数も無限に存在するので、コンピュータでは全ての実数を表現することはできない。
実数でもある決まったデータ長を割り当てて、有限な範囲の実数が扱えるようにしてある。
コンピュータ内部では実数も 2 進数で保存される。
実数は__**浮動小数点数**__で実現され、符号部（$\pm$）、指数部、仮数部を使って表現される。
一般的に使われている IEEE 754((IEEE 754 は浮動小数点数についての国際標準の規格。))では以下のように表現される。
\[
浮動小数点数=\pm(1+仮数部)\times 2^{(指数部-E_{0})}
\]
ここで $0\le 仮数部 < 1$ である。
$2$ の指数にある $E_{0}$ は定数で、その値は精度に応じて定められている。
浮動小数点数の精度は、使われるデータ長によって__**単精度**__、__**倍精度**__、__**四倍精度**__などといわれる。
四倍長精度以上を__**多倍長精度**__ということがある。
^ 精度 ^ データ長 ^ $E_{0}$ ^
| 単精度 | 32 ビット | $127$ |
| 倍精度 | 64 ビット | $1023$ |
| 四倍精度 | 128 ビット | $16383$ |
R では "numeric"（実数型）に**倍精度の浮動小数点数**が使われる。
各ビットは以下のように割り当てられる。
* 符号部：1 ビット
* 指数部：11 ビット
* 仮数部：52 ビット
各ビットの格納のイメージは下図である((実際の格納の順番は R が動いているサーバによって違うかもしれない。))。
{{ double.png?nolink |倍精度浮動小数点数のイメージ}}
ここで $b_{i}$ はそれぞれ $0$ か $1$ のどちらかである（2進数）。
仮数部は 10 進数で以下のように定義される。
\[
仮数部=\sum_{i=1}^{52}b_{i}\times 2^{-i}
\]
=== 仮数部の最大値 ===
仮数部の最大値は $b_{i}$ が全て $1$ の場合なので
\[
仮数部の最大値=\sum_{i=1}^{52}2^{-i}=\sum_{i=1}^{52}2^{-1}\times (2^{-1})^{i-1}=\frac{2^{-1}(1-2^{-52})}{1-2^{-1}}=1-2^{-52}
\]
になる。
ここで、等比数列（初項 $2^{-1}$、公比 $2^{-1}$）の和の公式を使った。
=== 仮数部の最小値 ===
仮数部の最小値は $b_{i}$ が全て $0$ の場合なので
\[
仮数部の最小値=0
\]
になる。
=== 仮数部の範囲 ===
仮数部の最大値と最小値より、仮数部の範囲は
\[
0\le 仮数部\le 1-2^{-52}
\]
になる。
=== 指数部の範囲 ===
指数部は 11 ビットの正の整数なので $0$ から $2047$ の範囲になるが、$0$ には特殊な意味を持たせるので（後述）、
利用できる範囲は
\[
1\le 指数部\le 2047
\]
である。
倍精度では $E_{0}=1023$ が使われるので、$2$ の指数の範囲は
\[
-1022\le 指数部-E_{0}\le 1023
\]
になる。
=== 正の浮動小数点数の範囲 ===
正の浮動小数点数で表現可能な範囲は $0$ を除いて
\[
(1+0)\times 2^{-1022}\le 正の浮動小数点数 \le (2-2^{-52})\times 2^{1023}
\]
になり、計算して四捨五入すると
\[
2.225074\times 10^{-308}\le 正の浮動小数点数 \le 1.797693\times 10^{308}
\]
になる。
ただし、この範囲では実数が $0$ の場合を表現できない。
そこで、仮数部と指数部がともに $0$ の場合を実数の $0$ と定義する。
===== 浮動小数点数の例 =====
例として 10 進数の小数 $9.3$ を倍精度の浮動小数点数にしてみよう。
----
==== 実数の2進数への変換 ====
まず、$9.3$ を 2 進数に変換する。
**① 整数部と小数部に分離**
$9.3$ の整数部は $9$、小数部は $0.3$ である。
\\
**② 整数部を 2 進数へ変換**
整数部 $9$ を 2 進数に変換する。
\[
9=(1001)_{2}
\]
\\
**③ 小数部を2進数へ変換**
小数部 $0.3$ を 2進数に変換する。
\[
0.3=(0.0\color{blue}{1001}1001\color{blue}{1001}1001\color{blue}{1001}1001\color{blue}{1001}1001\color{blue}{1001}1001\color{blue}{1001}1001\color{blue}{1001}\cdots)_{2}
\]
\\
$9.3$ は 2 進数で次のように表せる。
\[
9.3=(1001.0\color{blue}{1001}1001\color{blue}{1001}1001\color{blue}{1001}1001\color{blue}{1001}1001\color{blue}{1001}1001\color{blue}{1001}1001\color{blue}{1001}\cdots)_{2}
\]
これは無限に続く数なので、データ長さに合わせて__**丸めて**__コンピュータで扱えるようにしたのが__**浮動小数点数**__である。
----
==== 浮動小数点数への変換 ====
次にこれを倍精度の浮動小数点数にする。
まず、2 進数の小数の桁を 3 桁ずらして $(1+仮数部)$ になるようにする。これを__**正規化**__という。
\begin{align*}
9.3&=(1001.0\color{blue}{1001}1001\color{blue}{1001}1001\color{blue}{1001}1001\color{blue}{1001}1001\color{blue}{1001}1001\color{blue}{1001}1001\color{blue}{1001}\cdots)_{2}\\
&=(1.0010\color{blue}{1001}1001\color{blue}{1001}1001\color{blue}{1001}1001\color{blue}{1001}1001\color{blue}{1001}1001\color{blue}{1001}1001\color{blue}{1001}\cdots)_{2}\times 2^{3}\\
&=(1+0.0010\color{blue}{1001}1001\color{blue}{1001}1001\color{blue}{1001}1001\color{blue}{1001}1001\color{blue}{1001}1001\color{blue}{1001}1001\color{blue}{1001}\cdots)_{2}\times 2^{3}
\end{align*}
これを
\[
浮動小数点数=\pm(1+仮数部)\times 2^{(指数部-E_{0})}
\]
と比べると
\begin{align*}
指数部&=3+E_{0}=1026\\
仮数部&=(0.0010\color{blue}{1001}1001\color{blue}{1001}1001\color{blue}{1001}1001\color{blue}{1001}1001\color{blue}{1001}1001\color{blue}{1001}1001\color{blue}{1001}\cdots)_{2}
\end{align*}
である。
仮数部は倍精度で52ビットまでの数字しか扱えないので、
仮数部の小数点第53位を__**零捨一入**__((2進数なので四捨五入の代わりに零捨一入にする。))で丸めて小数点第52位までにする。
\begin{align*}
仮数部&=(0.0010\color{blue}{1001}1001\color{blue}{1001}1001\color{blue}{1001}1001\color{blue}{1001}1001\color{blue}{1001}1001\color{blue}{1001}1001\color{blue}{1}\!\!\!\!\!\!\!\!\!\backslash\backslash\backslash\cdots)_{2}\\
&\simeq(0.0010\color{blue}{1001}1001\color{blue}{1001}1001\color{blue}{1001}1001\color{blue}{1001}1001\color{blue}{1001}1001\color{blue}{1001}1010)_{2}
\end{align*}
\\
9.3 は倍精度の浮動小数点数で
\begin{align*}
9.3&\simeq(1 + 0.0010\color{blue}{1001}1001\color{blue}{1001}1001\color{blue}{1001}1001\color{blue}{1001}1001\color{blue}{1001}1001\color{blue}{1001}1010)_{2}\times 2^{3}\\
&\simeq(1.0010\color{blue}{1001}1001\color{blue}{1001}1001\color{blue}{1001}1001\color{blue}{1001}1001\color{blue}{1001}1001\color{blue}{1001}1010)_{2}\times 2^{3}\\
&\simeq(1001.0\color{blue}{1001}1001\color{blue}{1001}1001\color{blue}{1001}1001\color{blue}{1001}1001\color{blue}{1001}1001\color{blue}{1001}1010)_{2}
\end{align*}
となる。
===== Rにおける浮動小数点数 =====
R で直接 $9.3$ を入力してみよう。
> 9.3
[1] 9.3
ちゃんと $9.3$ が表現されているように見えるが、これは R ではデフォルトで 7 桁までの表示になっているためである。
オプション ''digits'' を使って表示桁数を最大の 22 桁に変更してみよう。
> print(9.3, digits = 22)
[1] 9.300000000000000710543
これより、浮動小数点数では $9.3$ を正確に表現できていないことが分かる。
$9.3$ をこちらが入力しても、コンピュータ内部では上述したように浮動小数点数による__**丸め**__が行われているので、$9.3$ ぴったりにはならない。
この丸めによって生じる誤差を__**丸め誤差**__という。
先ほど求めた $9.3$ の倍精度の浮動小数点数
\[
(1001.0\color{blue}{1001}1001\color{blue}{1001}1001\color{blue}{1001}1001\color{blue}{1001}1001\color{blue}{1001}1001\color{blue}{1001}1010)_{2}
\]
を 10 進数に変換して丸め誤差を確かめよう。
次のプログラム bin2dec.R は 2 進数の実数を 10 進数に変換するプログラムである。
''x'' に 2 進数の浮動小数点数を**文字列**で入力する必要がある。
Rmpfr は任意精度用のパッケージで、70 桁の精度まで丸めずに計算するために使っている。
# 2進数から10進数への変換
# 任意精度用のパッケージ
library(Rmpfr)
# 精度（ビット）
pbit <- 70
# 2進数（文字列）
x <- "1001.0100110011001100110011001100110011001100110011010"
# xを整数部と小数部に分ける
intdec <- unlist(strsplit(x, ".", fixed = TRUE))
# 整数部の桁数
idigit <- nchar(intdec[1])
# 整数部と小数部を結合して1文字ずつに分割
a <- unlist(strsplit(paste0(intdec), ""))
# 1桁ずつ計算して足し合わせ
r <- 0
for (i in 1 : length(a)) {
r <- mpfr(r + strtoi(a[i], 2) * 2 ^ (idigit - i), pbit)
}
# 10進数の出力
print(r, digits = 22)
結果
[1] 9.300000000000000710543
「''print(9.3, digits = 22)''」と同じ結果が得られることが分かる。
=== 正の浮動小数点数の最大値 ===
R で扱うことができる正の浮動小数点数の最大値、つまり R でどこまで大きい数を表現できるかは ''.Machine$double.xmax'' で確認できる。
> .Machine$double.xmax
[1] 1.797693e+308
''1.797693e+308'' は $1.797693\times 10^{308}$ のことである。
''1.797693e+308'' を超えると ''Inf''（無限大）として扱われる。
> 1.797693e+308
[1] 1.797693e+308
> 1.797694e+308
[1] Inf
\\
=== 正の浮動小数点数の最小値 ===
R で扱うことができる正の浮動小数点数の（$0$ より大きい）最小値、つまり R でどこまで数を小さく分けられるかは ''.Machine$double.xmin'' で確認できる。
> .Machine$double.xmin
[1] 2.225074e-308
''2.225074e-308'' は $2.225074\times 10^{-308}$ のことである。
===== 課題 =====
[[KMS>|Moodle Server（非公式）]]で第1回の課題を行いなさい。
締め切り：2026年4月15日（水）20時
講義を欠席した人は山本（{{:email.png?nolink|email}}）に連絡すること。