#define REAL_TYPE double // alternatively float
#define REAL_ONE 1.0
#define REAL_ZERO 0.0
// returns 1 if i is even, -1 if i is odd
REAL_TYPE oddSign(const int i)
{
return i%2 == 0 ? REAL_ONE : -REAL_ONE;
}
// computes the _x'th derivative of the p'th monomial of x (assuming p positive)
REAL_TYPE diffMono(const int p, const int _x, const REAL_TYPE x)
{
const int k = p - _x;
if(k < 0)
return REAL_ZERO;
REAL_TYPE ret = REAL_ONE;
if(k > 0)
ret = pown(x,k);
for(int i = 0; i < _x; i++)
ret *= (REAL_TYPE)(p-i);
return ret;
}
// Computes the _x'th derivative of the n'th Bernstein polynomial basis funstion of degree N at x
REAL_TYPE bernstein(const int N, // polynomial degree
const int n, // polynomial ID
const int _x, // derivative
const REAL_TYPE x) // position
{
const int m = N-n;
const REAL_TYPE y = REAL_ONE - x;
if(_x > N)
return REAL_ZERO;
REAL_TYPE ret = REAL_ZERO;
for(int i = 0; i <= _x; i++)
ret += oddSign(i)*binCoef(_x,i)*diffMono(n,_x-i,x)*diffMono(m,i,y);
return ret*binCoef(N,n);
}