/Users/mourrain/Devel/mmx/realroot/include/realroot/binomials.hpp

Go to the documentation of this file.

#ifndef realroot_SOLVE_SBDSLV_BINOMIALS_H
#define realroot_SOLVE_SBDSLV_BINOMIALS_H

#include <vector>
#include <assert.h>
//--------------------------------------------------------------------
namespace mmx {

#ifndef WITH_POW
#define WITH_POW
  template<class T>
  T pow(const T& a, int i)
  {
    assert(i>=0);
    if(i == 1) return a;
    if(i > 0)
      {
        if (a == T(1))  return a;
        if (a == T(-1)) return (i % 2  == 0) ? T(1) : T(-1);
        T y(a);
        T z(1);
        unsigned int n = i;
        unsigned int m;
        while (n > 0) {
          m = n / 2;
          if (n %2 ) 
            {
              z *= y;
            }
          y = y * y;
          n = m;
        }
        return z;
      }
    return T(1);
  }
#endif //WITH_POW
//----------------------------------------------------------------------
template <typename T>  T factorial(int n)
{
  T r=1;
  for (int i=2;i<=n;i++) r*=n;
  return(r);
}


#ifndef WITH_BINOMIAL
#define WITH_BINOMIAL

template <typename T>  T binomial(int n, int p)
{
  T n1=1;

  if(p!=0)
    {
      if((n-p)<p) p=n-p;
      int n2=n;
      for(int i=1;i<=p;i++,n2--) n1=n1*n2/T(i);
    }

  return(n1);
}
#endif //WITH_BINOMIAL
//----------------------------------------------------------------------
template<typename C>
struct binomials {
  static binomials default_;
  std::vector<int> offs;
  std::vector<C>   data;
  binomials(const binomials& b):offs(b.offs),data(b.data) {}
  binomials();
  // return the list [ binomial(d,i) ], for i = 0..d
  const C* operator[]( int d );
  const C& operator()( int i, int j ) { return (*this)[i][j]; };
};

// IMPLEMENTATION

template<class C>
binomials<C>::binomials() : offs(2), data(3) 
{
  offs[0] = 0;       data[0] = 1;
  offs[1] = 1; data[offs[1]] = 1; data[offs[1]+1] = 1;
};

template<class C> binomials<C> binomials<C>::default_;

template<class C> 
const C* binomials<C>::operator[]( int d )
{
  int s,j,x;
  if ( (s=offs.size()) > d ) return &*(data.begin())+offs[d];
  for ( x = 0, j = s; j <= d; x += j+1, j ++ ) {}
  offs.resize(d+1);
  data.resize(data.size()+x);
  for ( j = s; j <= d; j ++ )
    {
      offs[j] = offs[j-1] + j;
      for ( int i = 0; i < j-1; i ++ )
        data[offs[j]+i+1] = data[offs[j-1]+i] + data[offs[j-1]+i+1];
      data[offs[j]] = data[offs[j]+j] = 1;
    };
  return &*(data.begin())+offs[d];
};

} //namespace mmx

#endif