Hi,

I'm educating myself on noise generation at the moment by implementing some algorithms. I am following this explanation and most of it was straight forward. I have Perlin Noise with the different interpolation functions they describe working nicely in 1-3 Dimensions. But the part about gradient Noise confuses me.
I get the 1D case. I choose a gradient at each point and interpolate the points in between with:

f(x) = (g2-3*g2)*x*x*x + (2*g1-g2)*x*x + g1
where g1 and g2 are the gradients (differentials in 1D) at the sample points.

But I don't see How I am supposed to do the interpolation in the multidimensional case. In the 2D case, do I have to find a function that is quadric in both x and y with those gradients? There are 16 parameters and I only have 12 conditions (f = 0 in the 4 corners and 8 partial derivatives)?

In the 2D case of gradient noise, you calculate the value of the gradient at a given vertex as the dot product of the fractional part of the input coords and a gradient looked up in a table:

[source]

double grad_noise_2(double x, double y, int ix, int iy, unsigned int seed)
{
unsigned int hash=hash_coords_2(ix,iy,seed);
double *vec=&gradient2D_lut[hash][0];

double dx=x-(double)ix;
double dy=y-(double)iy;

return (dx*vec[0] + dy*vec[1]);
}
[/source]

In the above (as taken from my own implementation) the fractional part, (dx,dy) is calculated by taking (x,y) - (floor(x), floor(y)). This is then dotted against the looked-up gradient. The result is a value in the range (-1,1). Four of these, one for each corner of the enclosing grid square, are then interpolated. The interpolation is done using the interpolation factors which are, coincidentally, also equal to (dx,dy) as above. You can either linearly interpolate using these factors as-is, or you can modify them according to a cubic or quintic function to smooth out the function:

[source]

double cubicSpline(double t)
{
return (t*t*(3-2*t));
}

double quinticSpline(double t)
{
return t*t*t*(t*(t*6-15)+10);
}


// Typedef for a noise generator function (value, gradient, etc...) in 2 dimensions

typedef double (*worker_noise_2)(double, double, int, int, unsigned int);


// General 2D interpolation functions

double interp_X_2(double x, double y, double xs, int x0, int x1, int iy, unsigned int seed, worker_noise_2 noisefunc)
{
double v1=noisefunc(x,y,x0,iy,seed);
double v2=noisefunc(x,y,x1,iy,seed);
return lerp(xs,v1,v2);
}

double interp_XY_2(double x, double y, double xs, double ys, int x0, int x1, int y0, int y1, unsigned int seed, worker_noise_2 noisefunc)
{
double v1=interp_X_2(x,y,xs,x0,x1,y0,seed,noisefunc);
double v2=interp_X_2(x,y,xs,x0,x1,y1,seed,noisefunc);
return lerp(ys,v1,v2);
}


// 2D Gradient noise function

double anl::gradient_noise2D(double x, double y, unsigned int seed, interp_func interp)
{
int x0=fast_floor(x);
int y0=fast_floor(y);
int x1=x0+1;
int y1=y0+1;

double xs=quinticSpline(x-(double)x0); // or cubicSpline(x-(double)x0);
double ys=quinticSpline(y-(double)y0); // or cubicSpline(y-(double)y0);

return interp_XY_2(x,y,xs,ys,x0,x1,y0,y1,seed,grad_noise_2);
}
[/source]

In 2D, it requires 3 interpolations. The top left and top right values are interpolated, then the bottom left and bottom right values are interpolated, both using xs as the interpolant. Then those two results are interpolated using ys.

If you are interested, the noise library found at the link in my signature is downloadable, and includes implementations of value, gradient and simplex noise in 2D, 3D, 4D and 6D variants. You can take a look at them to see how they are implemented.

I got it working after reading the sources of libnoise yesterday. Also after reading yours and their code I came to realize that I seem to have developed an excessively template heavy coding style... .

