行き詰った

去年の暮れに、今年は Elixir をやる!と言っておきながらいまだにやってない。それどころか、C で MD5 ハッシュ値を計算するプログラムを作っていて行き詰ってしまった。
きっかけは、MD5 の計算アルゴリズムが意外にシンプルだとを知ったこと。具体的には次のページを見たことだ。

 cf. MD5 メッセージダイジェストアルゴリズム
 cf. MD5のアルゴリズム – 数学自由研究

確かに、考えていたよりは単純な計算で、それほど引っかかるようなところもなく MD5 ハッシュ値らしい文字列を出力するところまでは実装できた(もちろん、久しぶりの C なので些細な文法ミスで引っかかるところはあったけれど)。
が、結果が合わない。試しに sample.zip ファイルのハッシュ値を、いつも使っている md5make というソフトで求めてみる(GUI のソフトだけど、<ファイル名>.md5 というファイルにハッシュ値を書き出してくれる)と:

^o^ > cat sample.zip.md5
d7db18d8239cca197c66175a8cfb0e85

これに対して、オレの書いた md5hash プログラムの出力は:

^o^ > md5hash sample.zip
ece936c4aefc960c882f49b2015aac27

となる。

コードはこのエントリに最後に載せるけど、md5hash にはデバッグのためのオプション -d-t をつけてみた。
-d オプションは、対象のファイルから読み込んだブロック(MD5 アルゴリズムは 32 ビットを 1 ワードとして、16 ワードずつ処理するのでこれをブロックと呼んでいる)を表示するオプション。データ長を 16 ワードに合わせるためのパディングなども含めて出力する。
-t オプションは、計算に使う定数テーブルを出力するオプション。
どちらのオプションの出力も、見た限りでは間違っていないように見える。
となると、あとは実際の計算をしている部分に間違いがありそうだということになるんだけど、個々の関数やその呼び出し引数をチェックしてみても、間違ったところは見当たらない。
というわけで、行き詰っている。どうしよう。
とにかく、書いたコードを載せておく。しばらくして何か思いつくかもしれない。

#include
#include
#include

#include "md5.h"

void Usage(void);

int main(int argc, char *argv[])
{
    FILE *fp;
    WORD block[16];
    int p;
    int count = 0;
    int i;
    char hex[33];

    CONTEXT *context = NULL;
    unsigned long t[64];
    char md5hash[33];

    int opt, opt_dump = 0, opt_table = 0, option_index;
    struct option long_options[] = {
        {"dump", no_argument, NULL, 'd'},
        {"table", no_argument, NULL, 't'},
        {"help", no_argument, NULL, 'h'},
        {0, 0, 0, 0}
    };

    while ((opt = getopt_long(argc, argv, "dth", long_options, &option_index)) != -1) {
        switch(opt) {
        case 'd':
            opt_dump = 1;
            break;
        case 't':
            opt_table = 1;
            break;
        case 'h':
            Usage();
            exit(0);
        case '?':
            printf("Unkown option -%c\n", optopt);
            Usage();
            exit(1);
        }
    }

    if (opt_dump) {
        fp = fopen(argv[optind], "rb");
        while (1) {
            p = ReadBlock(fp, block);
            if (p == 64) {
                printf("block readed: %d\n", ++count);
                for (i = 0; i < 16; i++) {
                    HexEncode(block[i], hex);
                    printf(" %s\n", hex);
                 }
            } else {
                printf("rest %d bytes\n", p);
                for (i = 0; i < 16; i++) {
                    HexEncode(block[i], hex);
                    printf(" %s\n", hex); } break;
                }
            }
            fclose(fp);
        } else if (opt_table) {
            InitTable(t);
            for (i = 1; i<= 64; i++) {
                HexEncode(t[i - 1], hex);
                printf("%2d: 0x%s\n", i, hex);
            }
        } else {
            context = (CONTEXT *) malloc(sizeof(CONTEXT));
            InitContext(context);
            InitTable(t);
            fp = fopen(argv[optind], "rb");
            MD5Calc(fp, context, t, md5hash);
            fclose(fp);
            printf("%s\n", md5hash);
            free(context);
        }
    return 0;
}

void Usage(void)
{
    printf("Usage: md5hash [-d | -t | -h] \n");
    return;
}
typedef unsigned long WORD; // 32 bits = 4 bytes

