Cに書き直して満足して終了

Haskell vs F#
↑のエントリを見て、↓のエントリを書いた。今日はその続き。
HaskellとF#の解読の練習とちょっとした疑問でRubyに移植してみた

.NETerやHaskllerが高速化のために高い技術力を投下してる裏で、ワタクシおなじみの言語で書き換えたらどうなるかを一人でやってます。今日はコードをCに直してみました。やっぱ速い。

結果

$ gcc --version
i686-apple-darwin11-llvm-gcc-4.2 (GCC) 4.2.1 (Based on Apple Inc. build 5658) (LLVM build 2335.15.00)
Copyright (C) 2007 Free Software Foundation, Inc.
This is free software; see the source for copying conditions.  There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

$ gcc -O2 hoge.c
$ time ./a.out 
1000.000000
./a.out  0.06s user 0.00s system 82% cpu 0.070 total

はっ、速い！！同じマシン上でHaskellの0.13倍!!(昨日のエントリ参照)

移植遊びの感想

前回エントリと今回エントリを通したベンチマークでは、プログラミング言語の表現力の豊かさと、内部的に計算器へ値を渡すショートカットがいかに両方優れているかを競う感じになってて面白かった。そういう意味でF#もHaskellも凄い優れてると思った。

自分の中での収穫はF#が意外に良いって事が分かったこと。性能について聞いたことなかったけど、意外にもかなり良いみたいだし、記述も簡潔なんだなぁとかなり見直した。あまり、注目してなかったけど機会があれば是非使ってみたい言語だと思った。

それにしてもC速いといってもHaskellのコード量の3〜4倍くらいになってしまった・・・。こんだけ書いてやっとこさなので、徒労感が激しい。「人月の神話」によると、プログラミング言語に限らず、一人のエンジニアが一日に記述できるコード量は一定らしい。やっぱり計算機に優しいとか前例があるとかだけの理由で古くて表現力が乏しい言語を使い続けるのは問題がある。今回だとinduceBackward()内だけCで書くとかがいい。道具は選ぼう。やっぱCって工数的にはとても贅沢だよねぇと、そうとう月並みな感想。

使ったコード

無駄にmalloc()してるところあるので、減らせばもう少し速くなるかな?

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

#define BRANCH_SIZE 3

#define ITERATION   1000

typedef struct Branch_t
{
    int    k;
    double p;
} Branch;

typedef struct Node_t
{
    double df;
    Branch branches[BRANCH_SIZE];
} Node;


double* createValues(int size)
{
    double* result;
   
    if (!(result = malloc(size * sizeof (double)))) 
    {
        fprintf(stderr, "alloc error...\n");
        exit(1);
    }
    return result;
}

void releaseValues(double** values)
{
    if (*values)
        free(*values);
    *values = NULL;
}

void induceBackwardVerySlow(double* result,
                            Node*   nodes,  int nodes_length, 
                            double* values, int values_length)
{
    int i, j, n;
    
    n = values_length / 2;

    for (i = 0; i < nodes_length; i++)
    {
        Node*  node;
        double sum;

        node = nodes + i;
        sum  = 0.0;

        for (j = 0; j < BRANCH_SIZE; j++) {
            Branch* branch;

            branch = node->branches + j;
            sum += branch->p * values[n + branch->k] * node->df;
        }
        result[i] = sum;
    }
}

#define TEST_TREE_LENGTH 100
double* lastValues(int init) 
{ 
    double* values;
    int i;

    values = createValues(TEST_TREE_LENGTH * 2 + 1);

    for (i = 0; i < TEST_TREE_LENGTH * 2 + 1; i++)
        values[i] = (double) init;

    return values;
}

Node** testTree() 
{
    Node** tree;
    int i, j, m;
    
    if (!(tree = malloc(TEST_TREE_LENGTH * sizeof (Node*))))
    {
        fprintf(stderr, "alloc error...\n");
        exit(1);
    }

    for (i = TEST_TREE_LENGTH - 1, m = 0; i >= 0; i--, m++) 
    {
        int l;

        if (!(tree[m] = malloc((i * 2 + 1) * sizeof (Node))))
        {
            fprintf(stderr, "alloc error...\n");
            exit(1);
        }

        for (j = -i, l = 0; j <= i; j++, l++) 
        {
            tree[m][l].df = 1.0;
            tree[m][l].branches[0].k = j - 1;
            tree[m][l].branches[0].p = 1.0 / 6.0;
            tree[m][l].branches[1].k = j;
            tree[m][l].branches[1].p = 2.0 / 3.0;
            tree[m][l].branches[2].k = j + 1;
            tree[m][l].branches[2].p = 1.0 / 6.0;
        }
    }
    return tree;
}

double value(Node** tree, int tree_length, int i) 
{
    double  result;
    double* last_values;
    int j;

    last_values = lastValues(i);

    for (j = 0; j < tree_length; j++) 
    {
        Node*   nodes;
        int     nodes_length;
        double* result_values;
        int     values_length;
        
        nodes = tree[j];

        nodes_length  = 2 * tree_length - j * 2 - 1;
        values_length = nodes_length + 2;

        result_values = createValues(nodes_length);

        induceBackwardVerySlow(result_values, 
                               nodes, nodes_length, 
                               last_values, values_length);

        releaseValues(&last_values);
        last_values = result_values;
    }
    result = last_values[0];
    releaseValues(&last_values);
    return result;
}

int main(int argc, char* argv[]) 
{
    Node** tree;
    double max;
    int i;

    max = DBL_MIN;
    tree = testTree();

    for (i = 1; i <= ITERATION; i++)
    {
        double ret;

        ret = value(tree, TEST_TREE_LENGTH, i);

        if (ret > max)
            max = ret;
    }
    printf("%f\n", max);

    return 0;
}