ABACUS/src/COMBI/Combinatorics.cc

/****************************************************************

This software is part of J.-S. Caux's C++ library.

Copyright (c) 2006.

-----------------------------------------------------------

Combinatorics.cc

Defines all class related to combinatorics.

LAST MODIFIED:  04/09/06

******************************************************************/

#include "JSC.h"

using namespace std;

namespace JSC {

  Choose_Table::Choose_Table ()
    : Npower(0ULL), Npowerp1(1ULL), table(new unsigned long long int[1])
  {
    table[0] = 1ULL;
  }

  Choose_Table::Choose_Table (int Npower_ref)
    : Npower(Npower_ref), Npowerp1(Npower_ref + 1ULL)
  {
    dim = Npowerp1 * Npowerp1;

    // We can only go up to ULL_MAX:
    if (log(DP(ULLONG_MAX)) < DP(Npowerp1) * log(2.0))
      JSCerror("Choose_Table:  too large to contruct.");

    table = new unsigned long long int[dim];

    (*this).Fill_table();
  }

  void Choose_Table::Fill_table()
  {
    table[0] = 1ULL;
    int n,m;
    for (n = 0; n <= Npower; ++n) {
      table[Npowerp1 * n] = 1ULL;
      for (m = 1; m < n; ++m) {
	table[Npowerp1 * n + m] = table[Npowerp1 * (n-1) + m - 1] + table[Npowerp1 * (n-1) + m];
      }
      table[Npowerp1 * n + n] = 1ULL;
      for (m = n+1; m <= Npower; ++m)
	table[Npowerp1 * n + m] = 0ULL;
    }
  }

  int Choose_Table::power()
  {
    return(Npower);
  }

  unsigned long long int Choose_Table::choose(int N, int M)
  {
    if (N < 0 || N > Npower) JSCerror("N out of bounds in choose(N,M).");
    if (M < 0 || M > Npower) JSCerror("M out of bounds in choose(N,M).");

    return(table[Npowerp1 * N + M]);
  }

  std::ostream& operator<< (std::ostream& s, Choose_Table& Ref_table)
    {
      s << endl;
      for (int n = 0; n <= Ref_table.power(); ++n) {
	for (int m = 0; m <= Ref_table.power(); ++m)
	  s << Ref_table.choose(n, m) << " ";
	s << endl;
      }
      s << endl;
      return(s);
    }

  Choose_Table::~Choose_Table()
  {
    delete[] table;
  }


} // namespace JSC