Skip to content

第3回 データ形式、再帰関数、アルゴリズム3

本日の講義内容

  • データ形式
  • 再帰関数
  • 組み合わせ最適化

データ形式

これまで、ファイルを介したデータの入出力では、人間が読むことを前提としたテキストデータを扱ってきました。今回は、バイナリデータの読み書きを扱います。

バイナリデータという語を字義通りに解釈すると、2進数表現で記述されたデータ、を意味します。この意味では計算機が扱うすべてのデータは本質的にバイナリデータであると言えます。

広義にはテキストデータもバイナリデータの一種であるものの、通常、テキストデータとは、人間が直接読んで理解できる文字集合（ASCIIやUnicodeなど）を、UTF-8などの符号化方式で表現したデータを指します。

テキストデータとバイナリデータ

一方でバイナリデータは、特定の文字表現を前提とせず、バイト列そのものの構造によって意味が定まるデータのことを指すことが多いです。人間が一見して内容を読み取ることは困難ですが、計算機にとっては効率よく処理できる形式です。

特に、数値データの読み書きに着目すると、テキストファイルではなくバイナリファイルを用いることで、いくつかの重要な利点があります。

  • 文字列から数字への変換に伴う精度劣化の可能性がない
  • ファイルサイズが縮小される
  • 通常テキストを読み取るより速度面で速い

テキストデータの入出力

バイナリデータに先だって、テキストデータの入出力方法についておさらいしましょう。まずファイルを開いたり、閉じたりする関数として、fopen()やfclose()がありました。

#include <stdio.h>
FILE *fopen(char *name, char *mode);  // ファイルを開く
int fclose(FILE *fp);  // ファイルを閉じる

また、これまで標準入力および標準出力をベースにした入出力関数としてgetcharやputchar、printfやscanfが紹介されてきました。これらの関数にはそれぞれのFILE型ポインタバージョン、char*バージョンと呼べるものが存在します。

標準入出力バージョン FILE*バージョン char*バージョン
一文字ずつ入力 getchar fgetc / getc 配列アクセス
一文字ずつ出力 putchar fputc / putc 配列アクセス
一行ずつ入力 gets fgets 強いて言えば (strtok)
char*を出力 puts fputs 強いて言えば (strcpy)
フォーマット入力 scanf fscanf sscanf
フォーマット出力 printf fprintf sprintf

FILE型ポインタを用いた入力例を以下に示します。

readfile.c
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
#include <stdio.h>
#include <stdlib.h>
#include <string.h>


int main(int argc, char **argv) {
    const size_t bufsize = 500;
    char buf[bufsize];

    if (argc != 2){
        fprintf(stderr,"usage: %s filename\n",argv[0]);
        exit(EXIT_FAILURE);
    }

    FILE *fp;

    if ((fp = fopen(argv[1], "r")) == NULL) {
        fprintf(stderr, "error: can't open %s\n", argv[1]);
        exit(EXIT_FAILURE);
    }

    while (fgets(buf, bufsize, fp) != NULL) {
        size_t len = strlen(buf) - 1;
        printf("%zd\n", len);
    }

    fclose(fp);

    return EXIT_SUCCESS;
}

演習

  • 上記のプログラムを写経し、コンパイルしてみましょう。
  • テキストファイルを作成し、それを引数として実行してみましょう。

バイナリデータの入出力

続いて、バイナリデータの入出力方法について確認してみましょう。ファイルを開いたり、閉じたりする関数についてはテキスト、バイナリ共通で、同様です。

#include <stdio.h>
FILE *fopen(char *name, char *mode);  // ファイルを開く
int fclose(FILE *fp);  // ファイルを閉じる

ただし、ファイルを開く際の読み書きのモードについて、バイナリを扱う場合にはbを追加してrb（読み込みの場合）またはwb（書き込みの場合）としておく方がよいです。全てのOSで適切に動作する可搬性の高いコードという観点から

  • テキストの読み書きの場合はbをつけない
  • バイナリの読み書きの場合はbをつける

と覚えておきましょう。

Mac OSやLinuxではbの有無によって結果が変わることはありませんが、Windowsの場合は注意が必要です。Windows の場合、bが存在しない場合はテキストモードと判断し、改行の置換と0x1aをEOFとみなすという処理が入ります。bをつけるとバイナリモードとなり、変換なしでファイルを読み書きします。このようにWindowsの場合のみbの有無によってファイルの読み書きの挙動が異なるため注意しましょう。

コラム: Windowsのテキストファイル仕様

現在のMac OSやLinuxでは、改行コードは\nとなっているのでした。一方、 Windowsでは改行コードは\r\nとなっています。Windowsの場合テキストモードでは、この改行コードを\nへと変換します。なお、こうしたコードは、遡ると昔のタイプライタの動作 (carriage return (CR)、もともとタイプライタのヘッドを先頭に戻す機能とline feed (LF)、行送り機能) に対応して決まっているものと思います。こうした仕様を昔のOSが受け継いでいまに残っています。同様に、0x1aをEOFとみなすという処理も昔のOS由来の仕様です。興味のある方は、MS-DOSやCP/Mといったものについて調べてみるとよいでしょう。

バイナリファイルを直接読みとるときにポイントになるのは

  • どのファイルから読み取るか
  • ファイルから読み取ったデータをメモリのどこに書き込むか
  • どのサイズのデータをいくつ読み取るか

です。一方、変数や配列などのデータをファイルに書き込む場合は

  • どのファイルへ書き込むか
  • メモリのどこを読み取ってファイルに書き込むか
  • どれぐらいのサイズを読み取って書き込むか

の3点です。これらの情報をもとにバイナリの読み書きをするのが以下のfread()とfwrite()です。

#include <stdio.h>

// ファイルからsizeバイトのデータnitems個分を指定のアドレスに読み込む
size_t fread(void *data, size_t size, size_t nitems, FILE *fp); 

// 指定アドレスからsizeバイトのデータnitems個分を読み取り、ファイルに書き出す
size_t fwrite(void *data, size_t size, size_t nitems, FILE *fp);

たとえば16バイトのデータを読み取りたい場合、4バイトのデータが4つなのか、8バイトのデータが2つなのかによって、呼び出し方は変わりますが、読み出された後の先頭アドレスから16バイトのデータそのものには変化はありません。例えばその先頭アドレスを受け取るポインタの型が異なれば、異なるデータ型で異なる個数のデータとして解釈可能です。

なお、それぞれの関数の返り値は「読み込み/書き込みに成功した個数」となります。この情報を利用することで例えば「8バイトデータを3個ずつ順に読み取ろうとしたが最後だけ2個しか読み取れずファイル終端を迎えた」などの状況にも適切に対応することができます。

バイナリファイルの読み込み

まずは簡単な例としてテキストファイルをfread()を使って読み取りましょう。テキストファイルも単なるデータの並びとして扱うことができます。

readbinfile.c
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
#include <stdio.h>
#include <stdlib.h>
#include <errno.h>