I appended the result of my foray

Basic Usage:

//create a noise object:
PerlinNoise2d<double> noise(persistence, octave_count);
//or
PerlinNoise2d<double> noise(persistence, octave_count, lacunarity);

//get some noise:
double result = noise(Vector<double, 2>(x, y));

//the multidimensional generators use gradient noise with quintic interpolation by default
//but this can be changed with template parameters:
PerlinNoise2d<double, ValueNoise2d<double> > noise(persistence, octave_count);
//or even
PerlinNoise2d<double, ValueNoise2d<double, IntNoise1, CosineInterpolator> > noise(persistence, octave_count);
//etc.


sources:
[source lang="cpp"]
/*
* Perlin.h - Copyright © 2011 Jakob Progsch
*
* This software is provided 'as-is', without any express or implied
* warranty. In no event will the authors be held liable for any damages
* arising from the use of this software.
*
* Permission is granted to anyone to use this software for any purpose,
* including commercial applications, and to alter it and redistribute it
* freely, subject to the following restrictions:
*
* 1. The origin of this software must not be misrepresented; you must not
* claim that you wrote the original software. If you use this software
* in a product, an acknowledgment in the product documentation would be
* appreciated but is not required.
*
* 2. Altered source versions must be plainly marked as such, and must not be
* misrepresented as being the original software.
*
* 3. This notice may not be removed or altered from any source
* distribution.
*/

/*
* provides value and gradient noise functions in 1-3 dimensions
* including perlin noise and ridge multifractals.
*/


#ifndef PERLIN_H
#define PERLIN_H

#include <cstdlib>
#include <cmath>

#include "include/MathVector.h"

const int FACTOR1 = 1619;
const int FACTOR2 = 31337;
const int FACTOR3 = 6971;

//pseudorandom generator
template<class T>
struct IntNoise1 {
inline T operator()(int x)
{
x = (x<<13) ^ x;
return ( 1.0 - ( (x * (x * x * 15731 + 789221) + 1376312589) & 0x7fffffff) / 1073741824.0);
}
inline int integer(int x)
{
x = (x<<13) ^ x;
return (x * (x * x * 15731 + 789221) + 1376312589) & 0x7fffffff;
}
};


//interpolation functors
template<class T>
struct CosineInterpolator {
inline T operator()(T a, T b, T x)
{
T f = 0.5*(1.-std::cos(x*3.14159265358979323846));
return a+(b-a)*f;
}
};

template<class T>
struct LinearInterpolator {
inline T operator()(T a, T b, T x)
{
return a+(b-a)*x;
}
};

template<class T>
struct CubicInterpolator {
inline T operator()(T a, T b, T x)
{
T f = x*x*(3-2*x);
return a+(b-a)*f;
}
};

template<class T>
struct QuinticInterpolator {
inline T operator()(T a, T b, T x)
{
T f = x*x*x*(10-(15-6*x)*x);
return a+(b-a)*f;
}
};

//noise functors
template<class T, template<class> class TNoise = IntNoise1, template<class> class TInterpolator = QuinticInterpolator>
struct ValueNoise1d {
T operator()(const T pos)
{
int integerpos;
T fractional;

integerpos = int(pos);
fractional = pos - integerpos;

return Interpolator(integerpos, integerpos+1, fractional);
}
private:
TNoise<T> Noise;
TInterpolator<T> Interpolator;
};

