Click here to Skip to main content
Click here to Skip to main content
Go to top

Numerical integration at compile time

, 18 Feb 2013
Rate this:
Please Sign up or sign in to vote.
You have not lifted your finger off the mouse button yet when the calculation is served.

Introduction

It is possible to carry out far from trivial calculations at compile time. As a sample, in this article I will show how to calculate the area under any real function between arbitrary boundaries A and B. I will use standard and widespread algorithms, like trapezoidal rule and Simpson’s rule.

Computing at compile time means that the result of the operation is known by the time the process is loaded in memory, before entering the main function.

The magic of the main function

Have a look at this code:

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

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

double F(double x) {return cos(x);}  // Integrand
extern const double A = 0.0;         // Left boundary
extern const double B = PI / 2.0;    // Right boundary

 
int main(int argn, char *argc[])
{
    // Application of the Simpson's rule.
    // F is integrated from A to B with a partition size of 100.
    std::cout << simpsons_rule<F, A, B, 100>::result << std::endl;
 
    return 0;
} 

The program above simply displays the value stored in a static variable result belonging to the class simpsons_rule<>. This value is evaluated at compile time and stored in memory before main() is accessed. result stores the integral of cos(x) (function F) between 0.0 (A) and PI/2.0 (B). The exact result is 1. The displayed result is also 1. Thus, the approximation is almost perfect in this case (think about round-off errors).

If you now change F, A or B and recompile, you will get a different value for result.

So this is the crux and the magic of the matter: at runtime no computations are done; the only performed task is displaying the result.

Note that in order to use F, A and B as template arguments, they must have external linkage; they are globals and their values are known at compile time.

Now it is time to see the guts, but instead of diving into trapezoidal or Simpson’s rules implementations, it is much more illustrative to show a simpler (but analogous) case: the factorial.

Implementing the factorial (1)

To compute the factorial of an integer n (noted as n!), just apply the following formula:

The following code implements the factorial computation at compile time:

#ifndef FACTORIAL_H
#define FACTORIAL_H
 
 
#include <cstddef>

 
// [1]
template <size_t N>
struct factorial
{
    static const size_t result;
};
 
 
// [2]
template <size_t N>
const size_t factorial<N>::result = N * factorial<N - 1>::result;
 
 
// [3]
template <>
struct factorial<1>
{
    static const size_t result;
};
 
 
// [4]
const size_t factorial<1>::result = 1;
 
 
#endif

If you browse in other sources (ref. [3], for example), you will probably find less verbose implementations. Please, stay with this for the moment: it will help you find resemblances to trapezoidal and Simpson’s rules implementations later.

The code above declares a class template named factorial<> [1] and a specialization factorial<1> [3]. In [2] the static unsigned integer result is defined for factorial<>. In [4] the static unsigned integer result is defined for factorial<1>.

Now suppose you write this sentence in the main function:

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

As you can see in [2], the compiler resolves factorial<5>::result as:

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

But also in [2] factorial<4>::result is resolved as:

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

Following this reasoning in a recursive way we get:

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

factorial<1>::result is resolved in [4] as 1. We get then:

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

The compiler finally carries out the multiplication above, and factorial<5>::result is resolved at compile time as 120.

Implementing the factorial (2)

The following code is a not so known alternative, only valid since 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]
{};
 
 
// Thus, in the following specialization, Indices... corresponds
// to the sequence: {1, 2, 3, ..., N}.
template <size_t N, size_t... Indices>
struct factorial<N, 1, Indices...>                                 // [2]
{
    static const size_t result;
};
 
 
// 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...>                             // [3]
{
    static const size_t result;
};
 
 
// The indices in the sequence are recursively evaluated.
template <size_t N, size_t First, size_t... Rest>                  // [4]
const size_t factorial<N, 1, First, Rest...>::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>                                                // [5]        
struct factorial<N, 1>
{
    static const size_t result;
};
 
 
template <size_t N>                                                // [6]    
const size_t factorial<N, 1>::result = 1;
 
 
#endif 

How does the compiler resolve factorial<5>::result in this case? If you have a look at [1], you will see that factorial<5> can be resolved into factorial<5, 5, 5>, being Indices… an empty pack. Without taking your eyes off [1], you can follow the chain:

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}. 

Once the compiler gets factorial<5, 1, 1, 2, 3, 4, 5> it can be resolved by the specialization at [2]. It becomes then factorial<5, 1, Indices…>, being Indices… the pack: {1, 2, 3, 4, 5}.