int main(int argc, char **argv) {
    if (argc != 2) {
        fprintf(stderr, "usage: %s filename\n", argv[0]);
        return EXIT_FAILURE;
    }
    FILE *fp;
    if ((fp = fopen(argv[1],"rb")) == NULL) {
        // fopen は失敗した場合 errno の変数をセットする
        // perror は errno の番号に対応したメッセージを返してくれる
        perror(argv[1]); 
        return EXIT_FAILURE;
    }
    // 100バイトのchar配列（スタック領域）
    char buf[100];

    // 1バイトを100個分読み出す。rsize には読み出した個数が格納される
    size_t rsize = fread(buf, sizeof(char), 100, fp); 

    printf("%zu Byte read\n", rsize);

    // 文字列だと解釈したい場合は終端記号\0を追加
    // この作業は必須ではない
    if (rsize < 100) {
        buf[rsize] = '\0';
        printf("%s",buf);
    }    

    return EXIT_SUCCESS;
}

fread()は先頭アドレスを引数にとります。このメモリ領域は事前に確保されている必要があります。上記の例ではスタック領域に確保されています。動的に読み込むサイズを変更する場合はmallocを使ってヒープ領域に確保してから使用します。

演習

上記のプログラムを写経・コンパイルし、以下のhelloworld.txtをファイルとして保存したものを引数として実行してみましょう。

helloworld.txt
Hello,World

バイナリファイルの書き込み

以下のwritebinaryfile.cは1000万個のdoubleデータを用意しテキスト、バイナリの2通りに出力するプログラムです。コンパイルして実行し、できた2つのファイルサイズを比較してみましょう。

writebinaryfile.c
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
#include <stdio.h>
#include <stdlib.h>

int main(int argc, char **argv) {
    // 1000万要素作成したdoubleをテキスト/バイナリファイルで出力してみる
    // 先頭に個数情報を出力しておく
    // バイナリの場合はsizeof(size_t)分を最初に出力する
    if (argc != 3) {
        fprintf(stderr, "usage: %s <txt filename> <binary filename>\n",
                argv[1]);
        return EXIT_FAILURE;
    }

    size_t size = 10000000;
    double *d = (double *)malloc(sizeof(double) * size);  // とりあえず1000万確保

    srand(100);
    for (int i = 0; i < size; i++) 
        d[i] = 0.5423 * rand();

    // テキストに出力してみる
    FILE *fp;
    if ((fp = fopen(argv[1], "w")) == NULL) {
        fprintf(stderr, "%s: cannot open file.\n", argv[1]);
        return EXIT_FAILURE;
    }

    fprintf(fp, "%zu\n", size);
    for (int i = 0; i < size; i++) 
        fprintf(fp, "%f\n", d[i]);
    fclose(fp);

    // "wb"の'b'はwindowsは必須, linux/macはなくてもよいがここでは"wb"としておく
    if ((fp = fopen(argv[2], "wb")) == NULL) {
        fprintf(stderr, "%s: cannot open file.\n", argv[1]);
        return EXIT_FAILURE;
    }
    // 最初にsizeof(size_t) 1個分をサイズ情報として出しておく
    fwrite(&size, sizeof(size_t), 1, fp);
    // dポインタを先頭にからsize個のdoubleデータを出力
    fwrite(d, sizeof(double), size, fp);
    fclose(fp);
    free(d);

    return EXIT_SUCCESS;
}

演習

上記のプログラムを写経・コンパイルし、以下のように実行してみましょう。 

# 第1引数がテキストファイル名、第2引数がバイナリファイル名
$ ./writebinaryfile data.txt data.dat

テキストファイルとバイナリファイルをそれぞれ出力することができたら、読み取りについても手を動かして試してみましょう。

演習

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

  • テキストファイルdata.txtを読み取り、double型の配列にセットするプログラムを書いてみましょう。
  • バイナリファイルdata.datを読み取り、double型の配列にセットするプログラムを書いてみましょう。
    • ヒント: fread()を使って読み取り、先頭の個数情報を利用して適切にmallocする
  • 実行時間を比較してみましょう
    • 実行時間計測にはプログラムの前にtimeとつけます。UNIXにおけるtimeコマンドはプログラムの実行時間を計測するプログラムです。 
      $ time ./readtextfile data.txt
こたえ

readtextfile.c
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
#include <stdio.h>
#include <stdlib.h>

int main(int argc, char **argv){

    if (argc != 2){
        fprintf(stderr, "usage: %s <filename>\n", argv[0]);
        return EXIT_FAILURE;
    }

    FILE *fp;
    if ( (fp = fopen(argv[1],"r")) == NULL ){
        perror(argv[1]);
        return EXIT_FAILURE;
    }

    size_t size = 0;
    char buf[256];

    // サイズ情報の読み取り
    fgets(buf,256,fp);
    sscanf(buf,"%zu%*1[\n]", &size);

    // 読み取ったサイズ情報に基づいてmalloc
    double *data = (double*)malloc(sizeof(double)*size);
    size_t i = 0;

    // 1行ずつ読み取って値をセット
    while(fgets(buf,256,fp) != NULL){
        sscanf(buf,"%lf%*1[\n]",data + i++);
    }

    fclose(fp);
    free(data);

    return EXIT_SUCCESS;
}
readbinaryfile.c
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
#include <stdio.h>
#include <stdlib.h>

int main(int argc, char **argv){
    if (argc != 2 ){
        fprintf(stderr, "usage: %s filename\n", argv[0]);
        exit(1);
    }

    FILE *fp;
    if ((fp = fopen(argv[1],"rb")) == NULL){
        perror(argv[1]);
        exit(1);
    }

    size_t datasize = 0;
    fread(&datasize, sizeof(size_t),1,fp);

    double *data = (double*)malloc(sizeof(double)*datasize);
    fread(data, sizeof(double),datasize,fp);
    fclose(fp);
    free(data);

    return EXIT_SUCCESS;
}

バイナリファイル読み書きの具体例

実際に世の中に存在するファイル形式では「ヘッダ部分」と「データ部分」が存在することがほとんどです。たとえば音声ファイル形式の一つであるWAVファイルは、通常先頭44バイトにファイルサイズやチャンネル数などの情報が規定の方法で存在し、その後実際の音声データが所定のバイナリ（short型整数やfloat型の実数）で保存されています。それぞれのファイルの入出力を自作する場合は

  • ヘッダ部分を読み取り、その後の実データを読み取る準備をする（mallocによるメモリ領域確保や構造体へのヘッダ情報の代入など）
  • 実データ部分を読み取り、確保したメモリ領域にセットする

の機能を満たすように設計することになります。

演習

以下の音声ファイルをバイナリファイルとして読み込み、44バイトのヘッダ部分の情報を標準出力してみましょう。

🎵 voice.wavをダウンロード

  • クレジット: VOICEVOX:ずんだもん

  • 環境上ブラウザからのダウンロードが面倒な場合はコマンドラインから、 

    $ curl -O https://eeic-software2.github.io/2025/assets/voice.wav

  • 以下のような構造体を使い、それぞれの内容を標準出力してみましょう。

Wavheader
#include <stdio.h>
#include <stdlib.h>
#include <errno.h>

// WAVヘッダ用構造体
typedef struct {
    char riff[4];  // "RIFF"というヘッダ文字列
    int filesize;  // ファイルサイズ - 8 (byte)
    char wave[4];  // "WAVE"という文字列

    char fmt[4];  // "fmt"という文字列
    int fmtsize;  // リニアPCMなら16 (byte)       
    short format;         
    short channels;  // モノラルは1、ステレオは2
    int samplerate;  // サンプリング周波数 (Hz)
    int byterate;
    short blockalign;
    short bitspersample;

    char data[4];  // "data" という文字列
    int datasize;  // データサイズ (byte)
} WavHeader;
こたえ
readvoice.c
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
#include <stdio.h>
#include <stdlib.h>
#include <errno.h>