template<class T, template<class> class TNoise = IntNoise1, template<class> class TInterpolator = QuinticInterpolator>
struct ValueNoise2d {
template<class V>
T operator()(const V &pos)
{
Vector<int, 2> integerpos;
Vector<T, 2> fractional;
for(int d = 0;d<2;++d)
{
integerpos[d] = int(std::floor(pos[d]));
fractional[d] = pos[d] - integerpos[d];
}
T v00 = Noise(FACTOR1*(integerpos[0]) + FACTOR2*integerpos[1] );
T v10 = Noise(FACTOR1*(integerpos[0]+1) + FACTOR2*integerpos[1] );
T v01 = Noise(FACTOR1*(integerpos[0]) + FACTOR2*(integerpos[1]+1));
T v11 = Noise(FACTOR1*(integerpos[0]+1) + FACTOR2*(integerpos[1]+1));

T i0 = Interpolator(v00, v01, fractional[1]);
T i1 = Interpolator(v10, v11, fractional[1]);

return Interpolator(i0, i1, fractional[0]);
}
private:
TNoise<T> Noise;
TInterpolator<T> Interpolator;
};

template<class T, template<class> class TNoise = IntNoise1, template<class> class TInterpolator = QuinticInterpolator, unsigned table_size = 256>
struct GradientNoise2d {
GradientNoise2d() : gradients(table_size)
{
for(int i = 0;i<table_size;++i)
{
gradients = Vector<T, 2>(Noise(2*i), Noise(2*i+1));
gradients.normalize();
}
}

template<class V>
T operator()(const V &pos)
{
Vector<int, 2> integerpos;
Vector<T, 2> fractional;
for(int d = 0;d<2;++d)
{
integerpos[d] = int(std::floor(pos[d]));
fractional[d] = pos[d] - integerpos[d];
}
T v00 = GradNoise(pos, integerpos);
T v10 = GradNoise(pos, integerpos+Vector<int, 2>(1,0));
T v01 = GradNoise(pos, integerpos+Vector<int, 2>(0,1));
T v11 = GradNoise(pos, integerpos+Vector<int, 2>(1,1));

T i0 = Interpolator(v00, v01, fractional[1]);
T i1 = Interpolator(v10, v11, fractional[1]);

return Interpolator(i0, i1, fractional[0]);
}
private:
T GradNoise(Vector<T, 2> p, Vector<int, 2> i)
{
int n = Noise.integer(FACTOR1*i[0] + FACTOR2*i[1])%table_size;
Vector<T, 2> v;
for(int d = 0;d<2;++d)
v[d] = p[d] - i[d];
return 2.*dot(gradients[n], v);
}

TNoise<T> Noise;
TInterpolator<T> Interpolator;
std::vector< Vector<T, 2> > gradients;
};

template<class T, template<class> class TNoise = IntNoise1, template<class> class TInterpolator = QuinticInterpolator>
struct ValueNoise3d {
template<class V>
T operator()(const V &pos)
{
Vector<int, 3> integerpos;
Vector<T, 3> fractional;
for(int d = 0;d<3;++d)
{
integerpos[d] = int(std::floor(pos[d]));
fractional[d] = pos[d] - integerpos[d];
}
T v000 = Noise(FACTOR1*(integerpos[0]) + FACTOR2*integerpos[1] + FACTOR3*(integerpos[2]));
T v100 = Noise(FACTOR1*(integerpos[0]+1) + FACTOR2*integerpos[1] + FACTOR3*(integerpos[2]));
T v010 = Noise(FACTOR1*(integerpos[0]) + FACTOR2*(integerpos[1]+1) + FACTOR3*(integerpos[2]));
T v110 = Noise(FACTOR1*(integerpos[0]+1) + FACTOR2*(integerpos[1]+1) + FACTOR3*(integerpos[2]));
T v001 = Noise(FACTOR1*(integerpos[0]) + FACTOR2*integerpos[1] + FACTOR3*(integerpos[2]+1));
T v101 = Noise(FACTOR1*(integerpos[0]+1) + FACTOR2*integerpos[1] + FACTOR3*(integerpos[2]+1));
T v011 = Noise(FACTOR1*(integerpos[0]) + FACTOR2*(integerpos[1]+1) + FACTOR3*(integerpos[2]+1));
T v111 = Noise(FACTOR1*(integerpos[0]+1) + FACTOR2*(integerpos[1]+1) + FACTOR3*(integerpos[2]+1));

T i1 = Interpolator(v000, v001, fractional[2]);
T i2 = Interpolator(v010, v011, fractional[2]);
T i3 = Interpolator(v100, v101, fractional[2]);
T i4 = Interpolator(v110, v111, fractional[2]);

T j1 = Interpolator(i1, i2, fractional[1]);
T j2 = Interpolator(i3, i4, fractional[1]);

return Interpolator(j1, j2, fractional[0]);
}
private:
TNoise<T> Noise;
TInterpolator<T> Interpolator;
};


