数値積分のインターフェース

子細合ってヤバい形の積分を解く羽目になったので数値積分をすることにした。 別に数値積分の複雑なアルゴリズムの解説や異常なまでの最適化などを紹介するわけではなく、非常に基本的な台形積分を使い、そこまで最適化も考えずにやる。 だがそれだけだと何の面白みもないので、使い回すことを考えた時に数値積分を実行する関数のインターフェースをどうするか少し考えてみたい。

いくつかのパラメータを取る関数 { \displaystyle
f(x; \alpha, \beta, \gamma)
} がある。 xについて積分した値のパラメータ依存性を見たい。

そのまま書くと以下のようになる。

double f(double x, double a, double b, double c)
{
    /* calculate ... */;
}

double integrate(double xmin, double xmax, std::size_t N,
                 double a, double b, double c)
{
    const double dx = (xmax - xmin) / N;
    double sum = 0.0;
    for(std::size_t i=0; i<N; ++i)
    {
        const double lower = xmin + dx*i;
        const double upper = lower + dx;
        sum += dx * (f(lower, a, b, c) + f(upper, a, b, c)) / 2;
    }
    return sum;
}

int main()
{
    const double xmin = 0.0;
    const double xmax = 10.0;
    const std::size_t Nx = 10000;
    const double minalpha = 0.0;
    const double maxalpha = 10.0;
    const std::size_t Na = 100;

    const double dalpha = (maxalpha - minalpha) / Na;
    for(std::size_t i=0; i<Na; ++i)
    {
        const double alpha = minalpha + dalpha * i;
        const double F =
            integrate(xmin, xmax, Nx, alpha, beta, gamma);
        std::cout << alpha << ' ' << F << '\n';
    }
    return 0;
}

まあ用は足りるだろう。少しずつ汎用にしていくとすると、最も単純な拡張はtemplateだ。より高精度な計算が要求される可能性があるので、long doubleや多倍長で、あるいは精度が必要ないのでfloatで計算したくなるかもしれない。そのための機能がtemplateなので、これは素直な拡張だ。特に難しいTMPもいらない。

template<typename realT>
realT f(realT x, realT a, realT b, realT c)
{
    static_assert(std::is_floating_point<realT>::value, "");
    /* calculate ... */;
}

template<typename realT>
realT integrate(realT xmin, realT xmax, std::size_t N,
                realT a, realT b, realT c)
{
    static_assert(std::is_floating_point<realT>::value, "");

    const realT dx = (xmax - xmin) / N;
    realT sum = 0.0;
    for(std::size_t i=0; i<N; ++i)
    {
        const realT lower = xmin + dx*i;
        const realT upper = lower + dx;
        sum += dx * (f(lower, a, b, c) + f(upper, a, b, c)) / 2;
    }
    return sum;
}

いやまあ、精度のためなら多倍長有理数を受け入れられるようにしたほうがよいのではとかはあるが、とりあえず今回は浮動小数点数だけを見ておく。有理数型などが来た時でも、適切に演算子オーバーロードがされていればstatic_assertのチェックをゆるくするだけで済むだろう。

次に気になるのは、integrate関数がfが取るパラメータを全て受け取ってしまっていることだ。これは被積分関数に関係するパラメータであって、積分に関係するものではない。いざ積分されるとなったとき、関数fが持っていなければならない情報だ。意味的なモヤモヤ以外に、このままではintegrate関数はパラメータをreal型で3つ以下しか受け取らないような関数にしか適用できない。これでは汎用とは言い難いし、ライブラリで提供できるようなインターフェースではない。

C言語的解決策は、パラメータをまとめたstructを作り、void*をパラメータとして受け取った後、パラメータ構造体にキャストして中身を取り出すというものだ。これはGNU Scientific Libraryで用いられている手段である。 GNU Scientific Library – Reference Manual: Providing the function to solve 以下のようなコードを用意することになる。

struct params{double a; double b; double c;};
double f(double x, void* p)
{
    struct params *ps = (struct params *)p;
    double a = ps->a;
    double b = ps->b;
    double c = ps->c;
    /*calculate...*/
}

double integrate(double xmin, double xmax, size_t N,
                double (*f)(double, void*), void*);

これによって、ユーザーが呼び出す関数とパラメータ構造体の関係を適切に保っていれば、積分を実行する関数はそれを単に呼び出すだけで済むので、一般化が可能になるというわけだ。 数値積分を実行するだけのintegrate関数は、渡されたvoid*型が実際にはどんな型であるのかは気にしなくていい。その中身を取り出すのはその型を知っているはずのユーザーが(恐らく適切に)書いた関数であって、integrateではないからだ。

いかにもCという感じがする。ところで我々が書いているのはC++だ。C++では、何かわからない構造体を受け取るならtemplateが使える。void*はお役御免だ。