// WAVヘッダ用構造体
typedef struct {
    char riff[4];  // "RIFF"というヘッダ文字列
    int filesize;  // ファイルサイズ - 8 (byte)
    char wave[4];  // "WAVE"という文字列

    char fmt[4];  // "fmt"という文字列
    int fmtsize;  // リニアPCMなら16 (byte)       
    short format;         
    short channels;  // モノラルは1、ステレオは2
    int samplerate;  // サンプリング周波数 (Hz)
    int byterate;
    short blockalign;
    short bitspersample;

    char data[4];  // "data" という文字列
    int datasize;  // データサイズ (byte)
} WavHeader;

int main(int argc, char **argv) {
    if (argc != 2) {
        fprintf(stderr, "usage: %s filename\n", argv[0]);
        return EXIT_FAILURE;
    }

    FILE *fp;
    if ((fp = fopen(argv[1], "rb")) == NULL) {
        perror(argv[1]);
        return EXIT_FAILURE;
    }

    // 読み込み
    WavHeader header;
    if (fread(&header, sizeof(WavHeader), 1, fp) != 1) {
        perror("fread");
        fclose(fp);
        return EXIT_FAILURE;
    }

    // ヘッダ内容を表示
    printf("RIFF header      : %.4s\n", header.riff);
    printf("File size field  : %d byte\n", header.filesize);
    printf("Format           : %.4s\n", header.wave);

    printf("fmt chunk ID     : %.4s\n", header.fmt);
    printf("fmt chunk size   : %d byte\n", header.fmtsize);
    printf("Audio format     : %d\n", header.format);
    printf("Channels         : %d\n", header.channels);
    printf("Sample rate      : %d Hz\n", header.samplerate);
    printf("Byte rate        : %d\n", header.byterate);
    printf("Block align      : %d\n", header.blockalign);
    printf("Bits per sample  : %d\n", header.bitspersample);

    printf("data chunk ID    : %.4s\n", header.data);
    printf("Data size        : %d byte\n", header.datasize);

    fclose(fp);
    return EXIT_SUCCESS;
}

バイトオーダ

バイナリデータの読み書きに関連して、計算機やネットワークにおけるビット表現、バイト表現がどのようになっているかについて少し学んでおきましょう。以下のプログラムはunsigned char型（1バイト）の変数に1を代入し、そのビット表現を出力しています。

print_bit.c
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
#include <stdio.h>

int main(void){
    unsigned char a = 1;
    for (int b = 0 ; b < 8 ; b++) {
        printf("%d", (a & (1 << (7-b))) ? 1 : 0);
    }
    printf("\n");
    return 0;
}

次にunsigned short型（2バイト）の変数に1を代入したあと、そのアドレス領域をunsigned char型2つ分と思って読むとどうなるでしょうか?

print_bit2.c
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
#include <stdio.h>

int main(void){
    unsigned short a = 1;

    // unsigned char型ポインタを用意する
    // 上記の先頭アドレスをキャストして代入
    unsigned char *p = (unsigned char*)&a;

    for (int i = 0 ; i < sizeof(unsigned short) ; i++) {
        printf("%p\n", p+i);
        for (int b = 0 ; b < 8 ; b++){
            printf("%d", (p[i] & 1 << (7-b))?1:0);
        }
        printf("\n");
    }
    return 0;
}

演習

上記のプログラムそれぞれを写経して、コンパイル・実行してみましょう。

一つ目のプログラムの結果は右端に1ビット存在し、整数1の2進数表現として妥当な結果が得られているかと思います。2つ目のプログラムの結果はどうでしょう。この結果はCPUに依存しますが、おそらくほとんどの皆さんが使っていると考えられるAMDやIntelのCPU（x86）環境、最近のAppleのCPU (ARM) 環境では、最初のバイトの方に1が立ち、2番目のバイトが0になっているかと思います。このように1バイトを超える多バイトのデータをどの順番で格納したり伝送するかをバイトオーダと言います。多くのプロセッサの場合は、リトルエンディアンと呼ばれる方式で、最下位ビットの属するバイト（Least Significant Byte: LSB）が若いアドレスに格納されます。逆に最上位ビットの属するバイト（Most Significant Byte: MSB）から順に若いアドレスに格納する方式をビッグエンディアンといいます。TCP/IPのヘッダなど、ネットワーク周りではビッグエンディアンが採用されていることが多いです。

エンディアン

バイトオーダは普段のプログラミングで意識することは少ないですが、ネットワークデータを扱ったり、ファイル形式によっては特定のエンディアンを前提にしているケースがあります。計算機へのデータの格納方式の違いとして頭の片隅に入れておくとよいでしょう。

浮動小数点数の形式

浮動小数点数を様々な形式で表示してみましょう。C言語の浮動小数点は実は様々な形式で指定することができます。C言語の浮動小数点はIEEE754と呼ばれる標準規格に基づいています。特に16進数浮動小数点表記を用いるとそのビット構造との対応が非常にわかりやすくなります。10進数表記では厳密には表せないものを表記可能です。

print_float.c
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
#include <stdio.h>

int main(void){

    // IEEE 754 単精度(float, 32bit) の浮動小数点数、内部では近似値になる
    float f = 0.0000102411234;

    // %f:
    // 固定小数点表記（小数点以下6桁がデフォルト）
    printf("%%f: %f\n", f);

    // %.8f:
    // 小数点以下8桁まで表示
    printf("%%.8f: %.8f\n", f);

    // %e:
    // 指数表記（科学技術表記）
    printf("%%e: %e\n", f);

    // %.8e:
    // 指数表記で仮数部を8桁表示
    printf("%%.8e: %.8e\n", f);

    // %a:
    // 16進浮動小数点表記
    // IEEE 754 の内部表現（仮数×2^指数）をそのまま表示
    printf("%%a: %a\n", f);

    // %.8a:
    // 16進浮動小数点表記で仮数部の表示桁数を指定
    printf("%%.8a: %.8a\n", f);

    return 0;
}

演習

上記のプログラムをコンパイル・実行してみましょう。

再帰関数

再帰呼び出し

ある手続きの中で自分自身を呼び出すことを再帰呼び出しと言います。C言語の場合は関数で手続きを表現するので、引数を変えながら関数内で自分自身を呼び出すのが再帰です。以前の木構造やペイントソフトでも、この再帰呼び出しを活用していました。一例として、以下のような階乗を求める関数を考えます。

factorial.c
int factorial(int n){
    if ( n == 0 ) return 1;

    return n * factorial(n-1);
}

この関数の実行過程は以下のようになります。関数のreturn部分で、新しい引数を利用して関数が呼ばれていきます。終端で値が確定すると、次々と呼び出し元に戻り、最終的に最初の呼び出しの値が確定します。 階乗の実行過程

再帰の特徴

再帰呼び出しは基本的には等価なループ処理で書くことが可能です。また実行速度面では繰り返しで記述した方が高速な場合が多いです。しかし一方、再帰呼び出しを用いることでアルゴリズムを簡潔に実装できる場合があります。たとえば以下はユークリッドの互除法で最大公約数を計算するものです。再帰関数を用いることで簡潔な記述を実現できていることがわかります。

// while version of euclid gcd
int euclid(int m, int n) {
    int remainder;

    while ( (remainder = m%n) != 0) {
        m = n;
        n = remainder;
    }
    return n;
}