typedef struct tagContext {
    WORD A;
    WORD B;
    WORD C;
    WORD D;
} CONTEXT;

void InitContext(CONTEXT *context);
void InitTable(unsigned long *t);
void HexMD5(CONTEXT *context, char *hex);
int ReadBlock(FILE *fp, WORD *block);
int CalcPaddingSize(int rest);
void CharToBlock(unsigned char *b, WORD *block, int s, unsigned long long filesize);
WORD ULongToWORD(unsigned long n);
void HexEncode(WORD n, char *hex);
void HexEncode2(WORD n, char *hex);

WORD F(WORD x, WORD y, WORD z);
WORD G(WORD x, WORD y, WORD z);
WORD H(WORD x, WORD y, WORD z);
WORD I(WORD x, WORD y, WORD z);
WORD BitRotate(WORD x, int n);
WORD FF(WORD a, WORD b, WORD c, WORD d, WORD x, int s, unsigned long t);
WORD GG(WORD a, WORD b, WORD c, WORD d, WORD x, int s, unsigned long t);
WORD HH(WORD a, WORD b, WORD c, WORD d, WORD x, int s, unsigned long t);
WORD II(WORD a, WORD b, WORD c, WORD d, WORD x, int s, unsigned long t);
void MD5Calc(FILE *fp, CONTEXT *context, unsigned long *t, char *hex);
void RoundBlock(CONTEXT *context, WORD *block, unsigned long *t);
/*
 * MD5 hash calculation library
 *
 * cf. http://www.ipa.go.jp/security/rfc/RFC1321JA.html
 * cf. http://azisava.sakura.ne.jp/math/md5.html
 *
 */
#include
#include
#include
#include

#include "md5.h"

void InitContext(CONTEXT *context)
{
    context->A = 0x67452301;
    context->B = 0xEFCDAB89;
    context->C = 0x98BADCFE;
    context->D = 0x10325476;

    return;
}

void HexMD5(CONTEXT *context, char *hex)
{
    char hex_a[9], hex_b[9], hex_c[9], hex_d[9];
    int i;

    HexEncode2(context->A, hex_a);
    HexEncode2(context->B, hex_b);
    HexEncode2(context->C, hex_c);
    HexEncode2(context->D, hex_d);

    for (i = 0; i < 8; i++) {
        hex[i] = hex_a[i];
    }
    for (i = 8; i < 16; i++) {
        hex[i] = hex_b[i - 8];
    }
    for (i = 16; i < 24; i++) {
        hex[i] = hex_c[i - 16];
    }
    for (i = 24; i < 32; i++) {
        hex[i] = hex_d[i - 24];
    }
    hex[32] = '\0';
    return;
}

 int ReadBlock(FILE *fp, WORD *block)
{
    int r; unsigned char b[64];
    unsigned long long filesize;

    if ((r = fread(b, 1, 64, fp)) < 64) {
        filesize = ftell(fp);
        CharToBlock(b, block, r, filesize * 8);
    } else {
        CharToBlock(b, block, r, 0);
    }

    return r;
 }

int CalcPaddingSize(int rest)
{
    int p;

    if (rest < 56) {
        p = 56 - rest;
    } else if (rest == 56) {
        p = 64;
    } else {
        p = 64 - rest + 56;
    }

     return p;
}

void CharToBlock(unsigned char *b, WORD *block, int s, unsigned long long filesize)
{
    int i, j, k = 0;
    int p;
    unsigned long bk[16];
    unsigned char c[64];
    unsigned long fs_u, fs_l;

    if (s == 64) {
        for (i = 0; i < 64; i += 4) {
            bk[k] = 0;
            for (j = 0; j < 4; j++) {
                bk[k] <<= 8;
                bk[k] = bk[k] + b[i + j];
             }
            k++;
        }
    } else {
        for (i = 0; i < s; i++) {
            c[i] = b[i];
        }
        p = CalcPaddingSize(s);
        c[i] = 0x80;
        i++;
        for (; i < s + p; i++) {
            c[i] = 0;
        }
        fs_l = (unsigned long) filesize;
        filesize >>= 32;
        fs_u = (unsigned long) filesize;
        for (i = 63; i >= 60; i--) {
            c[i] = fs_u;
            fs_u >>= 8;
        }
        for (i = 59; i >= 56; i--) {
            c[i] = fs_l;
            fs_l >>= 8;
        }

        for (i = 0; i < 64; i += 4) {
            bk[k] = 0;
            for (j = 0; j < 4; j++) {
                bk[k] <<= 8;
                bk[k] = bk[k] + c[i + j];
            }
            k++;
        }
    }
     for (i = 0; i < 16; i++) {
        block[i] = ULongToWORD(bk[i]);
    }

     return;
}

