ソフトウェアII 第2回（2023/12/07）

本日のメニュー

  • 構造体再確認
    • 構造体の基本
    • 構造体ポインタ
    • 構造体のサイズとアライメント
  • 物理シミュレーション: 多体問題
  • 課題について

本日の講義のSlidoは講義Slackで提示します。

準備

ダウンロード

まずは今日のサンプルをダウンロードしましょう。ダウンロードするファイルは今回共通ですが、そのあとにやる作業が環境ごとに異なります。

# wget の場合
wget https://www.gavo.t.u-tokyo.ac.jp/~dsk_saito/lecture/software2/resource/soft2-lec02-src.tar.gz
# curl の場合
curl -O https://www.gavo.t.u-tokyo.ac.jp/~dsk_saito/lecture/software2/resource/soft2-lec02-src.tar.gz

ダウンロード後は解凍して、そのディレクトリに入ってください。 今回の解凍先は、できればGoogle Driveのマウントポイント等ではなく、ローカルが望ましいです。

tar xvzf soft2-lec02-src.tar.gz
cd soft2-lec02-src

シンボリックリンクの作成

ここでは少しUnixコマンド的な話です。動作確認のための異なる環境毎のライブラリが同じディレクトリに同梱されています。まずはシンボリックリンクをはる という操作をしてみましょう。シンボリックリンクをはるとは、あるファイルに対する別名を設定する操作です。ファイルをコピーすると同じ内容の別ファイルが作られますが、シンボリックリンクは、「この名前の時はここにアクセスしてね」という情報を付与します。Mac で右クリックした場合の「エイリアスを作成」に相当します。Windows の場合はショートカットに似ていますが、微妙に違います。まずはディレクトリにあるファイルを確認してみましょう。

# .a のファイルになにがあるか表示する
ls *.a 
libgravity-apple-arm64.a	libgravity-cygwin.a		libgravity-mingw.a
libgravity-apple-x86_64.a	libgravity-linux.a

今回、1種類のライブラリ（libgravity-）が5通りの環境（Linux, Cygwin, MinGW, IntelMac, M1Mac）毎に用意されています。プログラムをコンパイルする際は同じファイル名だと解釈したいので、ここではコピーではなくシンボリックリンクを使ってみます。シンボリックリンクを作る一般的なコマンドは以下です。

ln -s [ターゲットになるファイル名] [その別名]

たとえば同じディレクトリ 内にあるA.txt というファイルにB.txt という別名をつけたい場合は以下になります。

ln -s A.txt B.txt

ls -al をすると

total 104
drwxr-xr-x   4 dsk_saito  wheel    128 11 25 11:48 .
drwxrwxrwt  21 root       wheel    672 11 25 11:48 ..
-rw-r--r--   1 dsk_saito  wheel  50000 11 25 11:51 A.txt
lrwxr-xr-x   1 dsk_saito  wheel      5 11 25 11:48 B.txt -> A.txt

のようにB.txt が A.txt の別名であることが示されています。B.txt のファイル容量は5B で、これは対象ファイル名のバイト数に対応します。A.txt 自体は50000B です。この操作は一つ上にあるディレクトリ等でも使えます。

ln -s ../A.txt B.txt

この時はB.txt は一つ上の階層のA.txtを指します。異なるディレクトリにあるファイルを指定する場合に、同じファイル名のシンボリックリンクを作る場合は、後ろの引数を省略できます。

# 今いるディレクトリにC.txt がシンボリックリンクとして作られる
ln -s ../C.txt

さて、今回のディレクトリ構成の場合だと、環境毎に異なるファイルを対象に、同じシンボリックリンク名にします。自身の環境に応じて実行してください。

# Intel Mac
ln -s libgravity-apple-x86_64.a libgravity.a
# Apple M1 Mac
ln -s libgravity-apple-arm64.a libgravity.a
# Linux (WSL2, GCS含む)
ln -s libgravity-linux.a libgravity.a
# Cygwin
ln -s libgravity-cygwin.a libgravity.a
# MinGW64
ln -s libgravity-mingw.a libgravity.a

もしコマンドを誤って貼り間違えた場合は、別名の方を削除します。

rm libgravity.a
Let's try

以下をやってみましょう。

  • 適当なテキストファイルを作ってシンボリックリンク作成を試してみましょう
  • ls -al で正しくリンクが貼られているか確認してみましょう

ここまでで環境準備ができたので、一旦授業のトピックに行きましょう。

構造体再確認

ソフトウェアIでも扱いましたが、構造体について再確認しましょう。

構造体とは