// recursion version of euclid gcd
int euclid(int m, int n) {
    int remainder = m % n;
    return (!remainder) ? n : euclid(n, remainder);
}

再帰関数の停止条件

再帰を実装する場合に注意すべき点として、必ず再帰を停止する部分を書き、かつ再帰呼び出しのたびに入力が停止条件へ近づくようにするということがあります。このような停止条件を書かない場合、無限に関数が呼び出されることになります。

試しに止まらない再帰を書いてみましょう。どうなるでしょうか。

recursive.c
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
#include <stdio.h>

void test (int *a){
    printf("%d\n", (*a)++);
    test(a);
}

int main(void){
    int a = 0;
    test(&a);
    return 0;
}

演習

上記のプログラムを写経・コンパイルし実行してみましょう。

上記のプログラムは呼び出した回数を表示しながら最終的にSegmentation faultで終了すると思います。これは前回の講義で説明したスタック領域のサイズに関係します。再帰関数を呼び出す際には、呼び出す際の変数をコピーし、かつ呼び出した関数の終了後に戻ってくるアドレスを順次記録していきます。この情報はスタック領域に積み上げられていきます。

スタックオーバーフロー

正常に停止条件に辿り着いた場合には、関数の呼び出しが終了するたびに、それぞれの関数呼び出しに対応するメモリ領域が解放されていき、スタックから除かれます。一方、停止条件が記述されていない再帰では、いずれスタック領域の限界を超えてしまう、スタックオーバーフローする、ために、Segmentation faultで終了することになります。

組み合わせ最適化

今日は組み合わせ最適化問題を扱います。組み合わせ最適化問題とは、解が「順序」や「割り当て」などの組み合わせ的な構造をもつ離散的な解として表され、それらの中からある評価基準に基づいて最良の解を探索する問題をいいます。連続最適化がなめらかな山を登る問題だとすると、組み合わせ最適化は無数にある選択肢の中から一番良い組み合わせを探す問題です。

ナップサック問題

ナップサック問題は容量の決まったナップサックと、重さと価値が決まった品物のリストが与えられたときに、ナップサックの限界の重さの制約のもとで品物の価値の合計を最大化する問題です。リスト中の品物がそれぞれ1つだけで、それらを入れるかどうかを決める問題を0-1ナップサック問題といい、以下の式で定式化されます。

\[ {\rm max} \sum_{i=1}^n v_i x_i ~~~ {\rm subject\ to}~ \sum_{i=1}^n w_ix_i \le W, x_i \in (0,1)\]

ナップサック問題は計算の複雑性を議論するときに扱われる非常に有名かつ重要な問題です。

全探索による解法

品物の価値や重さを整数とした整数計画問題として考えると、ナップサック問題には効率的な解法がいくつか知られています。しかし、ここではそれらにトライする前に、再帰関数を用いながら全探索するナイーブなアプローチをとってみます。すなわち個々の品物の「入れる」「入れない」の全ての組み合わせを考えてみます。

item.h
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
#pragma once
#include <stddef.h>  // size_tのため必要

typedef struct item Item;
typedef struct itemset Itemset;


// Itemsetを初期化し、そのポインタを返す関数
// 乱数シードを第2引数にとる
Itemset *init_itemset(size_t number, int seed);

// Itemsetの動的確保された領域を解放する
void free_itemset(Itemset *list);

// Itemsetの内容を標準出力に表示する関数
void print_itemset(const Itemset *list);

// Itemsetの当該indexのアイテムへのポインタを取得
Item *get_item(const Itemset *list, size_t index);

// Itemの個数を取得
size_t get_nitem(const Itemset *list);

// ItemのWeightを取得
double get_itemweight(const Item *item);

// ItemのValueを取得
double get_itemvalue(const Item *item);

// ファイルからItemsetを設定 [未実装] 
Itemset *load_itemset(char *filename);

// Itemsetのパラメータを記録したバイナリファイルを出力する関数 [未実装]
void save_itemset(char *filename);
search.h
1
2
3
4
5
6
7
8
#pragma once
#include "item.h"

// ソルバー関数: 指定された設定でナップサック問題をとく [未完成]
double solve(const Itemset *list, double capacity);

// solve内で呼ぶ再帰的な探索関数: 指定されたindex以降の組み合わせで最適な価値の総和を返す
double search(int index, const Itemset *list, double capacity, unsigned char *flags, double sum_v, double sum_w);
util.h
1
2
3
4
5
#pragma once

// エラー判定付きの読み込み関数
int load_int(const char *argvalue);
double load_double(const char *argvalue);
item.c
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
#include <stdio.h>
#include <stdlib.h>

#include "item.h"

// 以下は構造体定義

// 構造体 Item
// 価値valueと重さweightが格納されている
struct item {
    double value;
    double weight;
};

// 構造体 Itemset
// Itemポインタをmallocする必要あり
struct itemset {
    size_t number;
    Item *item;
};

// 構造体をポインタで確保する作法を確認してみましょう
Itemset *init_itemset(size_t number, int seed) {
    Itemset *list = (Itemset *)malloc(sizeof(Itemset));

    Item *item = (Item *)malloc(sizeof(Item) * number);

    srand(seed);
    for (int i = 0; i < number; i++) {
        item[i].value = 0.25 * (rand() % 200);
        item[i].weight = 0.25 * (rand() % 200 + 1);
    }
    *list = (Itemset){.number = number, .item = item};
    return list;
}

// Itemsetのfree関数
void free_itemset(Itemset *list) {
    free(list->item);
    free(list);
}

// 表示関数
void print_itemset(const Itemset *list) {
    size_t n = list->number;
    const char *format = "v[%d] = %4.1f, w[%d] = %4.1f\n";
    for (int i = 0; i < n; i++) {
        Item *item = get_item(list, i);
        printf(format, i, get_itemvalue(item), i, get_itemweight(item));
    }
    printf("----\n");
}

Item *get_item(const Itemset *list, size_t index) {
    return &(list->item[index]);
}

size_t get_nitem(const Itemset *list) { return list->number; }

double get_itemweight(const Item *item) { return item->weight; }

double get_itemvalue(const Item *item) { return item->value; }
search.c
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
#include "search.h"

#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>

#include "item.h"

// ソルバーは search を index = 0 で呼び出すだけ
double solve(const Itemset *list, double capacity) {
    size_t nitem = get_nitem(list);

    // 品物を入れたかどうかを記録するフラグ配列 =>
    // !!最大の組み合わせが返ってくる訳ではない!!
    unsigned char flags[nitem];
    return search(0, list, capacity, flags, 0.0, 0.0);
}