WORD ULongToWORD(unsigned long n)
{
    unsigned char b;
    WORD w = 0;
    int i;

    for (i = 0; i < 4; i++) {
        b = (unsigned char) n;
        w = (w << 8) + b;
        if (i < 3) {
            n = n >> 8;
        }
    }

    return w;
}

void HexEncode(WORD n, char *hex)
{
    int r, i, j;
    char t[16] = {'0', '1', '2', '3', '4', '5', '6', '7','8', '9', 'a', 'b', 'c', 'd', 'e', 'f'};

    i = 0;
    hex[0] = '\0';
    while (1) {
        r = n % 16;
        for (j = i; j >= 0; j--) {
             hex[j + 1] = hex[j];
        }
        hex[0] = t[r];
        i++;
        n = n / 16;
        if (n == 0) {
            break;
        }
    }

    if (i < 8) {
        for (; i < 8; i++) {
            for (j = i; j >= 0; j--) {
                hex[j + 1] = hex[j];
            }
            hex[0] = '0';
        }
    }

    return;
}

void HexEncode2(WORD n, char *hex)
{
    char h[9];

    HexEncode(n, h);
    hex[0] = h[6];
    hex[1] = h[7];
    hex[2] = h[4];
    hex[3] = h[5];
    hex[4] = h[2];
    hex[5] = h[3];
    hex[6] = h[0];
    hex[7] = h[1];
    hex[8] = '\0';

    return;
}

/*
 * Initialize constants table
 */
void InitTable(unsigned long *t)
{
    double p;
    int n;

    p = pow(2.0, 32.0);

    for (n = 1; n <= 64; n++) {
        t[n - 1] = (unsigned long) (p * fabs(sin(n)));
    }

    return;
}

 /*
 * MD5 calcrate function
 */
void MD5Calc(FILE *fp, CONTEXT *context, unsigned long *t, char *hex)
{
    WORD block[16];
    int p;

    while (1) {
        p = ReadBlock(fp, block);
        RoundBlock(context, block, t);
        if (p < 64) {
            break;
        }
    }
    HexMD5(context, hex);

    return;
}