複数のデータ型をまとめあげ、新たなデータ型を作る仕組みです。まとめあげられたものを構造体と呼びます。 例えば以下は学生の名簿プログラムのようなものを考えた場合に、学生の必要情報をまとめあげた構造体です。

// id: 学生証番号相当の整数
// name: 名前
// age: 年齢。256歳でも大丈夫
// height: 身長
// weight: 体重

struct student {
    int id;
    char name[100];
    int age;
    double height;
    double weight; 
}; // 最後のセミコロンを忘れない

構造体の中の変数はメンバと呼びます。構造体のメンバへのアクセスには. を使います。

struct point {
    int x;
    int y;
};

int main() {
    struct point p1;
    
    p1.x = 10;
    p1.y = 20;
    
    return 0;
}

構造体の初期化や代入についてはソフトウェアIでも扱ったと思いますが、以下のstruct_init.c にまとめてあるので、確認しましょう。

// 構造体の初期化についてまとめました
#include <stdio.h>

// typedef と組み合わせてPoint型の宣言と定義を行います。
typedef struct point {
    double x;
    double y;
} Point;

int main(int argc, char **argv) {
    // 古いタイプの初期化
    Point p1 = { 1.0 , 2.0};
    
    printf("p1: (%3.1f %3.1f)\n", p1.x, p1.y);
    
    // 初期化指示子を用いた場合（こちらが可読性が高くおすすめ）
    Point p2 = {.x = 3.0, .y = 4.0};
    
    printf("p2: (%3.1f %3.1f)\n",p2.x, p2.y);
    
    // 複合リテラル: C99 から。
    // 構造体定数を(type-name){ initializer } で作成できる
    Point p3;
    
    // 一度宣言したあとは、p3 = { .x = 5.0, .y = 6.0}; だとエラー
    // (typename){} を使ってリテラルを代入
    p3 = (Point){ .x = 5.0, .y = 6.0};
    
    printf("p3: (%3.1f %3.1f)\n", p3.x, p3.y);
    
    // p3 に p2 を代入
    p3 = p2;
    
    printf("p3: (%3.1f %3.1f)\n",p3.x, p3.y);
    
    // 通常の配列（可変長でない）は、以下で初期化できる
    Point ps1[2] = { { .x = 1, .y = 2}, {.x = 3, .y = 4} };
    // **静的配列**の場合は、以下で配列サイズを取得可能
    size_t num_ps = sizeof ps1 / sizeof(Point);
    
    for (int i = 0 ; i < num_ps ; i ++)
        printf("ps1[%d]: (%3.1f %3.1f)\n", i,ps1[i].x, ps1[i].y);
    
    
    // 可変長配列の場合は、配列形式初期化が使えない
    size_t num = 2;
    // Point ps2[num] = { {.x=1,.y=1}, {.x=2,.y=2} }; // これはダメ
    Point ps2[num];
    ps2[0] = (Point){ .x = 1, .y = 2};
    ps2[1] = (Point){ .x = 3, .y = 4};
    
    for (int i = 0 ; i < num ; i ++)
        printf("ps2[%d]: (%3.1f %3.1f)\n", i,ps1[i].x, ps1[i].y);
    
    return 0;
}

構造体の宣言、定義とポインタ

構造体変数も単なる変数なので、その変数のアドレスを表すポインタを定義できます。

struct point {
    int x;
    int y;
};

int main(){
    struct point a;
    struct point *p = &a; // 構造体変数のアドレスをポインタに代入

    (*p).x = 10;
    p->x = 10; // アロー演算子（上と同じ意味でよく使う）

    return 0;
}

最後のアロー演算子は構造体ポインタを使う場合頻出の書き方です。よく覚えておきましょう。

構造体ポインタの使い方をいくつか見ておきましょう。代表的なものは関数呼び出しへの使用です。

void increment_age(struct student *s){
    s->age += 1;
}

void some_function() {
    // 構造体自体は、使う関数側で自動変数として確保
    struct student a =
      {
       .id = 1,
       .name = "Mike",
       .age = 21,
       .height = 175,
       .weight = 72 
      };
    // 構造体aが格納されているメモリのアドレスを渡す
    increment_age(&a);
}

もちろん構造体へのポインタではなく構造体を直接渡すことも可能です。ただしC言語の関数は全て値渡し（Call by Value）である という原則に従い、構造体のサイズ分のデータコピーが発生します。そのため大きな構造体を渡すときには注意が必要です。大きな構造体を関数に渡すときは原則ポインタの使用を考えましょう。

