数値計算とは、人間が解析的に解くのが難しかったり、解くのに時間がかかりすぎる問題をコンピュータ（計算機）に解かせることである。 コンピュータの発達に伴い、数値計算の手法も大きく発達してきた。 AI でも数値計算の手法が至るところで使われている。 本講義では R を使って数値計算の基本的な手法を学ぶ¹⁾。

• 数の表現と誤差
• 代数方程式
• 数値積分
• 乱数とモンテカルロ法
• 数値微分
• 常微分方程式の数値解法

Moodle（ムードル）とは Web ブラウザでアクセスするオンラインの学習管理システムである。 このクラスでは皆さんが小テストを受験したり、課題を提出するために利用する。

端末や RStudio のアカウントとは異なるので注意。

Moodle へのログインにはユーザ名とパスワードが必要である。 ユーザ名と初期パスワードはメールで配布する。 山本の他の講義で Moodle を使っていた人はそのまま使えるのでメールは送らない。

次のリンクをクリックして Moodle にログインしてみよう。

Moodle Server（非公式）

• パスワードは10文字(以上)に
• パスワードは大文字小文字・数字と記号を交ぜるように
• パスワードは自分は読めるが他人は読めないものに
• 他の場所のパスワードと使い回しをしないように

Moodle の使い方（学生編）

本講義では R を使ったプログラミングを行うため、 R を Web ブラウザ上で使える IDE(統合開発環境) として RStudio を利用する。

• R と RStudio

学内に RStudio サーバを用意してある。ただし、学内ネットワークからのみアクセスできる。

• 学内 RStudio サーバ1（学内限定）
• 学内 RStudio サーバ2（学内限定）

学外から利用したい場合はVPN 接続でアクセスする。

与えられたユーザ名とパスワードで学内 RStudio サーバにログインする。 ユーザ名は大学の端末と同じユーザ名にしてある。

プロジェクト機能とは、特定の作業ディレクトリでスクリプトファイルなどを管理する機能である。

「数値情報処理a」用のプロジェクトを作成するには以下の手順で行う。

1. 右上のプロジェクトボタンをクリックして [New Project] を選ぶ、もしくはメニューの [File]→[New Project] を選ぶと「Create Project」のウィンドウが開く。
2. [New Directory] を選ぶ。
3. [New Project] を選ぶ。
4. 「Directory name:」にプロジェクト名「数値情報処理a」を入力する²⁾。
5. 「Create project as subdirectory of:」 が「~」になっていなければ [Browse] をクリックして [Home] ボタンをクリックした後 [Choose] ボタンをクリックする。
6. [Create Project] をクリックする。

これでプロジェクトが作成できる。

他の講義でこの学内 RStudio サーバを使う場合は、ファイルが混ざってごちゃごちゃにならないように講義ごとにプロジェクトを作成し、プロジェクトを切り替えて使うと良い。

別のプロジェクトへの切り替えは右上のプロジェクトボタンから行う。

R で実行させるコマンドを記述するファイルをスクリプトファイルという。 プログラムでいうソースコードである。

スクリプトファイルの拡張子は .R になる。

左上の [+] ボタンをクリックして [R Script] を選ぶ。 もしくは、メニューの [Files]→[New File]→[R Script] を選ぶ。

左上のペインに 「Untitled1」 という空のスクリプトのタブが作成される。

スクリプトのタブでフロッピーディスクのボタンをクリック（ショートカットキーは Ctrl+S）でファイルとして保存できる。 保存先は RStudio サーバになる。

スクリプトのタブに次のコードを入力して factorial.R というファイル名で保存しよう。

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 を使わなければ以下のコードになる。

factorial2.R

# 階乗の計算
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 進数の状態で保存している。

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 文字で一つの桁を表すことになっている。

例えば 10 進数の整数で $51966$ は、2 進数と 16 進数では以下のように表される。 \[ 51966=(1100\ 1010\ 1111\ 1110)_{2}=(\mathrm{cafe})_{16} \] かっこの下付きの数字は、かっこの中の数字の進数を表す。

\(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}

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 進数の整数 $(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 進数に変換するには $N$ で割った余りを使う。 $N$ で割って商と余りを求め、その商をさらに $N$ で割って商と余りを求めることを繰り返す。 商が $0$ になると終了する。

10進数の $6$ を 2 進数に変換する。

余りを下から上に並べると $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}

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 進数への変換は sprintf() か as.hexmode() を使う。

> sprintf("%x", 51966)
[1] "cafe"
> as.hexmode(51966)
[1] "cafe"

10 進数から 2 進数への変換は intToBits() を使う方法がある（省略）。

整数は符号部（$\pm$：プラスかマイナス）と正の整数で表現される。

\[ 整数=\pm(正の整数) \]

整数は無限に存在するので、有限な容量しか持たないコンピュータでは全ての整数を表現することはできない。 そこで、整数にある決まったデータ長を割り当てて、十分大きな有限の範囲の整数が扱えるようにしてある。

R の “integer”（整数型）では一つの整数に32ビット（4バイト）のデータ長があてがわれてメモリに格納される。 各ビットは以下のように割り当てられる。

• 符号部：1 ビット
• 正の整数：31 ビット

各ビットの格納のイメージは下図である⁵⁾。

ここで $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” になる。

\(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 進数で $1$ より小さい小数を 2 進数に変換するには $2$ をかけた整数部を使う。 $2$ をかけて整数部と小数部に分け、その小数部にさらに $2$ をかけて整数部と小数部に分けることを繰り返す。 小数部が $0.0$ になると終了する。

10 進数の小数 $0.625$ を 2 進数に変換する。

整数部を上から下に並べると $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 \]

10 進数の小数 $0.3$ を 2 進数に変換する。

すぐに分かるが、この例では手順をいくら繰り返しても小数部が $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 進数の数の桁を 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} \]

逆に 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} \]

実数も無限に存在するので、コンピュータでは全ての実数を表現することはできない。 実数でもある決まったデータ長を割り当てて、有限な範囲の実数が扱えるようにしてある。

コンピュータ内部では実数も 2 進数で保存される。 実数は浮動小数点数で実現され、符号部（$\pm$）、指数部、仮数部を使って表現される。 一般的に使われている IEEE 754⁶⁾では以下のように表現される。 \[ 浮動小数点数=\pm(1+仮数部)\times 2^{(指数部-E_{0})} \] ここで $0\le 仮数部 < 1$ である。 $2$ の指数にある $E_{0}$ は定数で、その値は精度に応じて定められている。

浮動小数点数の精度は、使われるデータ長によって単精度、倍精度、四倍精度などといわれる。 四倍長精度以上を多倍長精度ということがある。

R では “numeric”（実数型）に倍精度の浮動小数点数が使われる。 各ビットは以下のように割り当てられる。

• 符号部：1 ビット
• 指数部：11 ビット
• 仮数部：52 ビット

各ビットの格納のイメージは下図である⁷⁾。

ここで $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$ を倍精度の浮動小数点数にしてみよう。

まず、$9.3$ を 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} \] これは無限に続く数なので、データ長さに合わせて丸めてコンピュータで扱えるようにしたのが浮動小数点数である。

次にこれを倍精度の浮動小数点数にする。

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 で直接 $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 桁の精度まで丸めずに計算するために使っている。

bin2dec.R

# 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)」と同じ結果が得られることが分かる。

Moodle Server（非公式）で第1回の課題を行いなさい。

締め切り：2026年4月15日（水）20時

講義を欠席した人は山本（）に連絡すること。