C++ source

C++ template code for Gaussian numerical integration

over n-dimensional cube [-1 .. 1] developed as C++ template metaprogram

Template parameters:

int ORD — order of the quadrature (must be >1)
int DIM — space dimension (must be >0)
typename FUN — integrand (may be a function or a functional object)

Control parameters:

CYCLE_INSTEAD_OF_STD_ALGO — if defined then cycles are used else std.algorithms
UNROLL_LOOPS_LIMIT — minimal array size to process by the cycle/std.algo (else by unrolled cycle)
NESTED_CYCLE_FOR_DIM_LESS_THAN_4 — if defined then ordinary nested cycles are used for dimensions 1-3

Other parameters:

real_t — real numbers type synonym

Source files:

Header (contains almost all the code)
cpp-file (contains the iterative calculation of the Gaussian quadrature coefficients that is necessary for orders >5)

How to use — a sample:

#include <iostream>
#include "GaussIntegr.h"

// integrand examples:
real_t fun3(real_t x, real_t y, real_t z) {return x*y*z;}
real_t fun3arr(const real_t* xyz) {return xyz[0]*xyz[1]*xyz[2];}
struct ffun3 {real_t operator()(real_t x, real_t y, real_t z)const {return x*y*z;} };

//integrands utilizing a vector data type:
struct vect
{real_t x,y,z;
vect(real_t v=0.): x(v), y(v),z(v) {}
real_t operator[](int idx)const {return idx==0 ? x : idx==1 ? y : z;}
real_t& operator[](int idx) {return idx==0 ? x : idx==1 ? y : z;}
vect& operator+=(const vect& v) {x+=v.x; y+=v.y; z+=v.z; return *this;}
vect& operator*=(real_t v) {x*=v; y*=v; z*=v; return *this;}
vect operator+(const vect& v)const { return (vect(*this))+=v;}
vect operator*(real_t v)const { return (vect(*this))*=v;}

real_t vfun3(vect v) {return v[0]*v[1]*v[2];}
vect vfun3v(const vect& v) {vect result(v); result.x+=v.x; result.y+=v.y; result.z+=v.z; return result;}


int main()
// creation of an integrating object for 3d-space and 4th order of Gauss quadrature:
static const fIntegrate<3,4> integr;

real_t result(0.);

std::cout<<integr.ByPlusAssgn_ArrArg(fun3arr,result)<<'\n'; //this method is absent if defined NESTED_CYCLE_FOR_DIM_LESS_THAN_4


std::cout<<integr.ByPlus_ArrArg<real_t>(fun3arr)<<'\n'; //is absent if defined NESTED_CYCLE_FOR_DIM_LESS_THAN_4


<< vect v;
<< integr.ByPlusAssgn<vect>(fvunv,v);

v = integr.ByPlus<vect,vect>(fvunv);

return 0;

Some test results (for the Gauss order=3):