编译时数值积分

• 下载 compile_time_integration_v.2.0.zip - 6.6 KB

简要前言

这是我第二次尝试实现基本的编译时数值积分方法。在第一篇文章的版本中，积分是在 main() 函数之前完成的，但实际上并非在编译时完成。我非常感谢迄今为止收到的所有建设性的意见和反馈。我希望这次一切都会顺利…… 

引言

自 C++11 以来，我们可以在编译时执行相当复杂的计算。在本篇文章中，我将展示如何计算任意实函数在任意边界下的面积。我将使用标准且广泛的算法，例如梯形法则和辛普森法则。

编译时计算意味着，在程序加载到内存时，远在进入 main 函数之前，操作的结果就已经知道了。

main 函数的魔力

看看这段代码（仅自 C++11 起有效）

#include <cmath>
#include <iostream>
#include "simpsons_rule.h"


constexpr double PI=4.0*atan(1.0);    // PI number


struct my_function
{
    constexpr static double f(double x) {return cos(x);}
};


struct left_boundary
{
    constexpr static double value=0.0;
};


struct right_boundary
{
    constexpr static double value=PI/2.0;
};


int main(int argn, char *argc[])
{
    constexpr double partition_size=100;

    // Application of the Simpson's rule.
    // my_function::f is integrated from left_boundary::value to
    // right_boundary::value with a partition size of partition_size.
    // A recursive technique is used.

    constexpr double simpsons_rule_result=
        simpsons_rule<
            my_function,
            left_boundary,
            right_boundary,
            partition_size
            >::result;

    std::cout << "Simpson's rule: " << simpsons_rule_result << std::endl;

    static_assert(fabs(simpsons_rule_result-1.0)<0.001, "Failure"); // [1]


    return 0;
} 

上面的程序只是显示了属于 simpsons_rule<> 类的静态常量表达式 (constexpr) result 的值。这个值在编译时计算，并在访问 main() 之前存储在内存中。result 存储的是函数 my_function::f 在 0.0 (left_boundary::value) 和 PI/2.0 (right_boundary::value) 之间的积分。精确结果是 1。显示的结果也是 1。因此，在这种情况下，近似值几乎是完美的（考虑舍入误差）。

如果您现在更改 my_function::f、left_boundary::value 或 right_boundary::value 并重新编译，您将获得 result 的不同值。

所以，这就是问题的关键和魔力所在：运行时不进行任何计算；唯一执行的任务就是显示结果。要证明这一点，只需像这样更改表达式 [1]

static_assert(fabs(simpsons_rule_result-1.0)>0.001, "Failure"); // [1] Changed "<" into ">".

如果您现在尝试重新编译，编译将失败，因为积分误差实际上小于 0.001。

现在是时候看看内部机制了，但与其深入研究梯形法则或辛普森法则的实现，不如展示一个更简单（但类似）的案例：阶乘。

实现阶乘 (1)

要计算整数 n 的阶乘（表示为 n!），只需应用以下公式

下面的代码实现了编译时阶乘计算（仅自 C++11 起有效）

#ifndef FACTORIAL_H
#define FACTORIAL_H


#include <cstddef>


template <size_t N>     // [1]
struct factorial
{
    constexpr static size_t result=N*factorial<N-1>::result;;
};

template <>
struct factorial<1>     // [2]
{
    constexpr static size_t result=1;
};


#endif

上面的代码声明了一个名为 factorial<> 的类模板 [1] 以及一个特化版本 factorial<1> [2]。在 [1] 中，为 factorial<> 定义了静态 constexpr result。在 [2] 中，为 factorial<1> 定义了静态 constexpr result。

现在假设您在 main 函数中写下这句话

std::cout << factorial<5>::result << std::endl;

如您在 [1] 中所见，编译器将 factorial<5>::result 解析为

factorial<5>::result = 5 * factorial<4>::result