template<class T, template<class> class TNoise = IntNoise1, template<class> class TInterpolator = QuinticInterpolator, unsigned table_size = 256>
struct GradientNoise3d {
GradientNoise3d() : gradients(table_size)
{
for(int i = 0;i<table_size;++i)
{
gradients = Vector<T, 3>(Noise(3*i), Noise(3*i+1), Noise(3*i+2));
gradients.normalize();
}
}

template<class V>
T operator()(const V &pos)
{
Vector<int, 3> integerpos;
Vector<T, 3> fractional;
for(int d = 0;d<3;++d)
{
integerpos[d] = int(std::floor(pos[d]));
fractional[d] = pos[d] - integerpos[d];
}
T v000 = GradNoise(pos, integerpos);
T v100 = GradNoise(pos, integerpos+Vector<int, 3>(1,0,0));
T v010 = GradNoise(pos, integerpos+Vector<int, 3>(0,1,0));
T v110 = GradNoise(pos, integerpos+Vector<int, 3>(1,1,0));
T v001 = GradNoise(pos, integerpos+Vector<int, 3>(1,0,1));
T v101 = GradNoise(pos, integerpos+Vector<int, 3>(1,0,1));
T v011 = GradNoise(pos, integerpos+Vector<int, 3>(0,1,1));
T v111 = GradNoise(pos, integerpos+Vector<int, 3>(1,1,1));

T i1 = Interpolator(v000, v001, fractional[2]);
T i2 = Interpolator(v010, v011, fractional[2]);
T i3 = Interpolator(v100, v101, fractional[2]);
T i4 = Interpolator(v110, v111, fractional[2]);

T j1 = Interpolator(i1, i2, fractional[1]);
T j2 = Interpolator(i3, i4, fractional[1]);

return Interpolator(j1, j2, fractional[0]);
}
private:
T GradNoise(Vector<T, 3> p, Vector<int, 3> i)
{
int n = Noise.integer(FACTOR1*i[0] + FACTOR2*i[1] + FACTOR3*i[2])%table_size;
Vector<T, 2> v;
for(int d = 0;d<3;++d)
v[d] = p[d] - i[d];
return 2.*dot(gradients[n], v);
}

TNoise<T> Noise;
TInterpolator<T> Interpolator;
std::vector< Vector<T, 3> > gradients;
};

//perlin noise and ridge multifractal functors

template<class T, class TNoise = GradientNoise3d<T> >
struct RidgeMultifractal3d {
RidgeMultifractal3d(T p, int o) : persistence(p), lacunarity(2), octaves(o) { }
RidgeMultifractal3d(T p, int o, T l) : persistence(p), lacunarity(l), octaves(o) { }

template<class V>
T operator()(const V &pos)
{
T result = 0;
T frequency = 1, amplitude = 1;
T scaler = 0;
for(int i=0;i<octaves;++i)
{
result += (1-2*std::abs(noise(pos*frequency)))*amplitude;
scaler += amplitude;
frequency *= lacunarity;
amplitude *= persistence;
}
return result/scaler;
}
private:
T persistence;
T lacunarity;
int octaves;
TNoise noise;
};