void print_age(const struct student *s){
    printf("age = %d\n", s->age);
}

上記の先頭にあるconst は「ポインタを渡すけどその参照先は書き換えちゃダメ」ということを示す修飾子です。これにより例えばs->age = 28; のような変更が関数内で書かれていた場合はコンパイラに怒られます。

冒頭の構造体のプログラムでは宣言と定義 を同時に行っています。実は定義（具体的にどういうメンバ変数を組み合わせるか）を除いた以下の記述も問題なく動きます。

struct point; // 宣言のみ。どういう構造体かは定義しない

typedef struct point Student; // typedefと組み合わせて Student型を宣言

定義のない構造体宣言に何の意味があるのかと思われますが、実は構造体のメンバに直接アクセスしてほしくない場合には宣言だけしておくことがあります。この場合には構造体を指すポインタだけを宣言し、このポインタを介して実際のメンバ変数とやりとりします。詳細は分割コンパイル を取り扱う際に述べますが、頭の片隅に置いておいて下さい。以下の例だと Student s; としてStudent型の変数を宣言するとエラーになりますが、Student *s; とStudent型へのポインタは問題なく宣言することができます。

typedef struct student Student;

int main(){
    Student s; // この宣言はStudent型にどれだけメモリを割り当てるかわからないのでダメ
    return 0;
}


// コンパイルに成功する例
typedef struct student Student;

int main(){
    Student *s; // この宣言はStudent型のアドレスを格納するポインタの宣言なのでOK
    return 0;
}

構造体のサイズとアライメント

構造体は複数の変数を組み合わせたものです。構造体を構成した場合は、メンバ変数は連続したメモリ領域に確保されます。このとき構造体のサイズがどうなるかを確認してみましょう。sizeof演算子は型名に適用した場合、その型のサイズを返します。下記のプログラムの実行結果はどのようになるでしょうか？

#include <stdio.h>

struct point1{
    int x;
    int y;
};

struct point2{
    int x;
    int y;
    char c;
};

int main() {
   printf ( "%zu\n", sizeof(struct point1));
   printf ( "%zu\n", sizeof(struct point2));
   return 0;
}
Let's try

上記プログラムを写経して実行してみましょう。

なぜこのような結果になるでしょうか。構造体のメモリ配置には変数のアライメントという概念が関係してきます。CPUがメモリからデータを読み出す際に、まとまった単位（チャンク : 64bitマシンでは通常8B）でやりとりします。そのために特定のメモリアドレスの倍数にデータが配置されるように調整することで性能が改善し、この処理をアライメント といいます。この処理はコンパイラが自動的に行いますが、その過程で、構造体に配置された変数がそのまま連続して並べると条件を満たさない場合に、パディング という処理で空のデータ領域を途中に挿入します。なおx86アーキテクチャ（ほとんどのみなさんが普段使っているCPUアーキテクチャ）やARMアーキテクチャの場合、データサイズと同じ単位（sizeof で得られるサイズ）でデータのアライメントが取られます。なお構造体のアライメントを考える場合はその構造体が配列として連続してメモリ確保されるケースを想定するとよいです。struct point2 の場合、char は1バイトアライメントなので、どのアドレスに設置しても良さそうですが、この構造体の配列を確保した場合に、次のint の先頭アドレスがint のアライメント（4バイト）からずれます。これが起きないようにコンパイラはchar のあとに3バイトのパディングを挿入しています。

Let's try

以下をやってみましょう。

  • 構造体の初期化・代入のサンプルコードであるstruct_init.c を参考に、1000個の要素を含むstruct point の配列から座標平均（=重心）を計算するプログラムを書いてみましょう。各点の座標はランダムに初期化してみてください。
  • 基本型の変数（char int double float など）を複数種類、宣言し、そのアドレスを表示してみてください。アドレスの末尾に着目してみましょう。
  • struct_alignment.c を修正し、定義された構造体のサイズを調べてみましょう。

物理シミュレーション: 多体問題

前回の講義ではライフゲームを進行するためのルール（条件）を実装し、シミュレーションを行いました。今回は物理法則を条件とする計算機シミュレーションを行ってみましょう。扱う問題は多体問題です

万有引力と多体系

電気系の学生にとっては、クーロン力の方が基本的な力であるとは思いますが、物理としてやはり考慮しなければいけないのは万有引力です。質量を有する2つの物体間には以下に示す力が働きます。

$$ F_1 = F_2 = G\frac{m_1m_2}{r^2} $$

\(m_i\)は物体\( i \)の質量、\(r\)は物体間の距離、\( G\) は重力定数です。クーロン力の定義に似てますね :-)