同样在 [1] 中，factorial<4>::result 被解析为

factorial<4>::result = 4 * factorial<3>::result

沿着这个思路递归下去，我们得到

factorial<5>::result = 5 * 4 * 3 * 2 * factorial<1>::result

factorial<1>::result 在 [2] 中被解析为 1。我们因此得到

factorial<5>::result = 5 * 4 * 3 * 2 * 1

编译器最终执行了上面的乘法，factorial<5>::result 在编译时被解析为 120。

实现阶乘 (2)

下面的代码是一个不太为人所知的替代方案（仅自 C++11 起有效）

#ifndef FACTORIAL_H
#define FACTORIAL_H


#include <cstddef>


// With the following definition, factorial<N> is recursively
// resolved as factorial<N, 1, 1, 2, 3, ..., N>.
template <size_t N, size_t M=N+1, size_t... Indices>
struct factorial : public factorial<N, M-1, M-1, Indices...>    // [1]
{};

// The sequence above is split into the First index {1} plus
// the Rest...: {2, 3, ..., N}.
template <size_t N, size_t First, size_t... Rest>
struct factorial<N, 1, First, Rest...>                          // [2]
{
    constexpr static size_t result=
        First*factorial<N, 1, Rest...>::result;
        // (Second factor evaluated at [2])
};

// And finally, this specialization matches the case in which
// all indices have been consumed.
template <size_t N>
struct factorial<N, 1>                                          // [3]
{
    constexpr static size_t result=1;
};


#endif

在这种情况下，编译器如何解析 factorial<5>::result？如果您看一下 [1]，您会发现 factorial<5> 可以解析为 factorial<5, 5, 5>，其中 Indices… 是一个空包。在不离开 [1] 的情况下，您可以沿着链条进行

factorial<5> -> factorial<5, 5, 5>, with Indices... empty.
factorial<5, 5, 5> -> factorial<5, 4, 4, 5>, with Indices... = {5}.
factorial<5, 4, 4, 5> -> factorial<5, 3, 3, 4, 5>, with Indices... = {4, 5}.
factorial<5, 3, 3, 4, 5> -> factorial<5, 2, 2, 3, 4, 5>, with Indices... = {3, 4, 5}.
factorial<5, 2, 2, 3, 4, 5> -> factorial<5, 1, 1, 2, 3, 4, 5>, with Indices... = {2, 3, 4, 5}. 

一旦编译器获得 factorial<5, 1, 1, 2, 3, 4, 5>，它就可以通过 [2] 中的特化版本进行解析。它变成了 factorial<5, 1, First, Rest…>，其中 First=1，Rest… 是包：{ 2, 3, 4, 5}。

#ifndef FIBONACCI_H
#define FIBONACCI_H


#include <cstddef>


template <
    size_t SequenceSize,            // 1, 2, 3,...
    size_t CurrentSequenceSize = 0,
    size_t M = 1,
    size_t N = 0,
    size_t... Series
    >
struct fibonacci :
public fibonacci<SequenceSize, CurrentSequenceSize + 1, M + N, M, N, Series...>
{};

template <size_t SequenceSize, size_t M, size_t N, size_t First, size_t... Rest>
struct fibonacci<SequenceSize, SequenceSize, M, N, First, Rest...>
{
    constexpr static size_t value=First;
};


#endif

// trapezoidal_rule<> integrates a function F::f(x) between
// boundaries A::value and B::value with a partition size of
// N by means of the trapezoidal rule. See:
//
// http://en.wikipedia.org/wiki/Trapezoidal_rule
//
// Computations are intended to be done at compile time.
//
// Legend:
// F::f: integrand (function to be integrated)
// A::value: left boundary
// B::value: right boundary
// N: partition size
// I: initially equal to N-1. It is then
//    recursively decreased down to 0.
//


#ifndef TRAPEZOIDAL_RULE_H
#define TRAPEZOIDAL_RULE_H