// 再帰的な探索関数
double search(int index, const Itemset *list, double capacity, unsigned char *flags, double sum_v, double sum_w) {
    size_t max_index = get_nitem(list);
    assert(index >= 0 && sum_v >= 0 && sum_w >= 0);

    // 必ず再帰の停止条件を明記する (最初が望ましい)
    if (index == max_index) {
        const char *format_ok = ", total_value = %5.1f, total_weight = %5.1f\n";
        const char *format_ng = ", total_value = %5.1f, total_weight = %5.1f NG\n";

        bool is_not_exceed = (sum_w <= capacity);  // 重さ評価のtrue/false
        for (int i = 0; i < max_index; i++) {
            printf("%d", flags[i]);
        }
        printf((is_not_exceed) ? format_ok : format_ng, sum_v, sum_w);
        return (is_not_exceed) ? sum_v : 0;
    }

    // 以下は再帰の更新式: 現在のindexの品物を使う or 使わないで分岐し、index
    // をインクリメントして再帰的にsearch()を実行する

    // index番目を使わない場合
    flags[index] = 0;
    const double v0 = search(index + 1, list, capacity, flags, sum_v, sum_w);

    // index番目を使う場合
    flags[index] = 1;
    Item *item = get_item(list, index);
    sum_v += get_itemvalue(item);
    sum_w += get_itemweight(item);
    const double v1 = search(index + 1, list, capacity, flags, sum_v, sum_w);

    // 使った場合の結果と使わなかった場合の結果を比較して返す
    return (v0 > v1) ? v0 : v1;
}
util.c
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
#include "util.h"

#include <errno.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

int load_int(const char *argvalue) {
    long nl;
    char *e;
    errno = 0;  // errno.h で定義されているグローバル変数を一旦初期化
    nl = strtol(argvalue, &e, 10);
    if (errno == ERANGE) {
        fprintf(stderr, "%s: %s\n", argvalue, strerror(errno));
        exit(1);
    }
    if (*e != '\0') {
        fprintf(stderr, "%s: an irregular character '%c' is detected.\n",
                argvalue, *e);
        exit(1);
    }
    return (int)nl;
}

double load_double(const char *argvalue) {
    double ret;
    char *e;
    errno = 0;  // errno.h で定義されているグローバル変数を一旦初期化
    ret = strtod(argvalue, &e);
    if (errno == ERANGE) {
        fprintf(stderr, "%s: %s\n", argvalue, strerror(errno));
        exit(1);
    }
    if (*e != '\0') {
        fprintf(stderr, "%s: an irregular character '%c' is detected.\n",
                argvalue, *e);
        exit(1);
    }
    return ret;
}
main.c
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
#include <assert.h>
#include <stdio.h>
#include <stdlib.h>

#include "item.h"
#include "search.h"
#include "util.h"

// main関数
// プログラム使用例: ./knapsack 10 20
// 10個の品物を設定し、キャパ20 でナップサック問題をとく

int main(int argc, char **argv) {
    /* 引数処理: ユーザ入力が正しくない場合は使い方を標準エラーに表示して終了 */
    if (argc != 3) {
        fprintf(stderr, "usage: %s <the number of items (int)> <max capacity (double)>\n", argv[0]);
        exit(1);
    }

    // 個数の上限はあらかじめ定めておく
    const int max_items = 100;

    const int n = load_int(argv[1]);
    assert(n <= max_items);  // assertで止める
    assert(n > 0);           // 0以下も抑止する

    const double W = load_double(argv[2]);
    assert(W >= 0.0);

    printf("max capacity: W = %.f, # of items: %d\n", W, n);

    // 乱数シードを1にして、初期化 (ここは変更可能)
    int seed = 1;
    Itemset *items = init_itemset(n, seed);
    print_itemset(items);

    // ソルバーで解く
    double total = solve(items, W);

    // 表示する
    printf("----\nbest solution:\n");
    printf("value: %4.1f\n", total);

    free_itemset(items);
    return 0;
}

演習

上記のサンプルプログラムをコンパイルして実行してみましょう。.hファイルをincludeディレクトリに、.cファイルをsrcディレクトリに配置して、binディレクトリ以下にknapsackという名前で実行ファイルを作成してみてください。

# 実行: 第1引数が品物の個数 第2引数が重さの制約
$ ./knapsack 10 70.0

結果として全ての解候補が表示され、そのうち重さ制限に違反するものにNGが表示されているはずです。

全体構造

  • include以下の3つのヘッダファイル
    • item.h
      • 重さと価値をメンバにもつItem構造体のtypedef
      • Item構造体の動的配列とその個数を保持するItemset構造体のtypedef
      • Item/Itemset構造体に関連する関数のプロトタイプ宣言
    • search.h
      • 探索に関するsolve()およびsearch()のプロトタイプ宣言
    • util.h
      • 引数からの整数および実数値読み取り関数（strtolおよびstrtodのラッパー）
  • src以下
    • 上記の実装とmain.c

のようになっています。

main関数

まずmain関数の構造としては

  1. コマンドライン引数で品物の数と容量を決める
  2. 個々の品物の価値、重さはランダムに決める
  3. solve()関数を呼び出す

という流れで実行されています。

solve関数

solve関数はsearch関数に初期値を与えて呼び出す関数となっています。

  • 0番目の品物からスタート
  • 合計価値、合計重量は最初はゼロ
  • 品物の組み合わせを一時記録するflagsはここで定義

search関数

search関数では「i番目以降の品物の「入れない」/「入れる」の全ての組み合わせを考えた場合の、品物の合計価値の最大値」を返します。現在着目している品物を入れるか入れないかの2通りについて異なる引数パターンで再帰的にsearch関数を呼び出すことで全探索を実現しています。このとき本来重さ制限のためにすでに探索する必要のない候補も探索しています。

再帰木

演習

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

  • 枝刈り: 重量制限に違反する解を探索しないよう改良しましょう
    • ヒント: 再帰探索の際、品物をいれると重量オーバする場合にはそちら側の再帰呼び出しを実行しないようにする
  • 最適解の表示: 最後に最適解を表示するように改良しましょう
    • ヒント: searchの返り値を最適値との組み合わせを表す構造体に変更する
  • ビットを用いた全探索: 品物を入れるか入れないかを記録するflagsを、適当なサイズの整数に置き換えてビット演算で実現してみましょう
    • ヒント: 再帰を使わない場合は記録している整数を0からfor文で回せば全探索を実現できます

枝刈り

こたえ

枝刈りは品物のフラグを1にセットした際にナップサックの容量と比較して、すでに重量オーバの場合にsearchを呼び出さないようにして実現します。三項演算子を用いると簡潔に書けます。

search.c
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
#include "search.h"

#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>

#include "item.h"

// ソルバーは search を index = 0 で呼び出すだけ
double solve(const Itemset *list, double capacity) {
    size_t nitem = get_nitem(list);

    // 品物を入れたかどうかを記録するフラグ配列 =>
    // !!最大の組み合わせが返ってくる訳ではない!!
    unsigned char flags[nitem];
    return search(0, list, capacity, flags, 0.0, 0.0);
}

// 再帰的な探索関数
double search(int index, const Itemset *list, double capacity, unsigned char *flags, double sum_v, double sum_w) {
    size_t max_index = get_nitem(list);
    assert(index >= 0 && sum_v >= 0 && sum_w >= 0);

    // 必ず再帰の停止条件を明記する (最初が望ましい)
    if (index == max_index) {
        const char *format_ok = ", total_value = %5.1f, total_weight = %5.1f\n";
        const char *format_ng = ", total_value = %5.1f, total_weight = %5.1f NG\n";

        bool is_not_exceed = (sum_w <= capacity);  // 重さ評価のtrue/false
        for (int i = 0; i < max_index; i++) {
            printf("%d", flags[i]);
        }
        printf((is_not_exceed) ? format_ok : format_ng, sum_v, sum_w);
        return (is_not_exceed) ? sum_v : 0;
    }

    // 以下は再帰の更新式: 現在のindexの品物を使う or 使わないで分岐し、index
    // をインクリメントして再帰的にsearch()を実行する

    // index番目を使わない場合
    flags[index] = 0;
    const double v0 = search(index + 1, list, capacity, flags, sum_v, sum_w);

    // index番目を使う場合
    flags[index] = 1;
    Item *item = get_item(list, index);
    sum_v += get_itemvalue(item);
    sum_w += get_itemweight(item);

    // 重量オーバの場合は呼ばないように変更
    bool is_not_exceed = (sum_w <= capacity);
    const double v1 = (is_not_exceed) ? search(index + 1, list, capacity, flags, sum_v, sum_w) : 0;

    // 使った場合の結果と使わなかった場合の結果を比較して返す
    return (v0 > v1) ? v0 : v1;
}