さてさらに多くの物体がある場合はどうでしょうか。多くの物体がある場合に、一つの物体に着目して他の物体からの影響を記述します。粒子\( i\) が粒子\( j\) から受ける重力のベクトルは

$$ {\boldsymbol F}_{ij} = G \frac{m_im_j }{|{\boldsymbol r}_j - {\boldsymbol r}_i|^2} \frac{({\boldsymbol r}_j - {\boldsymbol r}_i)}{|{\boldsymbol r}_j - {\boldsymbol r}_i|} $$

最終的に距離をまとめると

$$ {\boldsymbol F}_{ij} = G \frac{m_im_j }{|{\boldsymbol r}_j - {\boldsymbol r}_i|^3} ({\boldsymbol r}_j - {\boldsymbol r}_i) $$

となります。他の物体全てからを考慮して

$$ {\boldsymbol F} _{ij} = \sum _{j \ne i} {\boldsymbol F} _{ij}$$

で算出できます。

運動方程式

物理未履修者もいるので、確認しておくと物体の動きを微分方程式で記述したのが運動方程式です。位置の時間微分が速度ですが、速度と力の関係が微分方程式で記述されます。

$$ \frac{d\boldsymbol r _i}{dt} = \boldsymbol v _i$$

$$ m_i\frac{d\boldsymbol v _i}{dt} = \boldsymbol F _i $$

この右辺に先ほどの他の物体からの力を代入すると、速度の微分は

$$ \frac{d\boldsymbol v _i}{dt} = G \sum _{j \ne i} \frac{m_j }{|{\boldsymbol r}_j - {\boldsymbol r}_i|^3} ({\boldsymbol r}_j - {\boldsymbol r}_i)$$

となります。この微分方程式をコンピュータで計算してみましょう。

オイラー法

コンピュータは人間のような微分方程式の解き方はできません。ここでは微分方程式を近似して数値計算で解く方法のひとつであるオイラー法を紹介します。

微分の定義に戻ると関数\( x(t) \) の微分について以下が成り立ちます。

$$ x'(t) = \lim_{\Delta t \rightarrow 0} \frac{x(t+\Delta t)-x(t)}{\Delta t}$$

極限を取る前の式を変形すると

$$ x(t+\Delta t) \simeq x(t) + x'(t) \Delta t $$

この式から出発して時間が\( \Delta t\) 毎に進行すると仮定すると、前の時刻の加速度、速度をもとにした更新式を導けます。 基本的な手順は以下のようになります。

  • 微小時間\( dt\) に対する加速度の更新式に基づき加速度の値を更新

$$ a_i^{(t)} = G \sum _{j \ne i} \frac{m_j}{|r_j^{(t)}-r_i^{(t)}|^3}(r_j^{(t)}-r_i^{(t)})$$

  • 加速度と速度、微小時間の関係から速度の値を更新

$$v_i^{(t+1)} = v_i^{(t)} + a_i^{(t)}\Delta t$$

  • 速度と位置、微小時間の関係から位置を更新

$$r_i^{(t+1)} = r_i^{(t)} + v_i^{(t)}\Delta t$$

サンプルのコンパイル

サンプルプログラムをコンパイル、実行してみましょう。

# サンプルプログラム
gcc -o gravity -Wall main.c -L. -lgravity -lm
./gravity

コンパイルオプションはほぼ前回と同じ構成ですが、数学関数を使うために最後に-lm をつけています。

全体像

main.c の全体構造は以下のようになっています。

  • gravity.h: 今回の構造体の定義と関数プロトタイプ宣言
  • libgravity.a: プロトタイプ宣言で宣言されている関数の実装
  • main.c: main 関数の記述
Let's try

以下に取り組んでみましょう。

  • 物理パラメータや座標、速度を変化させてコンパイルし、遊んでみましょう。
  • メイン関数のあるソースを読んで全体像を把握してみましょう。
  • ヘッダファイルを読んでどのような関数が宣言されているか確認しましょう。
  • 全体像を理解した上で、実装すべき関数をmy_????? で実装しましょう。
  • ヘッダファイル内の使わなくなったプロトタイプ宣言の削除
  • -L. および -l???? を除いてコンパイル: -lm は残す

気を付けるポイントをあげます

  • プロトタイプ宣言はgravity.hに書いて問題ありませんが、関数の実装はmain.c に書いてください。
  • 構造体定義を変更するとライブラリが動かなくなります。