void RoundBlock(CONTEXT *context, WORD *X, unsigned long *T)
{
    WORD a = context->A;
    WORD b = context->B;
    WORD c = context->C;
    WORD d = context->D;

// round 1
    a = FF(a, b, c, d, X[0], 7, T[0]);
    d = FF(d, a, b, c, X[1], 12, T[1]);
    c = FF(c, d, a, b, X[2], 17, T[2]);
    b = FF(b, c, d, a, X[3], 22, T[3]);
    a = FF(a, b, c, d, X[4], 7, T[4]);
    d = FF(d, a, b, c, X[5], 12, T[5]);
    c = FF(c, d, a, b, X[6], 17, T[6]);
    b = FF(b, c, d, a, X[7], 22, T[7]);
    a = FF(a, b, c, d, X[8], 7, T[8]);
    d = FF(d, a, b, c, X[9], 12, T[9]);
    c = FF(c, d, a, b, X[10], 17, T[10]);
    b = FF(b, c, d, a, X[11], 22, T[11]);
    a = FF(a, b, c, d, X[12], 7, T[12]);
    d = FF(d, a, b, c, X[13], 12, T[13]);
    c = FF(c, d, a, b, X[14], 17, T[14]);
    b = FF(b, c, d, a, X[15], 22, T[15]);

// round 2
    a = GG(a, b, c, d, X[1], 5, T[16]);
    d = GG(d, a, b, c, X[6], 9, T[17]);
    c = GG(c, d, a, b, X[11], 14, T[18]);
    b = GG(b, c, d, a, X[0], 20, T[19]);
    a = GG(a, b, c, d, X[5], 5, T[20]);
    d = GG(d, a, b, c, X[10], 9, T[21]);
    c = GG(c, d, a, b, X[15], 14, T[22]);
    b = GG(b, c, d, a, X[4], 20, T[23]);
    a = GG(a, b, c, d, X[9], 5, T[24]);
    d = GG(d, a, b, c, X[14], 9, T[25]);
    c = GG(c, d, a, b, X[3], 14, T[26]);
    b = GG(b, c, d, a, X[8], 20, T[27]);
    a = GG(a, b, c, d, X[13], 5, T[28]);
    d = GG(d, a, b, c, X[2], 9, T[29]);
    c = GG(c, d, a, b, X[7], 14, T[30]);
    b = GG(b, c, d, a, X[12], 20, T[31]);

// round 3
    a = HH(a, b, c, d, X[5], 4, T[32]);
    d = HH(d, a, b, c, X[8], 11, T[33]);
    c = HH(c, d, a, b, X[11], 16, T[34]);
    b = HH(b, c, d, a, X[14], 23, T[35]);
    a = HH(a, b, c, d, X[1], 4, T[36]);
    d = HH(d, a, b, c, X[4], 11, T[37]);
    c = HH(c, d, a, b, X[7], 16, T[38]);
    b = HH(b, c, d, a, X[10], 23, T[39]);
    a = HH(a, b, c, d, X[13], 4, T[40]);
    d = HH(d, a, b, c, X[0], 11, T[41]);
    c = HH(c, d, a, b, X[3], 16, T[42]);
    b = HH(b, c, d, a, X[6], 23, T[43]);
    a = HH(a, b, c, d, X[9], 4, T[44]);
    d = HH(d, a, b, c, X[12], 11, T[45]);
    c = HH(c, d, a, b, X[15], 16, T[46]);
    b = HH(b, c, d, a, X[2], 23, T[47]);

// round 4
    a = II(a, b, c, d, X[0], 6, T[48]);
    d = II(d, a, b, c, X[7], 10, T[49]);
    c = II(c, d, a, b, X[14], 15, T[50]);
    b = II(b, c, d, a, X[5], 21, T[51]);
    a = II(a, b, c, d, X[12], 6, T[52]);
    d = II(d, a, b, c, X[3], 10, T[53]);
    c = II(c, d, a, b, X[10], 15, T[54]);
    b = II(b, c, d, a, X[1], 21, T[55]);
    a = II(a, b, c, d, X[8], 6, T[56]);
    d = II(d, a, b, c, X[15], 10, T[57]);
    c = II(c, d, a, b, X[6], 15, T[58]);
    b = II(b, c, d, a, X[13], 21, T[59]);
    a = II(a, b, c, d, X[4], 6, T[60]);
    d = II(d, a, b, c, X[11], 10, T[61]);
    c = II(c, d, a, b, X[2], 15, T[62]);
    b = II(b, c, d, a, X[9], 21, T[63]);

    context->A += a;
    context->B += b;
    context->C += c;
    context->D += d;

    return;
}

/*
 * MD5 basic functions
 */
WORD F(WORD x, WORD y, WORD z)
{
    return (x & y) | (~x & z);
}

WORD G(WORD x, WORD y, WORD z)
{
    return (x & z) | (y & ~z);
}

WORD H(WORD x, WORD y, WORD z)
{
    return x ^ y ^ z;
}

WORD I(WORD x, WORD y, WORD z)
{
    return y ^ (x | ~z);
}

/*
 * Bit rotation
 */
WORD BitRotate(WORD x, int n)
{
    return x << n | x >> (32 - n);
}

/*
 * MD5 round functions
 */
WORD FF(WORD a, WORD b, WORD c, WORD d, WORD x, int s, unsigned long t)
{
    WORD aa;
    aa = BitRotate(a + F(b, c, d) + x + t, s);
    return b + aa;
}

WORD GG(WORD a, WORD b, WORD c, WORD d, WORD x, int s, unsigned long t)
{
    WORD aa;
    aa = BitRotate(a + G(b, c, d) + x + t, s);
    return b + aa;
}

WORD HH(WORD a, WORD b, WORD c, WORD d, WORD x, int s, unsigned long t)
{
    WORD aa;
    aa = BitRotate(a + H(b, c, d) + x + t, s);
    return b + aa;
}

WORD II(WORD a, WORD b, WORD c, WORD d, WORD x, int s, unsigned long t)
{
    WORD aa;
    aa = BitRotate(a + I(b, c, d) + x + t, s);
    return b + aa;
}