情報地球科学演習 (2024年前学期) 第09回

余力がある人向きの課題: LibreOffice を使ってみよう

これまでの演習で扱ってきたような統計的な処理は当然ながら、Excel のような表計算ソフトでやることもできる (というか、普通はそうする)。そこで、余力のある人は LibreOffice というソフトウェアを使って、統計的な処理ができることを試してみよう。なおこんなことをすれば、わざわざ Fortran 90 でプログラミング演習した (させられた) 苦労が台無しになってしまいそうな気もするけれど、複数の異なる手段で同じ問題を解かせてみる ことは検算する上で非常に役に立つものだし、またそういう手段をいくつも手元に用意しておくことは戦略として重要である、と言い張って教員は自分を正当化することにする。

LibreOffice とは、いわゆる「オフィススイート」と呼ばれるソフトウェアの1つであり、

Microsoft Office とか市販の主だったオフィススイートと (ほぼ) 同様に使える
フリーソフトウェアとして自由にダウンロードし、使用・再配布することができる
Windows はもちろん Linux や Mac でも利用できる

といった (実にありがたい) 特徴をもっている。ということで (自他ともに認める貧乏症の) 亀山は、ン万円も出してまで私物のパソコン用に Microsoft Office を買おうとなんて思わない (し、本当に買わない)。また実際、普段の作業のほぼ全てを Linux 上で行っている亀山は、(事務的仕事を除いた) 自分の作業に Microsoft Office は全く使わないし、必要な場合でも代わりに Linux 上で LibreOffice を使えば十分に用は足りている。

この演習で使っている Linux にも (バージョンがかなり古いけど) LibreOffice がインストールされている。例えば端末の中で

$ libreoffice &

と実行してみると、LibreOffice の起動画面が立ち上がる (画面の上端にある「メニューバー」の「アプリケーション」->「オフィス」からでもOK)。起動画面からも分かるように、

Word の代わりになりそうな Writer
Excel の代わりになりそうな Calc
PowerPoint の代わりになりそうな Impress

という機能が LibreOffice には備わっている。実際、起動画面の左側の「新規作成」のところにある「Calc 表計算ドキュメント」をクリックすると、確かに表計算ソフトらしき画面が出てくることが確認できるだろう。

いま起動した LibreOffice の Calc の中に、これまでの演習で使ってきた matsuyama-kion.txt の中身を読み込ませてみよう。ちょっとややこしいが、次のような手順で読み込ませる (インポートする) ことができるはずだ。ちなみに、これとほとんど同じ手順をたどれば Excel に読み込ませることもできる (ことを2019年度の受講学生に教えてもらった)。

Calc のメニューバーから「ファイル」->「開く」として出てきた画面の中から、matsuyama-kion.txt を選択する。 matsuyama-kion.txt が候補に挙がってこない人は、「すべてのファイル」を候補に挙げる設定になっているかどうか確認せよ。
出てきた「テキストのインポート」のうち、「区切りのオプション」を「固定幅」と設定し、「フィールド」の欄でマウス操作によってデータの区切りの位置を指定する。
全ての区切りの位置が指定できたら「OK」をクリック。

また LibreOffice の Calc の形式で保存しておきたい人は、Calc のメニューバーから「ファイル」->「名前をつけて保存」を実行すればよい。その際、ファイルの種類は「ODF 表計算ドキュメント (.ods)」を選択するのが吉。ちなにみ ODF とは OpenDocument Format の頭文字をとったもので、ODF 形式のファイルは LibreOffice だけでなく (Microsoft Office を含む) 他のオフィススイートでも読み書きできるようになっている。

LibreOffice の Calc には、自分で式を入力して「表計算」することもできる。またもちろん average とか stdev をはじめ、いくつもの有用なマクロが標準で用意されている。各自でいろいろ試してみてほしい。

離散数学を体感してみよう

1限の授業でも述べた通り、計算機の中で登場する数には「無限」とか「連続」という (「数直線」から素直に想像される) イメージとは異なる性質を持っている。以下では、そのような例をいくつか体感してみよう。

以下に挙げた例の中から少なくとも1つ以上のプログラムを実行してみよ。そして本日の演習の出席確認のため、実際に試してみたものの中から (最も印象的だったもの?) 1つを選び、そのプログラムを末尾のフォーム経由で提出せよ。

整数の「上限」と「下限」

ふつうの数学では、整数 $m$ に対して必ず $m < m+1$ であるし、また必ず $m > m-1$ でもあったはずだ。しかし計算機の中では、常にそうなるとは限らない。