課題（締切: 12/20）

  1. struct_alignment.c の3つの構造体について、構造体中のメンバ変数に割り当てられているアドレスを表示するプログラムを作成せよ。

    • ファイル名はstruct_alignment.c のままでよい。
    • 講義でやった変数の型毎のアライメント等をもとにどのようなパディングが起きているかを考察しstruct_alignment.md に記して同梱せよ。

  2. まずはstruct_compare.tar.gz をダウンロードし展開する（以下参照）。このディレクトリにあるstruct_compare.c をコンパイルし、結果を確認せよ。このプログラムは構造体にある細工を行う関数trick を用いている（この関数はlibobj-*.a にある）。細工の前後で起きた変化を確認した上で以下に取り組むこと（考察はすべてstruct_compare.mdに記載し同梱する）。

    • 何が起きているかを説明せよ。
    • 上記の現象から構造体を比較する際になぜメンバを逐次比較する必要があるかについて記載せよ。
    • struct_compare.c にコードを追加して、trickが行った細工 を見つけよ。
    • コードを追加したstruct_compare.c を同梱すること。

  3. gravity.c をもとに今回の物理シミュレーションについて2次元空間を扱えるように拡張せよ。

    • 物体数を３つ以上に増やして動作を確認すること
    • 例えば円運動などを例示できるとよい
    • ファイル名は my_gravity1.c
    • 本課題についてはmy_???? への実装が終わってから取り組むこと。課題の実装にあたって構造体の定義を変更する必要があるため、中途半端にとりかかると事前に同梱したライブラリが機能せずバグります。

  4. 課題3を実装した上で、第一引数にオブジェクトの個数、第二引数に2次元座標と質量が記載されたファイル（以下参照）から物体の設定を読み取れるようにせよ。

  • 例えば ./my_gravity2 3 data.dat のように実行できるようにする。与えられた個数よりdata.dat 内の物体数が多い場合も与えられた個数分の物体だけ読み取るようにすればよい。逆に少ない場合は、0質量の物体を生成するなどして対応する。
  • ファイル名はmy_gravity2.c
  • 自身で試したオブジェクトファイルも同梱せよ。拡張子はdatとする。実体はテキストファイル。
  1. 物体同士の衝突（融合）現象を実装せよ。
    • 星同士の距離がある値以内になったら融合するとしてよい（運動量保存）
    • ファイル名はmy_gravity3.c
  2. ［発展課題］　本アプリケーションをさらに発展させよ。
    • ファイル名はmy_gravity4.c
    • 発展のさせ方は任意
      • 中点法やルンゲクッタ法等によるシミュレーションの精度向上
      • 惑星配置
      • 万有引力つきビリヤード（壁や反射係数の導入）
      • 衝突による破壊・散乱を考慮、電界・磁界を考慮
      • etc

課題2について

課題2にあたり、まずは以下の通りstruct_compare.tar.gz をダウンロードする。

# wget の場合
wget https://www.gavo.t.u-tokyo.ac.jp/~dsk_saito/lecture/software2/resource/struct_compare.tar.gz
# curl の場合
curl -O https://www.gavo.t.u-tokyo.ac.jp/~dsk_saito/lecture/software2/resource/struct_compare.tar.gz

その後展開し、環境に応じてシンボリックリンクを貼る。

# 展開
tar xvzf struct_compare.tar.gz
cd struct_compare
# 該当する環境のライブラリにシンボリックリンクを貼る
ln -s libobj-apple-arm64.a libobj.a
# コンパイル
gcc -o struct_compare struct_compare.c -L. -lobj
# 実行
./struct_compare

課題4ファイル形式

  • 先頭が# で始まる行はコメント行とする
  • 一行に一つの物体
  • 空白区切りで質量m、x 座標、y座標、x初速、y初速の順で記述
  • 必要情報のあとの # 以下もコメントとして記述可能
# m x y vx vy
60.0 0.0 -20.0 0.0 2.0 # 小さい玉
100000.0 0.0 1000.0 0.0 0.0 # 地球

提出方法

  • ITC-LMS にて提出する
  • 全てのプログラム/ファイルをまとめ、zip またはtar.gzで圧縮し提出
  • 必ずREADME.md を同梱し、やった内容を簡潔に記述する
  • ファイル名はSOFT-12-07-NNNNNNNN.zip または SOFT-12-07-NNNNNNNN.tar.gz とする
    • NNNNNNNN部分は学籍番号（ハイフンは除く）
    • JについてはJ??????? のようにする
  • 課題について
    • 基本課題は毎回提出
    • 発展課題は成績計算に全6回中上位3回分を採用する