最適解の表示

こたえ

最適解を表示する場合、再帰の末端でflagsをコピーするとともに、再帰の合流地点（品物を入れた場合と入れない場合の比較をするところ）で、使用しなかった方のflagパターンを適切にfreeします。

search.h
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
#pragma once
#include "item.h"

typedef struct ans{
    double value;
    unsigned char *flags;
} Answer;

// ソルバー関数: 指定された設定でナップサック問題をとく
Answer solve(const Itemset *list, double capacity);

// solve内で呼ぶ再帰的な探索関数: 指定されたindex以降の組み合わせで最適な価値の総和と組み合わせを返す
Answer search(int index, const Itemset *list, double capacity, unsigned char *flags, double sum_v, double sum_w);

search.c
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
#include "search.h"

#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>

#include "item.h"

// ソルバーは search を index = 0 で呼び出すだけ
Answer solve(const Itemset *list, double capacity) {
    size_t nitem = get_nitem(list);

    // 品物を入れたかどうかを記録するフラグ配列 =>
    // !!最大の組み合わせが返ってくる訳ではない!!
    unsigned char flags[nitem];
    return search(0, list, capacity, flags, 0.0, 0.0);
}

// 再帰的な探索関数
Answer search(int index, const Itemset *list, double capacity, unsigned char *flags, double sum_v, double sum_w) {
    size_t max_index = get_nitem(list);
    assert(index >= 0 && sum_v >= 0 && sum_w >= 0);
    const Answer invalid = (Answer){ .value = 0, .flags = NULL};  // 解なしを定義

    // 必ず再帰の停止条件を明記する (最初が望ましい)
    if (index == max_index) {
        const char *format_ok = ", total_value = %5.1f, total_weight = %5.1f\n";
        const char *format_ng = ", total_value = %5.1f, total_weight = %5.1f NG\n";

        bool is_not_exceed = (sum_w <= capacity);  // 重さ評価のtrue/false
        const char *format = (is_not_exceed)? format_ok : format_ng;
        unsigned char *arg = (is_not_exceed)? malloc(sizeof(unsigned char) * max_index) : NULL;

        if (is_not_exceed){
            for (int i = 0; i < max_index; i++) {
                arg[i] = flags[i];
                printf("%d", flags[i]);
            }
        }
        printf(format, sum_v, sum_w);
        return (is_not_exceed) ? (Answer){ .value = sum_v, .flags = arg } : invalid;
    }

    // 以下は再帰の更新式: 現在のindexの品物を使う or 使わないで分岐し、index
    // をインクリメントして再帰的にsearch()を実行する

    // index番目を使わない場合
    flags[index] = 0;
    const Answer a0 = search(index + 1, list, capacity, flags, sum_v, sum_w);

    // index番目を使う場合
    flags[index] = 1;
    Item *item = get_item(list, index);
    sum_v += get_itemvalue(item);
    sum_w += get_itemweight(item);
    bool is_not_exceed = (sum_w <= capacity);
    const Answer a1 = (is_not_exceed) ? search(index + 1, list, capacity, flags, sum_v, sum_w) : invalid;

    // 使った場合の結果と使わなかった場合の結果を比較して返す
    if ( a0.value > a1.value){
        free(a1.flags);
        return a0;
    } else {
        free(a0.flags);
        return a1;
    }
}
main.c
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
#include <assert.h>
#include <stdio.h>
#include <stdlib.h>

#include "item.h"
#include "search.h"
#include "util.h"

// main関数
// プログラム使用例: ./knapsack 10 20
// 10個の品物を設定し、キャパ20 でナップサック問題をとく

int main(int argc, char **argv) {
    /* 引数処理: ユーザ入力が正しくない場合は使い方を標準エラーに表示して終了 */
    if (argc != 3) {
        fprintf(stderr, "usage: %s <the number of items (int)> <max capacity (double)>\n", argv[0]);
        exit(1);
    }

    // 個数の上限はあらかじめ定めておく
    const int max_items = 100;

    const int n = load_int(argv[1]);
    assert(n <= max_items);  // assertで止める
    assert(n > 0);           // 0以下も抑止する

    const double W = load_double(argv[2]);
    assert(W >= 0.0);

    printf("max capacity: W = %.f, # of items: %d\n", W, n);

    // 乱数シードを1にして、初期化 (ここは変更可能)
    int seed = 1;
    Itemset *items = init_itemset(n, seed);
    print_itemset(items);

    // ソルバーで解く
    Answer ans = solve(items, W);

    // 表示する
    printf("----\nbest solution:\n");
    printf("value: %4.1f\n", ans.value);
    for (int i = 0 ; i < get_nitem(items); i++){
        printf("%d", ans.flags[i]);
    }
    printf("\n");

    free_itemset(items);
    return 0;
}

ビットを用いた全探索

こたえ

ビットで表現した場合は、整数をfor文で回すだけで全ての組み合わせパターンが列挙できます。最適解を表すフラグを整数だけで表現できるのがメリットです。

search.h
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
#pragma once
#include "item.h"

typedef struct ans{
    double value;
    unsigned int flags;  // intで管理
} Answer;

// ソルバー関数: 指定された設定でナップサック問題をとく
Answer solve(const Itemset *list, double capacity);

// solve内で呼ぶ再帰的な探索関数: 指定されたindex以降の組み合わせで最適な価値の総和と組み合わせを返す
Answer search(int index, const Itemset *list, double capacity, unsigned int flags, double sum_v, double sum_w);

search.c
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
#include "search.h"

#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>

#include "item.h"

// ソルバーは search を index = 0 で呼び出すだけ
Answer solve(const Itemset *list, double capacity) {
    unsigned int flags = 0; // 32ビット;
    return search(0, list, capacity, flags, 0.0, 0.0);
}

// 全探索関数
Answer search(int index, const Itemset *list, double capacity, unsigned int flags, double sum_v, double sum_w) {

    const char *format_ok = ", total_value = %5.1f, total_weight = %5.1f\n";
    const char *format_ng = ", total_value = %5.1f, total_weight = %5.1f NG\n";
    double max_v = 0;
    unsigned int max_arg = 0;
    size_t nitem = get_nitem(list);

    for (unsigned int i = 0 ; i < 1 << nitem ; i++) {
        double tmp_v = 0;
        double tmp_w = 0;
        for( unsigned int j = 0 ; j < nitem ; j++) {
            // 今回のリストの左からj番目の品物の有無を取り出す
            // ビット演算だけだとどの品物かの情報も含まれているのでboolに代入
            Item *item = get_item(list,j);
            bool b = i & 1 << (nitem-1-j);

            tmp_v += get_itemvalue(item) * b ;
            tmp_w += get_itemweight(item) * b;
            printf("%d", b);
        }
        if (tmp_w <= capacity){
            printf(format_ok,tmp_v,tmp_w);
            if ( tmp_v > max_v){
                max_v = tmp_v;
                max_arg = i;
            }
        } else {
            printf(format_ng,tmp_v,tmp_w);
        }
    }

    return (Answer){ .value = max_v, .flags = max_arg};
}
main.c
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
#include <assert.h>
#include <stdio.h>
#include <stdlib.h>