下の2つの例は、整数 $m$ の値を変えながら $m+1$ や $m-1$ を計算して $m$ との大小関係を比較し、$m > m+1$ あるいは $m < m-1$ になるまで動き続けるというプログラムである。常識的に考えればこのプログラムの実行が止まる (すなわち $m > m+1$ あるいは $m < m-1$ という条件が満たされる) ことはあり得ないはずなのだが、実際にこの2つのプログラムを入力・コンパイル・実行してみることにより、本当に $m > m+1$ あるいは $m < m-1$ という異常な条件が満たされてしまう場合がある ことを確認してみよ。またそのときの $m$ の値がいくつであったかを調べてみよ。なおこれを含め、本日の演習のプログラムをコンパイル・実行するには例えば gfortran enshu10a1.f90 として a.out を作った後に ./a.out とすれば OK。ただしこのプログラムの実行が終了するまでには (ちょっと長めの) 時間がかかる可能性もあるので、本当に終わるまで気長に待つ ようにしておいてください。

$m > m+1$ になり得る例 (コピペ可)

program enshu10a1
  implicit none
  integer :: m,m_tasu1
  m = 1
  do
     m_tasu1 = m+1
     if (m > m_tasu1) then
        write(*,*)"m > m+1 when m= ",m,m_tasu1
        exit
     end if
     m = m_tasu1
  end do
end program enshu10a1

$m < m-1$ になり得る例 (コピペ可)

program enshu10a2
  implicit none
  integer :: m,m_hiku1
  m = 0
  do
     m_hiku1 = m-1
     if (m < m_hiku1) then
        write(*,*)"m < m-1 when m= ",m,m_hiku1
        exit
     end if
     m = m_hiku1
  end do
end program enshu10a2

Fortran 90 における整数 (integer) 型の変数には4バイト (=32ビット) の容量が与えられている。ということで、取り扱い可能な整数は $-2^{31}=-2147483648$ から $2^{31}-1=2147483647$ までの (0を含む) $2^{32}=4294967296$ 個に限られる。上のプログラムでは、この範囲を超える整数を扱わねばならなくなった場面で「異常」が発生し、計算が正しく実行できなくなってしまったのであった。

もはや単なる昔話でしかないかも知れないが、「2ギガバイトの壁」なんてのがあったなぁ。例えば昔の Windows でよく使われていたハードディスクのフォーマットだと、「2ギガバイトより大きなファイルを作成できない」なんてことが普通にあったものでしたわ。

実数の表現誤差と丸め誤差

計算機の中で実数を扱う際にも、実数が「有限の桁数でしか表現できない」ことや「10進法で表現されていない」という制約が顔を出してくることもある。

下のプログラム (コピペ可) は、ある自然数 $n$ を与え (この例では $n=2$)、その $n$ に対して $\displaystyle\sum_{i=1}^{100n}\frac{1}{n}=100$ という計算を実行し、その結果を表示するものである。自然数 $n$ の値を1から20程度の範囲で変えてプログラムを実行し、求まった総和の値がどのくらい正確かを調べてみよ。 $n$ の値によっては、総和の値が厳密に100だと求まる場合と求まらなかった場合があるはずだが、総和の値が厳密に100と求まる $n$ の値にはどのような特徴があるか。

program enshu10b
  implicit none
  integer :: i,n
  real :: souwa,atai
  n = 2
  atai = 1.0/real(n)
  souwa = 0.0
  do i = 1,n*100
     souwa = souwa + atai
  end do
  write(*,*)n,atai,souwa,atai*real(n*100)
end program enshu10b

Fortran 90 における実数 (real) 型の変数には4バイト (=32ビット) の容量が与えられており、2進23桁浮動小数点数 として表示される。このような実数として $\frac{1}{n}$ の値が (表現誤差もなく) 厳密に表わすことのできる $n$ は2のべき乗に限られてしまう。上のプログラムで総和が厳密に求まる場合が非常に少ないのはこれが理由であった。しかも $\frac{1}{n}$ が厳密に表現できていない場合には、$\frac{1}{n}$ を足す演算を1回実行するたびに、結果の数値を有効な桁数の範囲に「丸める」操作が行われた結果、「丸め」による誤差がどんどん累積し、正解から大きくずれた総和の値が求まってしまう。

またこのような計算にあたって、2進法と10進法の違いは意外な「盲点」ともいえるものである。実際、$\frac{1}{10}$ を2進法で表示しようとすると無限小数になってしまい、有限の桁数で表示しようとすれば必然的に表現誤差を含んでしまうのであった。

加減算で発生する「情報落ち」

ふつうの数学では、正の実数 $a$ に対して必ず $1+a > 1$ であったはずだ。しかし計算機の中では、常にそうなるとは限らない。