template<class T, class TNoise = GradientNoise3d<T> >
struct PerlinNoise3d {
PerlinNoise3d(T p, int o) : persistence(p), lacunarity(2), octaves(o) { }
PerlinNoise3d(T p, int o, T l) : persistence(p), lacunarity(l), octaves(o) { }

template<class V>
T operator()(const V &pos)
{
T result = 0;
T frequency = 1, amplitude = 1;
T scaler = 0;
for(int i=0;i<octaves;++i)
{
result += noise(pos*frequency)*amplitude;
scaler += amplitude;
frequency *= lacunarity;
amplitude *= persistence;
}
return result/scaler;
}
private:
T persistence;
T lacunarity;
int octaves;
TNoise noise;
};


template<class T, class TNoise = GradientNoise2d<T> >
struct RidgeMultifractal2d {
RidgeMultifractal2d(T p, int o) : persistence(p), lacunarity(2), octaves(o) { }
RidgeMultifractal2d(T p, int o, T l) : persistence(p), lacunarity(l), octaves(o) { }

template<class V>
T operator()(const V &pos)
{
T result = 0;
T frequency = 1, amplitude = 1;
T scaler = 0;
for(int i=0;i<octaves;++i)
{
result += (1-2*std::abs(noise(pos*frequency)))*amplitude;
scaler += amplitude;
frequency *= lacunarity;
amplitude *= persistence;
}
return result/scaler;
}
private:
T persistence;
T lacunarity;
int octaves;
TNoise noise;
};

template<class T, class TNoise = GradientNoise2d<T> >
struct PerlinNoise2d {
PerlinNoise2d(T p, int o) : persistence(p), lacunarity(2), octaves(o) { }
PerlinNoise2d(T p, int o, T l) : persistence(p), lacunarity(l), octaves(o) { }

template<class V>
T operator()(const V &pos)
{
T result = 0;
T frequency = 1, amplitude = 1;
T scaler = 0;
for(int i=0;i<octaves;++i)
{
result += noise(pos*frequency)*amplitude;
scaler += amplitude;
frequency *= lacunarity;
amplitude *= persistence;
}
return result/scaler;
}
private:
T persistence;
T lacunarity;
int octaves;
TNoise noise;
};


template<class T, class TNoise = ValueNoise1d<T> >
struct PerlinNoise1d {
PerlinNoise1d(T p, int o) : persistence(p), lacunarity(2), octaves(o) { }
PerlinNoise1d(T p, int o, T l) : persistence(p), lacunarity(l), octaves(o) { }

T operator()(const T pos)
{
T result = 0;
T frequency = 1, amplitude = 1;
T scaler = 0;
for(int i=0;i<octaves;++i)
{
result += noise(pos*frequency)*amplitude;
scaler += amplitude;
frequency *= lacunarity;
amplitude *= persistence;
}
return result/scaler;
}
private:
T persistence;
T lacunarity;
int octaves;
TNoise noise;
};

#endif
[/source]

And the MathVector.h file in case someone actually wants to try it:
[source lang="cpp"]
/*
* Vector.h - Copyright © 2011 Jakob Progsch
*
* This software is provided 'as-is', without any express or implied
* warranty. In no event will the authors be held liable for any damages
* arising from the use of this software.
*
* Permission is granted to anyone to use this software for any purpose,
* including commercial applications, and to alter it and redistribute it
* freely, subject to the following restrictions:
*
* 1. The origin of this software must not be misrepresented; you must not
* claim that you wrote the original software. If you use this software
* in a product, an acknowledgment in the product documentation would be
* appreciated but is not required.
*
* 2. Altered source versions must be plainly marked as such, and must not be
* misrepresented as being the original software.
*
* 3. This notice may not be removed or altered from any source
* distribution.
*/

/*
* Vector.h provides a simple static vector template class with
* a basic expression template ansatz for vector operations.
*/

#ifndef VECTOR_H
#define VECTOR_H

#include <cmath>
#include <algorithm>
#include <ostream>

template<class T, unsigned D> class Vector;