#include <cstddef>


//
// Class: summation<>
//

template <typename F, typename A, typename B, size_t N, size_t I=N-1>
struct summation
{
    constexpr static double result=
        F::f(A::value+I*(B::value-A::value)/N)+summation<F, A, B, N, I-1>::result;
};

template <typename F, typename A, typename B, size_t N>
struct summation<F, A, B, N, 0>
{
    constexpr static double result=0.0;
};


//
// Class: trapezoidal_rule<>
//

template <typename F, typename A, typename B, size_t N>
struct trapezoidal_rule
{
    constexpr static double result=
        (B::value-A::value)/N*(0.5*(F::f(A::value)+F::f(B::value))+summation<F, A, B, N>::result);
};


#endif

上面的代码与阶乘的第一种实现方式类似。如果您想寻找与第二种实现方式的相似之处，这里是代码

// trapezoidal_rule_2<> yields the same result as trapezoidal_rule<>
// but it uses an indices unrolling technique instead of recursion.
//
// Legend:
// F::f: integrand (function to be integrated)
// A::value: left boundary
// B::value: right boundary
// N: partition size
// M: initially equal to N. Later it is unrolled
//      into the sequence: {1, 1, 2, 3, ..., N-1}.
//


#ifndef TRAPEZOIDAL_RULE_2_H
#define TRAPEZOIDAL_RULE_2_H


#include <cstddef>


//
// Class: summation_<>
//

// With the following definition, summation_<F, A, B, N> is recursively
// resolved as summation_<F, A, B, N, 1, 1, 2, 3, ..., N-1>.
template <typename F, typename A, typename B, size_t N, size_t M=N, size_t... Indices>
struct summation_ : public summation_<F, A, B, N, M-1, M-1, Indices...>    // [1]
{};

// The sequence above is split into the First index {1} plus the Rest...: {2, 3, ..., N-1}.
template <typename F, typename A, typename B, size_t N, size_t First, size_t... Rest>
struct summation_<F, A, B, N, 1, First, Rest...>                           // [2]
{
    constexpr static double result=
        F::f(A::value+First*(B::value-A::value)/N)+summation_<F, A, B, N, 1, Rest...>::result;
        // (Second addend evaluated at [2])
};

// And finally, this specialization matches the case in which all indices have been consumed.
template <typename F, typename A, typename B, size_t N>
struct summation_<F, A, B, N, 1>
{
    constexpr static double result=0.0;
};


//
// Class: trapezoidal_rule_2<>
//

template <typename F, typename A, typename B, size_t N>
struct trapezoidal_rule_2
{
    constexpr static double result=
      (B::value-A::value)/N*(0.5*(F::f(A::value)+F::f(B::value))+summation_<F, A, B, N>::result);
};


#endif

在我们理解了两种阶乘实现之后，我认为这两个最后一段代码是清晰、注释良好且自明了的。

辛普森法则（一种更精确的积分方法）的实现遵循相同的思路，但稍微冗长一些，因为它源于一个更复杂的公式。您可以在随附的 zip 文件中找到它。

源代码

梯形法则和辛普森法则的实现代码已附上，以及一个 main 示例代码。代码是用 Code::Blocks 编写的，并使用 MinGw gcc (v.4.8.1-4) 编译。

参考文献

1. 梯形法则。( http://en.wikipedia.org/wiki/Trapezoidal_rule)。
2. 辛普森法则。( http://en.wikipedia.org/wiki/Simpson's_rule)。
3. 编译时函数执行。( http://en.wikipedia.org/wiki/Compile_time_function_execution)。
4. Constexpr - C++11 中的通用常量表达式。( http://www.cprogramming.com/c++11/c++11-compile-time-processing-with-constexpr.html)。
5. “解包”元组以调用匹配的函数指针。( http://stackoverflow.com/questions/7858817/unpacking-a-tuple-to-call-a-matching-function-pointer)。