下の例は、正の実数 $a$ の値を変えながら $1+a$ を計算して 1 との大小関係を比較し、$1+a\le1$ になるまで動き続けるというプログラムである。常識的に考えればこのプログラムの実行が止まる (すなわち $1+a\le1$ という条件を満たす) ことはあり得ないはずなのだが、実際にこのプログラムを入力・コンパイル・実行してみることにより、本当に $1+a\le1$ という異常な条件を満たす場合がある ことを確認してみよ。またそのときの $a$ の値がいくつであったかを調べてみよ。

program enshu10c
  implicit none
  real :: a,ichi_tasu_a
  a = 1.0
  do
     ichi_tasu_a = 1.0+a
     if (ichi_tasu_a <= 1.0) then
        write(*,*)"1+a <= 1 when a= ",a,ichi_tasu_a
        exit
     else
        write(*,*)"1+a >  1 when a= ",a,ichi_tasu_a
     end if
     a = a*0.5
  end do
end program enshu10c

$a$ の値 (正確には $a$ の絶対値) が 1 と比べて十分大きいうちは問題ないのだが、$a$ の値が 1 と比べて小さくなり過ぎると、$1+a$ という計算の結果が「有効数字」の範囲からはみ出てしまうことになり、「$a$ を足す」という演算の情報が結果の数値からこぼれ落ちてしまうことになる。このような現象が「情報落ち」と呼ばれるものであった。

浮動小数点数において、「1より大きい最小の数」と1との差を「計算機イプシロン」あるいは「機械イプシロン (machine epsilon)」という。 IEEE 754 の (正規化された) 2進23桁浮動小数点数では、計算機イプシロンの値は $2^{1-(23+1)}\approx1.192\times10^{-7}$ である。上のプログラムでは、「$1+2a > 1$ だが $1+a\le1$ になってしまう $a$ の値」を求めることで、計算機イプシロン (のだいたいの値) を見積ろうとしていることに相当する。

計算機イプシロンは、次回以降の「数値シミュレーション」編の中でよく登場する「繰り返し計算」の手順を組み立てる上で重要な役割を持っている。例えば2つの実数 $x$ と $dx$ に対して、もし$\frac{|dx|}{|x|}$ の値が計算機イプシロンより小さくなっていたとしたら、 \[ \frac{|x+dx|}{|x|}\le\frac{|x|+|dx|}{|x|}=1+\frac{|dx|}{|x|}\approx1 \] であるから、もはや $x+dx$ と $x$ (の絶対値) の違いを比べる意味がなくなってしまうことになる。この考え方をさらに進めて、「$dx$ の絶対値と $x$ の絶対値の比 $\frac{|dx|}{|x|}$ が計算機イプシロン (の何倍か) よりも小さくとることができたら繰り返し計算を終える」というようなやり方がよく用いられている。

「単精度」の実数と「倍精度」の実数

ここまでの演習では実数として「ごく標準的な」ものを使用してきた。それらは実数型の変数に4バイト (=32ビット) の容量を与えるものであった。しかしこの「単精度 (single precision)」実数よりも大きな容量を用いて実数を表現することもある。例えば容量が8バイト (=64ビット) の「倍精度 (double precision)」実数も標準で用意されているし、容量が16バイト (=128ビット) の「4倍精度 (quadruple precision)」実数が用意されている環境もある。当然予想されるように、単精度の実数よりも倍精度の実数のほうが正確さが増すであろう。そこで上で用いたプログラム (enshu10c.f90) を再利用して、そのことを確認してみよう。

a と ichi_tasu_a という2つの単精度の実数型変数を宣言するには

あるいは (4バイトの容量がある、という意味で)

(

)

と書けばよい。これと同様にして、これら2つを倍精度の実数型変数として宣言するには、

(

)

あるいは4倍精度の実数型変数として宣言するには

(

)

とすればよい (説明はここ)。 enshu10c.f90 をこのように書き直した上で実行し、単精度実数の場合と倍精度実数 (あるいは4倍精度実数) で結果がどう違ってくるかを調べてみよ。

IEEE 754 の (正規化された) 倍精度実数は2進52桁浮動小数点数であり、この場合の計算機イプシロンは$2^{1-(52+1)}\approx2.220\times10^{-16}$ となる。この値は単精度実数の計算機イプシロンの値よりも8桁以上も小さくなっていることに注意。この意味を (大雑把に) いうと、倍精度の実数では単精度の実数に比べて有効数字の桁数が8桁ほど増えている ことになり、その分だけ計算の精度を高めることが可能になる。このような理由から科学技術計算、特に数値シミュレーション業界においては、「倍精度」の実数を標準的に使用することが多い (と思う)。

出席確認用プログラム提出フォーム

Moodle のコースページ以下の課題ファイル提出により、「離散数学」の実例の中で自分が最も印象に残った Fortran 90 プログラムを提出せよ。