template<typename realT>
struct params{realT a; realT b; realT c;};

template<typename realT>
realT f(realT x, const params<realT>& ps)
{
    realT a = ps.a;
    realT b = ps.b;
    realT c = ps.c;
    /*calculate...*/
}

template<typename realT, typename paramT>
realT integrate(realT xmin, realT xmax, std::size_t N,
                realT(*f)(realT, const paramT&), paramT ps)
{
//...
    f(x, ps);
//...
}

少しC++っぽくなった。void*がなくなったからだ。だが関数ポインタを使っている辺りがまだ微妙な感じがする。 そもそも、積分するときに必要なのは関数であって、関数とそのパラメータではない。パラメータを受け取らない関数もある。 integrateの中ではf(x)のように呼べた方がよい。そして、C++には関数に値を持たせる方法がいくつかある。一つがstd::bindだ。

template<typename realT>
realT f(realT x, realT a, realT b, realT c);

int main()
{
    using namespace std::placeholders;
    const auto f_ =
        std::bind(f<double>, _1, /*a=*/2.0, /*b=*/3.0, /*c=*/4.0);
    std::cout << f_(x) << std::endl;
    return 0;
}

ただしこのstd::bindの返り値はunspecifiedなので、変数で受けるときにはautoを使うかstd::function<double(double)>のようにするしかない。 ちなみにstd::functionは毎回空になっていないかチェックしたはずで、オーバーヘッドが入る点は微妙だ。実行時に入れる場合はチェックは必須だが、コンパイル時に型が決まるならautoにしたほうが速くなるだろう。

この方法を取ると、integrateは以下のようになる。

template<typename realT, typename Function>
realT integrate(realT xmin, realT xmax, std::size_t N,
                Function&& func)
{
// ...
    sum += func(x);
// ...
}

int main()
{
    using namespace std::placeholders;
    const double F = integrate(xmin, xmax, Nx,
                               std::bind(f<double>, _1, a, b, c));
}

integrateの宣言は明快になった。十分に汎化されたので、std::bind以外の「関数のようなもの」を呼び出せるようになった。関数そのものはもちろん、関数オブジェクトもだ。

template<typename realT>
realT f(realT x){return x * x;}

template<typename realT>
struct functor
{
    realT operator()(const realT x){return x * a;}
    realT a;
};

int main()
{
    {
        const double F =
            integrate(xmin, xmax, Nx, f<double>);
    }
    {
        const double F =
            integrate(xmin, xmax, Nx, functor<double>{1.0});
    }
}

ここまでで十分だろう。パラメータ設定済みの関数はそのまま渡し、実行時にパラメータを渡す関数はstd::bindでくっつけるか関数オブジェクトで渡す。これら3つで呼び出しのインターフェースが揃っていることが重要だ。templateはインターフェースが揃っているオブジェクトは大抵受け取れる。ダック・タイピングというわけだ。

若干蛇足気味だが、関数オブジェクトが使えるようになったのでまだ少し幅が広がる。N次元空間で積分するような関数を考えよう。 頭を使わない実装だと効率が大変なことになりそうだが、そちらは気にしないことにする。インターフェースだけ考えて、中では頭の良いアルゴリズムが使われていることにしておく。

ここでintegrateには関数オブジェクトを渡すことにして、その関数オブジェクトではいくつかの値を定義しておくことにする。

template<typename realT>
struct functor
{
    typedef realT variable_type;
    typedef realT result_type;
    constexpr static std::size_t dim = 1;

    result_type operator()(const variable_type x){return x * a;}
    realT a;
};

そして、これまで別々に渡していた積分区間と数値積分における幅をまとめた型を作る。

template<typename realT, std::size_t N>
struct range
{
    constexpr static std::size_t dim = N;
    typedef std::array<realT, N> variable_type;
    std::array<realT, N> vmin;
    std::array<realT, N> vmax;
    std::array<std::size_t, N> steps;
};
template<typename realT>
struct range<realT, 1>
{
    constexpr static std::size_t dim = N;
    typedef realT variable_type;
    realT vmin;
    realT vmax;
    std::size_t steps;
};

こうすれば、

template<typename Range, typename Function>
typename Function::result_type
integrate(Range&& rg, Function&& func)
{
    static_assert(Range::dim == Function::dim, "");
    static_assert(
        std::is_same<typename Range::variable_type,
            typename Function::variable_type>::value, "");
    /* ... */
}

int main()
{
    const auto F = integrate(
        range<double, 2>{{{0.0, 0.0}}, {{1.0, 2.0}}, {{100, 200}}},
        func{1.0, 2.0});
}

のようにして、「積分区間」と「パラメータ設定済み関数」を渡すことで積分が可能になる。