//base class for all expression templates
template<class T, unsigned D, class A>
class VectorExpr {
public:
inline operator const A&() const
{
return *static_cast<const A*>(this);
}
};

//better use macros instead of copy pasting this stuff all over the place
#define MAKE_VEC_VEC_EXPRESSION(NAME, EXPR, FUNCTION) \
template<class T, unsigned D, class A, class B> \
class NAME : public VectorExpr<T, D, NAME<T, D, A, B> > { \
public: \
NAME(const A& pa, const B& pb) : a(pa), b(pb) { } \
inline T operator[](unsigned i) const { return EXPR; } \
private: \
const A& a; \
const B& b; \
}; \
\
template<class T, unsigned D, class A, class B> \
inline NAME<T,D,A,B> \
FUNCTION(const VectorExpr<T,D,A> &a, const VectorExpr<T,D,B> &b) \
{ \
return NAME<T,D,A,B>(a, b); \
}

#define MAKE_VEC_SCAL_EXPRESSION(NAME, EXPR, FUNCTION) \
template<class T, unsigned D, class A> \
class NAME : public VectorExpr<T, D, NAME<T, D, A> > { \
public: \
NAME(const A& pa, const T& pb) : a(pa), b(pb) { } \
inline T operator[](unsigned i) const { return EXPR; } \
private: \
const A& a; \
const T& b; \
}; \
\
template<class T, unsigned D, class A> \
inline NAME<T,D,A> \
FUNCTION(const VectorExpr<T,D,A> &a, const T &b) \
{ \
return NAME<T,D,A>(a, b); \
}

#define MAKE_SCAL_VEC_EXPRESSION(NAME, EXPR, FUNCTION) \
template<class T, unsigned D, class A> \
class NAME : public VectorExpr<T, D, NAME<T, D, A> > { \
public: \
NAME(const T& pa, const A& pb) : a(pa), b(pb) { } \
inline T operator[](unsigned i) const { return EXPR; } \
private: \
const T& a; \
const A& b; \
}; \
\
template<class T, unsigned D, class A> \
inline NAME<T,D,A> \
FUNCTION(const T &a, const VectorExpr<T,D,A> &b) \
{ \
return NAME<T,D,A>(a, b); \
}

#define MAKE_VEC_EXPRESSION(NAME, EXPR, FUNCTION) \
template<class T, unsigned D, class A> \
class NAME : public VectorExpr<T, D, NAME<T, D, A> > { \
public: \
NAME(const A& pa) : a(pa) { } \
inline T operator[](unsigned i) const { return EXPR; } \
private: \
const A& a; \
}; \
\
template<class T, unsigned D, class A> \
inline NAME<T,D,A> \
FUNCTION(const VectorExpr<T,D,A> &a) \
{ \
return NAME<T,D,A>(a); \
}

//create actual functions and operators
MAKE_VEC_VEC_EXPRESSION (EMulExpr, a * b, multiply_elements)
MAKE_VEC_VEC_EXPRESSION (EDivExpr, a / b, divide_elements)
MAKE_VEC_VEC_EXPRESSION (AddExpr, a + b, operator+)
MAKE_VEC_VEC_EXPRESSION (SubExpr, a - b, operator-)
MAKE_VEC_SCAL_EXPRESSION(DivExpr, a / b, operator/)
MAKE_VEC_SCAL_EXPRESSION(MulExpr1, a * b, operator*)
MAKE_SCAL_VEC_EXPRESSION(MulExpr2, a * b, operator*)
MAKE_VEC_EXPRESSION (NegExpr, -a, operator-)

//static size assertion since the vector size ist also static
template<unsigned I, unsigned J>
struct STATIC_DIMENSION_MISMATCH_ASSERTION;

template<unsigned I>
struct STATIC_DIMENSION_MISMATCH_ASSERTION<I,I> { };

