// ************************************************************************* //
/// AUTHOR Alexis Angelidis
/// STARTED Thu May 26 15:35:24 NZST 2005
// ************************************************************************* //

// ************************************************************************* //
#include <iostream>

// ************************************************************************* //
/// DESCRIPTION: Computes q to the power k.
// ************************************************************************* //
template <class T>
T	pow(const T	q,
	    const int	k)
{
  T	ret(1), r(q);

  for (int i(1); i <= k; i *= 2, r *= r)
    if (k & i)
      ret *= r;
  return ret;
}

// ************************************************************************* //
/// DESCRIPTION: Computes factorial i.
// ************************************************************************* //
inline int factorial(const int i)
{
  int	ret(1);

  for (int k(1); k <= i; ++k)
    ret *= k;
  return ret;
}

// ************************************************************************* //
/// DESCRIPTION:
///  Bernstein polynomial of degree N in L variables.
///  L	- number of variables (bivariate, trivariate, ...)
///      1 == point // 2 == line // 3 = triangle // 4 = quad // ...
///  N	- +1 = number of control points (degree of polynomial)
/// CONTEXT: Sum(i_k) = N
///  i	- index of curve within lattice
///  u	- barycentric coordinates, Sum(u_k) = 1.0
// ************************************************************************* //
template <int L, int N>
double	bernstein(const int	i[L],
		  const double	u[L])
{
  double	ret(1.0);
  int		iret(1);

  for (int k(0); k < L; ++k)
    {
      ret *= pow(u[k], i[k]);
      iret *= factorial(i[k]);
    }
  return factorial(N) * ret / iret;
}

// ************************************************************************* //
template <int L, int N, class T>
void	recBezierSimplex(T		&tot,
			 int		&count,			 
			 const double	u[L],
			 const T	*const b,
			 int		wordBuffer[L],
			 const int	sumword = 0,
			 const int	i = 0)
{
  if (i == L)
    {
      if (sumword == N)
	{
	  tot += b[count] * bernstein<L, N>(wordBuffer, u);
	  ++count;
	}
      return;
    }
  for (int k(0); k <= N; ++k)
    {
      if (N < sumword + k)
	break;
      wordBuffer[i] = k;
      recBezierSimplex<L, N>(tot, count, u, b, wordBuffer, sumword + k, i + 1);
    }
}

// ************************************************************************* //
/// DESCRIPTION: Multivariate bezier simplex of degree N in L variables.
/// b	- interpolated value, at sites, listed in alphabetical order
///       Order example <l = 2, n = 2>: [02][11][20]
/// u	- barycentric coordinates within lattice (or spherical barycentric)
///  L	- number of variables (bivariate, trivariate, ...)
///      1 == point // 2 == line // 3 = surface // 4 = volume // ...
///  N	- + 1 = number of control points per side (degree of polynomial)
/// EXAMPLES OF CONTROL NETS:
/// L = 2 n = 2        l = 2, n = 3
/// b02 b11 b20        b03 b12 b21 b30
/// 
/// L = 3 n = 2        l = 3, n = 3
///       b002                 b003
///    b011  b101           b012  b021
/// b020  b110  b200     b021  b111 b012
///                    b030 b120 b210 b300
/// L = 4 n = 2 (Tetrahedron)
///  b0002 b0011 b0020
///        b1001      
///  b0101       b1010
///        b0110
///  b0200 b1100 b2000
// ************************************************************************* //
template <int L, int N, class T>
T	bezierSimplex(const double	u[L],
		      const T		*const b)
{
  int	wordBuffer[L];
  int	count(0);
  T	tot(0.0);

  recBezierSimplex<L, N>(tot, count, u, b, wordBuffer);
  return tot;
}

// ************************************************************************* //
int main()
{
  double u[3] = { 0.0, 1.0, 0.0 };
  double b[6] = { 2.2, 2.0, 2.0, 2.0, 2.0, 2.0 };

  std::cout << bezierSimplex<3, 3>(u, b) << std::endl;
  return -1;
}