#include "item.h"
#include "search.h"
#include "util.h"

// main関数
// プログラム使用例: ./knapsack 10 20
// 10個の品物を設定し、キャパ20 でナップサック問題をとく

int main(int argc, char **argv) {
    /* 引数処理: ユーザ入力が正しくない場合は使い方を標準エラーに表示して終了 */
    if (argc != 3) {
        fprintf(stderr, "usage: %s <the number of items (int)> <max capacity (double)>\n", argv[0]);
        exit(1);
    }

    // 個数の上限はあらかじめ定めておく
    const int max_items = 100;

    const int n = load_int(argv[1]);
    assert(n <= max_items);  // assertで止める
    assert(n > 0);           // 0以下も抑止する

    const double W = load_double(argv[2]);
    assert(W >= 0.0);

    printf("max capacity: W = %.f, # of items: %d\n", W, n);

    // 乱数シードを1にして、初期化 (ここは変更可能)
    int seed = 1;
    Itemset *items = init_itemset(n, seed);
    print_itemset(items);

    // ソルバーで解く
    Answer ans = solve(items, W);

    // 表示する
    printf("----\nbest solution:\n");
    printf("value: %4.1f\n", ans.value);
    for(int i = get_nitem(items)-1 ; i >=0 ;i--) {
        printf("%d", (ans.flags & (1 << i))?1:0);
    }
    printf("\n");

    free_itemset(items);
    return 0;
}

より効率的な解法

ナップサック問題は、単純な全探索では品物数に対して指数的な計算量が必要になります。しかし構造が比較的単純であるため、全探索よりも効率の良い厳密解法がいくつか知られています。

代表的な方法の一つが動的計画法です。物の重さとナップサックの容量が整数である場合、「先頭から何個の品物を使って、ある重さ以下で得られる最大価値」という部分問題を再利用することで、計算量を疑似多項式時間まで削減できます。この方法は実装も比較的分かりやすく、容量が小さい場合には実用上きわめて強力です。

もう一つの代表的な方法が分枝限定法です。「これ以上探索しても最適解を更新できない」と判断された部分問題を途中で打ち切ることで、実際に調べる組合せの数を大幅に減らします。

興味のある方はさまざまな解法を調べてみるとよいでしょう。

コラム: ILPソルバ

ナップサック問題は整数線形計画(Integer Linear Programming, ILP)として定式化されます。したがって、一般のILPソルバに解かせることが可能です。CPLEXや、SCIPといったソルバを活用したり、ソースコードを眺めてみると何か発見があるかもしれません。

巡回セールスマン問題

巡回セールスマン問題 (traveling salesman problem, TSP) は都市の集合と移動コストの情報が与えられたときに、全ての都市を巡回する最短経路を求める問題です。ただし同じ都市を2回以上訪れてはいけないという制約があります。

実装例

サンプルプログラムを以下からダウンロードできます。

eeic-software2/tsp

$ git clone https://github.com/eeic-software2/tsp.git

演習

サンプルプログラムを以下のようにコンパイルして実行してみましょう。 

$ cd tsp
# コンパイル（bin以下に gencity, tsp がコンパイルされる）
$ make
# 実行前にターミナルを少し大きくしておく
# ファイルから都市の座標を読み取る
$ bin/tsp data/city10seed3.dat

ここで、makeはプログラムのコンパイルと実行ファイルの生成を行うためのコマンドです。今までのコンパイル手順よりもずっと楽に複数ファイルのコンパイルが済んだと思います。この詳細については次回以降に学んでいきます。

現状は最短経路の計算は実装されておらず、ただ番号順に都市を回るプログラムになっています。gencity.cは乱数のシードと都市数を受け取り、都市座標を表すバイナリファイルを生成するプログラムです。引数なしで実行すると使い方が表示されます。今回のプログラムの検証に適宜用いてください。

  • include/src以下の3つのファイルセット
    • city.h/city.c
      • 町の構造体と関連関数（2都市間の距離計算およびファイルからの都市配置読み込み）
    • map.h/map.c
      • 描画に関する構造体と関連関数
    • solve.h/solve.c
      • TSPを解くための関数
  • src/main.c
    • メイン関数
  • data
    • 都市配置情報のバイナリファイルを置いているディレクトリ

のようになっています。

実際にはsolve関数の中で以下に示す都市の巡回に関する配列を更新し、巡回順が確定した時点で距離を計算するプログラムになっています。

// 都市iの位置を表す構造体: 距離を計算するときに用いる
city[i];
// j番目に訪れる都市の番号: route[0] = 0 とする（必ず都市0から出発）
route[j];
// 都市の巡回順を生成する際に、都市kが巡回済みかを記録する配列
visited[k];
例えば以下に示す配列になっている場合は0->2->4->1->3->0と巡回することになります。

route[0] = 0;
route[1] = 2;
route[2] = 4;
route[3] = 1;
route[4] = 3;

全探索による解法

巡回セールスマン問題の巡回候補を列挙することは、与えられた整数列の順列を求めることに相当します。下図は都市0を出発点として都市数が4の時の巡回パターンを示したものです。

TSP探索木

都市数が4の場合はそれほどの数ではありませんが、都市数が増えると爆発的に探索すべき候補が増えることがわかります。

全探索の前準備として、まずは[0,...,6]の整数列が与えられたときに、その順列を全て生成するプログラムを考えてみましょう。再帰を用いて記述できるはずです。

演習

以下のpermutation.cのpermutation関数を完成させて [0,...,6] の順列を全て生成するプログラムを書いてみましょう

permutation.c
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
#include <stdio.h>
#include <stdlib.h>

void permutation(int *pattern, int *used, size_t number);

int main(int argc, char **argv) {
    size_t number = 7;

    int *pattern = (int *)calloc(number, sizeof(int));
    int *used = (int *)calloc(number, sizeof(int));

    for (int i = 0; i < number; i++) {
        // まだパターンが確定してないことのフラグを立てておく
        pattern[i] = -1;
    }

    permutation(pattern, used, number);

    free(used);
    free(pattern);

    return 0;
}

void permutation(int *pattern, int *used, size_t number) {
    int start = -1;

    // patternをfor文で確認し，未確定の最初のインデックスをstartに入れる

    /* ここにコードをかく */

    // パターンが確定した場合（再帰の終端）
    if (start == -1) {
        printf("[ ");
        for (int i = 0; i < number; i++) {
            printf("%d ", pattern[i]);
        }
        printf("]\n");
        return;
    }

    // パターンが確定してない場合の処理を以下にかく
    // usedとpatternを更新しながらpermutationを呼び出す

    /* ここにコードをかく */

    // 終了
    return;
}

こたえ
permutation.c
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
#include <stdio.h>
#include <stdlib.h>

void permutation(int *pattern, int *used, size_t number);

