Bernstein polynomials

I’ve done a bit of searching and it looks like C implementations of the Bernstein basis polynomials that evaluate arbitrary derivatives and do not use recursion are very hard to obtain.

This seems quite useful to have, so I’m posting my own implementation with the hope that it might save someone else the work and possibly provide me with some improvement tips


#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);
}

I’ve tested it for sporadic x values with N = 2…5 and various derivatives and it all checks out. But if someone feels like verifying that would offcourse be appreciated.