The technique above is very interesting because it allows us to get a pack of indices for the operands involved in our computation. In the case of the factorial, indices and operands (numbers to be multiplied) are the same, but this is just a particular case.

In [3] our indices pack is split in two: the first index and the rest. You can find in [4] the definition of result for the specialization in [3]. It is self-explanatory.

The termination condition is reached when Rest… is empty. Such case is implemented in [5] and [6].

After following the previous steps, the compiler finally resolves factorial<5>::result as 120. As you can imagine, both methodologies come together to the same binary output.

An esoteric Fibonacci sequence

As another illustrative example, you can use the technique above to compute an element in a Fibonacci sequence. Be indulgent with me if it looks weird to your eyes:

#ifndef FIBONACCI_H
#define FIBONACCI_H


#include <cstddef>


template <
          size_t SequenceSize,
          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...>
{
    static const size_t value;
};


template <size_t SequenceSize, size_t M, size_t N, size_t First, size_t... Rest>
const size_t fibonacci<SequenceSize, SequenceSize, M, N, First, Rest...>::value = First;


#endif 

If you type (for instance) fibonacci<5>::value in the main function, it will be resolved as 3 (the fifth element in the series). I agree this is not the most intuitive way to compute a Fibonacci sequence, of course. 

The trapezoidal rule

Observe that factorial of n can be written schematically as:

where f(x)=x and {i1, i2,…, in} is a collection of indices with values ranging from 1 to n.

The trapezoidal rule can be arranged into a similar pattern:

where: 

Definition of g()

In the code below, summation<> class template carries out the summation of the g()’s:

// trapezoidal_rule<> integrates a function F(x) between
// boundaries A and B 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: integrand (function to be integrated)
// A: left boundary
// B: 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>
 
 
typedef double (*integrand_type)(double);
typedef const double& boundary_type;
 
 
//
// Class: summation<>
//

template <integrand_type F, boundary_type A, boundary_type B, size_t N, size_t I = N - 1>
struct summation
{
    static const double result;
};
 
 
template <integrand_type F, boundary_type A, boundary_type B, size_t N, size_t I>
const double summation<F, A, B, N, I>::result =
       F(A + I * (B - A) / N) + summation<F, A, B, N, I - 1>::result;
 
 
template <integrand_type F, boundary_type A, boundary_type B, size_t N>
struct summation<F, A, B, N, 0>
{
    static const double result;
};
 
 
template <integrand_type F, boundary_type A, boundary_type B, size_t N>
const double summation<F, A, B, N, 0>::result = 0.0;
 
 
//
// Class: trapezoidal_rule<>
//

template <integrand_type F, boundary_type A, boundary_type B, size_t N>
struct trapezoidal_rule
{
    static const double result;
};
 
 
template <integrand_type F, boundary_type A, boundary_type B, size_t N>
const double trapezoidal_rule<F, A, B, N>::result =
       (B - A) / N * (0.5 * (F(A) + F(B)) + summation<F, A, B, N>::result);
 
 
#endif 

The code above bears resemblance to the first implementation of factorial. If you look for resemblance to the second implementation, here you are the code:

// trapezoidal_rule_2<> yields the same result as trapezoidal_rule<>
// but it uses an indices unrolling technique instead of recursion.
//
// Legend:
// F: integrand (function to be integrated)
// A: left boundary
// B: 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>

 
typedef double (*integrand_type)(double);
typedef const double& boundary_type;
 
 
//
// 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 <integrand_type F, boundary_type A, boundary_type 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]
{};
 
 
// Thus, in the following specialization, Indices... corresponds to the
// sequence: {1, 2, 3, ..., N - 1}.
template <integrand_type F, boundary_type A, boundary_type B, size_t N, size_t... Indices>
struct summation_<F, A, B, N, 1, Indices...>                                 // [2]
{
    static const double result;
};
 
 
// The sequence above is split into the First index {1} plus the Rest...: {2, 3, ..., N - 1}.
template <integrand_type F, boundary_type A, boundary_type B, size_t N, size_t First, size_t... Rest>
struct summation_<F, A, B, N, 1, First, Rest...>                             // [3]
{
    static const double result;
};
 
 
// The indices in the sequence are recursively evaluated.
template <integrand_type F, boundary_type A, boundary_type B, size_t N, size_t First, size_t... Rest>
const double summation_<F, A, B, N, 1, First, Rest...>::result =
       F(A + First * (B - A) / 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 <integrand_type F, boundary_type A, boundary_type B, size_t N>
struct summation_<F, A, B, N, 1>
{
    static const double result;
};
 
 
template <integrand_type F, boundary_type A, boundary_type B, size_t N>
const double summation_<F, A, B, N, 1>::result = 0.0;
 
 
//
// Class: trapezoidal_rule_2<>
//

template <integrand_type F, boundary_type A, boundary_type B, size_t N>
struct trapezoidal_rule_2
{
    static const double result;
};
 
 
template <integrand_type F, boundary_type A, boundary_type B, size_t N>
const double trapezoidal_rule_2<F, A, B, N>::result =
       (B - A) / N * (0.5 * (F(A) + F(B)) + summation_<F, A, B, N>::result);
 
 
#endif 

Once we have understood both factorial implementations, I think these two last pieces of code are clear, well commented and self-explanatory.

The implementation for Simpson’s rule (a more precise integration method) follows the same ideas but it is a little more verbose, because it arises from a bit more complicated formula. You can find it in the enclosed zip.

Source code

The trapezoidal and Simpson’s rules implementation code has been attached, along with a main sample code. The code has been written with Code::Blocks and compiled with MinGw (mingw-get-inst-20120426).

References

  1. Trapezoidal rule. (http://en.wikipedia.org/wiki/Trapezoidal_rule).
  2. Simpson’s rule. (http://en.wikipedia.org/wiki/Simpson's_rule).
  3. Compile time function execution. (http://en.wikipedia.org/wiki/Compile_time_function_execution).
  4. “Unpacking” a tuple to call a matching function pointer. (http://stackoverflow.com/questions/7858817/unpacking-a-tuple-to-call-a-matching-function-pointer).

License

This article, along with any associated source code and files, is licensed under The Code Project Open License (CPOL)

Share

About the Author

David Serrano Martínez
Systems Engineer
Spain Spain
I work as a senior industrial engineer for Public Administration in Spain. I have experience in developing software for ballistic computations. I like maths and programming and, above all, riding my mountain bike. I run my own website on computation and numerical calculus at rumigaculum.com. Contact me here.

Comments and Discussions

 
QuestionCompiler limit on recursions Pinmvpgeoyar28-Feb-13 16:48 
AnswerRe: Compiler limit on recursions PinmemberDavid Serrano Martínez28-Feb-13 21:02 
GeneralRe: Compiler limit on recursions Pinmvpgeoyar1-Mar-13 5:34 
GeneralMy vote of 5 PinmvpMichael Haephrati מיכאל האפרתי20-Feb-13 8:38 
GeneralRe: My vote of 5 PinmemberDavid Serrano Martínez22-Feb-13 8:59 
GeneralIntegration not at compile time Pinmemberarhimede19-Feb-13 5:45 
GeneralRe: Integration not at compile time PinmemberDavid Serrano Martínez19-Feb-13 7:50 
GeneralRe: Integration not at compile time Pinmemberarhimede19-Feb-13 9:04 
GeneralRe: Integration not at compile time PinmemberDavid Serrano Martínez22-Feb-13 8:58 
GeneralRe: Integration not at compile time Pinmemberarhimede22-Feb-13 20:26 
GeneralRe: Integration not at compile time PinmemberKevin Drzycimski25-Feb-13 1:35 
GeneralRe: Integration not at compile time PinmemberDavid Serrano Martínez25-Feb-13 4:01 
GeneralRe: Integration not at compile time Pinmember.:floyd:.10-Mar-14 3:47 
GeneralRe: Integration not at compile time PinmemberDavid Serrano Martínez10-Mar-14 4:32 
GeneralRe: Integration not at compile time Pinmemberjohn morrison leon28-Aug-14 2:36 
GeneralMy vote of 5 PinmemberKevin Drzycimski18-Feb-13 23:50 
GeneralRe: My vote of 5 PinmemberDavid Serrano Martínez19-Feb-13 3:32 
GeneralMy vote of 5 PinmemberJohnWallis4218-Feb-13 17:08 
GeneralRe: My vote of 5 PinmemberDavid Serrano Martínez19-Feb-13 3:35 

General General    News News    Suggestion Suggestion    Question Question    Bug Bug    Answer Answer    Joke Joke    Rant Rant    Admin Admin   

Use Ctrl+Left/Right to switch messages, Ctrl+Up/Down to switch threads, Ctrl+Shift+Left/Right to switch pages.

| Advertise | Privacy | Mobile
Web02 | 2.8.140916.1 | Last Updated 18 Feb 2013
Article Copyright 2013 by David Serrano Martínez
Everything else Copyright © CodeProject, 1999-2014
Terms of Service
Layout: fixed | fluid