#define ASSERT_DIMENSION(I, J) sizeof(STATIC_DIMENSION_MISMATCH_ASSERTION<I,J>)

//actual vector class
template<class T, unsigned D>
class Vector : public VectorExpr<T, D, Vector<T,D> > {
public:
static const unsigned Dim = D;
typedef T Type;

Vector()
{
std::fill(data, data+Dim, T());
}

Vector(const T &p1)
{
ASSERT_DIMENSION(1, Dim);
data[0] = p1;
}

Vector(const T &p1, const T &p2)
{
ASSERT_DIMENSION(2, Dim);
data[0] = p1;
data[1] = p2;
}

Vector(const T &p1, const T &p2, const T &p3)
{
ASSERT_DIMENSION(3, Dim);
data[0] = p1;
data[1] = p2;
data[2] = p3;
}

Vector(const T &p1, const T &p2, const T &p3, const T &p4)
{
ASSERT_DIMENSION(4, Dim);
data[0] = p1;
data[1] = p2;
data[2] = p3;
data[3] = p4;
}

template<class A>
Vector(const VectorExpr<T, D, A>& a)
{
const A& ao ( a );
for(unsigned i = 0;i<Dim;++i)
data = ao;
}

T& operator[] (unsigned i) { return data; }
const T& operator[] (unsigned i) const { return data; }

const Vector& operator*=(const T &b)
{
for(unsigned i = 0;i<Dim;++i)
data *= b;
return *this;
}

const Vector& operator/=(const T &b)
{
for(unsigned i = 0;i<Dim;++i)
data /= b;
return *this;
}

template<class A>
const Vector& operator+=(const VectorExpr<T, D, A>& a)
{
const A& ao ( a );
for(unsigned i = 0;i<Dim;++i)
data += ao;
return *this;
}

template<class A>
const Vector& operator-=(const VectorExpr<T, D, A>& a)
{
const A& ao ( a );
for(unsigned i = 0;i<Dim;++i)
data -= ao;
return *this;
}

template<class A>
Vector& operator=(const VectorExpr<T, D, A>& a)
{
const A& ao ( a );
for(unsigned i = 0;i<Dim;++i)
data = ao;
return *this;
}

Vector& normalize()
{ *this /= abs(*this); return *this; }
private:
T data[Dim];
};

template<class T, unsigned D, class A>
inline Vector<T, D> eval(const VectorExpr<T, D, A>& a)
{
return Vector<T, D>(a);
}

//"reduction" functions that don't return expression templates
template<class T, unsigned D, class A>
inline T sum(const VectorExpr<T, D, A>& a)
{
const A& ao ( a );
T res = 0;
for(unsigned i = 0;i<D;++i)
res += ao;
return res;
}

template<class T, unsigned D, class A, class B>
inline T dot(const VectorExpr<T, D, A>& a, const VectorExpr<T, D, B>& b)
{
return sum(multiply_elements(a, b));
}

template<class T, unsigned D, class A>
inline T squared_norm(const VectorExpr<T, D, A>& a)
{
const A& ao ( a );
T res = 0;
for(unsigned i = 0;i<D;++i)
{
T tmp = ao;
res += tmp*tmp;
}
return res;
}

template<class T, unsigned D, class A>
inline T abs(const VectorExpr<T, D, A>& a)
{
return std::sqrt(squared_norm(a));
}

template<class T, unsigned D, class A>
std::ostream& operator<<(std::ostream& out, const VectorExpr<T, D, A>& a)
{
const A& ao ( a );
T res = 0;
out << '(' << ao[0];
for(unsigned i = 1;i<D;++i)
{
out << ", " << ao;
}
out << ')';

return out;
}

#undef MAKE_VEC_VEC_EXPRESSION
#undef MAKE_VEC_SCAL_EXPRESSION
#undef MAKE_SCAL_VEC_EXPRESSION
#undef MAKE_VEC_EXPRESSION

#endif

[/source]