int main(int argc, char **argv) {
    size_t number = 7;

    int *pattern = (int*)calloc(number, sizeof(int));
    int *used = (int*)calloc(number, sizeof(int));
    for (int i = 0 ; i < number ; i++) {
        // まだパターンが確定してないことのフラグを立てておく
        pattern[i] = -1;
    }

    permutation(pattern, used, number);

    free(used);
    free(pattern);

    return 0;
}

void permutation(int *pattern, int *used, size_t number){
    int start = -1;

    for (int i = 0 ; i < number ; i++) {
        if (pattern[i] == -1) {
            start = i;
            break;
        }
    }
    // パターンが確定した場合（再帰の終端）
    if (start == -1) {
        printf("[ ");
        for (int i = 0 ; i < number ; i++) {
            printf("%d ", pattern[i]);
        }

        printf("]\n");
        return;
    }

    // それ以外の場合は一つパターンを確定させて再帰するのを残りの数だけ繰り返す
    for (int i = 0 ; i < number ; i++) {
        if (!used[i]) {
            pattern[start] = i;
            used[i] = 1;
            permutation(pattern, used, number);
            pattern[start] = -1;
            used[i] = 0;
        }
    }
    return;
}

演習

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

  • solve.c内のsolve関数とsolve.hを修正し、全ての巡回パターンから最短経路を求めるように書き換えましょう。
    • solve関数から、再帰的な探索関数searchを呼び出すようにする
    • searchは最小距離と巡回パターンをセットにした構造体Answerを返すようにする
    • 巡回パターンrouteが重複なく埋まるようにする
    • 各都市がすでに巡回に含まれているかどうかをvisitedで管理する
    • 普通に順列を生成すると、時計回り、反時計回りは別のものとして生成され、同じ距離になる。これを枝刈りするには、例えば特定の2都市の出現順序を固定するなどするとよい（このとき間に別の都市が入るのは構わない）。たとえば、0->1->2->3->4と0->4->3->2->1は同じ経路の循環方向が逆のものになる。1は必ず2の前に通るといった制約を設けることで、この2つの重複を避けることができる。
こたえ

solve.h
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
#pragma once
#include "city.h"

typedef struct ans{
    double dist;
    int *route;
} Answer;

double solve(const City *city, int n, int *route, int *visited);

Answer search(const City *city, int n, int *route, int *visited);
solve.c
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
#include "solve.h"

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include "city.h"
#include "map.h"

double solve(const City *city, int n, int *route, int *visited) {
    route[0] = 0;  // 循環した結果を避けるため、常に0番目からスタート
    visited[0] = 1;

    Answer ans = search(city, n, route, visited);

    memcpy(route, ans.route, sizeof(int) * n);
    free(ans.route);
    return ans.dist;
}

Answer search(const City *city, int n, int *route, int *visited) {
    int start = 0;
    // 訪問した個数を数える
    for (int i = 1; i < n; i++) {
        if (!route[i]) {
            start = i;
            break;
        }
    }
    // 全て訪問したケース（ここが再帰の終端条件）
    if (start == 0) {
        double sum_d = 0;
        for (int i = 0; i < n; i++) {
            const int c0 = route[i];
            const int c1 = route[(i + 1) % n];  // nは0に戻る
            sum_d += distance(get_city(city, c0), get_city(city, c1));
        }
        int *retarg = (int *)malloc(sizeof(int) * n);
        memcpy(retarg, route, sizeof(int) * n);

        return (Answer){.dist = sum_d, .route = retarg};
    }

    // 特定の分岐における最小の巡回経路を調べる
    Answer min = {.dist = 10000000000, .route = NULL};
    for (int i = 1; i < n; i++) {
        // 未訪問なら訪れる
        if (!visited[i]) {
            if (i == 2 && !visited[1]) continue;  // 逆順の巡回経路を抑制

            route[start] = i;
            visited[i] = 1;

            Answer tmp = search(city, n, route, visited);

            // 最小の巡回経路かどうか確認
            if (tmp.dist < min.dist) {
                free(min.route);
                min = tmp;
            } else {
                free(tmp.route);
            }

            route[start] = 0;
            visited[i] = 0;
        }
    }

    return min;
}

山登り法

組み合わせ最適化問題に対する解法は、一般に厳密解法、近似解法、発見的解法（ヒューリスティック）の三つに大別されます。

厳密解法

理論的に最適解を必ず求めることができる解法です。さきほどナップサック問題に対して実施した全探索などはこれにあたります。問題サイズが小さい場合には有効ですが、多くの組み合わせ最適化問題では計算量が問題サイズに応じて急激に増大するため、大規模問題には適用が困難です。

近似解法

最適解そのものではなく、最適解に近い解を効率よく求めることを目的とした解法です。解の品質について理論的な保証が与えられる点が特徴です。

発見的解法（ヒューリスティック）

厳密な最適性保証や近似保証を持たない代わりに、実装が容易で現実的な計算時間で「十分に良い解」を得ることを目的とする解法です。問題固有の構造や直感的な改善規則を利用することが多く、組み合わせ最適化問題に広く用いられています。

ここでは発見的解法としての山登り法を考えましょう。山登り法は以下のプロセスで解候補を修正し、探索していきます。

  1. 初期解を何らかの方法で生成
  2. 解を微修正し、近傍の中でより良い解に更新
  3. 解の改良ができなくなるまで 2 に戻り繰り返す

実装は簡単ですが、最適解が求まる保証はなく、局所最適解にはまる可能性があります。山登り法で最適に近い解を求めたい場合、例えば多数の初期解から設定して得られた複数の局所最適解を最終的に比較するようなアプローチをとることになります。

巡回セールスマンの場合、どのように巡回パターンの近傍を定義するかが問題になります。ここでは都市の巡回パターンのうち2つの都市を入れ替えたものを考え、これを近傍とします。すなわち列として見たときに、ハミング距離が2の状態です。

これを踏まえ以下の方針で探索を行うことを考えます。

  1. 複数の初期解をランダムに生成する
  2. それぞれの初期解について
    1. あるパターンが与えられたときに全ての近傍パターン（2都市の入れ替え）を調べる
    2. 距離が最小だったパターンを採用。2-aに戻る。
    3. 全ての近傍より小さいパターンを見つけたら終了。局所最適解とする
  3. 2で得られた全ての局所最適解を比較し最適なものを解とする。

なお、巡回セールスマンにおいて、この近傍定義における山登り法のことを2-opt法と呼びます。これは今決まっている巡回経路のパス（枝）から\(k\)本を削除して再接続をおこなう\(k\)-opt法において、\(k=2\)の場合です。実用的には、計算量と解の品質のバランスから、2-opt法か3-opt法が用いられます。

レポート

本日はここまでです。宿題として以下をお願いします。

  1. Slackで配付した宿題リンクから、レポートを提出してください。
    • リンクは公開しないでください
    • 締切は今回は次の授業翌日の夜23:59までです（今回の場合、2025/12/26, 23:59）
    • リポジトリを更新するかたちで提出してください
      • リポジトリのトップ階層、またはsrcディレクトリ以下など、わかりやすい場所に作成したコードを配置してください
      • 自動採点的な機能はありません
    • レポートの内容は以下になります

week3_1

巡回セールスマン問題を山登り法で解くプログラムを実装してみましょう。

  • 都市数が20の場合で解が求められることを確認すること
  • デバッグなどで別パターンを利用したい場合、tsp/bin/gencityで都市のパターンは生成できます 
    # 都市数20、シード12345で都市のパターンを生成
$ ./gencity 20 12345 city20seed12345.dat