#include <cassert>
#include <iomanip>
#include <iostream>
#include <limits>
#include <list>

using namespace std;

typedef list<double> State_Description;

static State_Description operator+(State_Description lhs, const double &rhs) {
  lhs.push_back(rhs);
  return lhs;
}

static State_Description & operator+=(State_Description &lhs, const double &rhs) {
  lhs.push_back(rhs);
  return lhs;
}

static ostream & operator<<(ostream &os, const State_Description &sd) {
  os << '<';
  if(!sd.empty()) {
    auto st = sd.begin();
    os << *st;
    while(++st != sd.end())
      os << ", " << *st;
  }
  os << '>';
  return os;
}

typedef list<pair<State_Description, double>> SD_Values;

static SD_Values operator+(SD_Values lhs, const pair<State_Description, double> &rhs) {
  lhs.push_back(rhs);
  return lhs;
}

static SD_Values & operator+=(SD_Values &lhs, const pair<State_Description, double> &rhs) {
  lhs.push_back(rhs);
  return lhs;
}

static double sum(const SD_Values &lhs) {
  double sum = 0;
  for(auto lt = lhs.begin(); lt != lhs.end(); ++lt)
    sum += lt->second;
  return sum;
}

static double operator*(const SD_Values &lhs, const SD_Values &rhs) {
  double sum = 0;
  assert(lhs.size() == rhs.size());
  auto lt = lhs.begin();
  auto rt = rhs.begin();
  while(lt != lhs.end() && rt != rhs.end()) {
    assert(lt->first == rt->first);
    sum += lt++->second * rt++->second;
  }
  return sum;
}

static ostream & operator<<(ostream &os, const SD_Values &sdv) {
  os << '{';
  if(!sdv.empty()) {
    auto st = sdv.begin();
    os << st->first << ':' << st->second;
    while(++st != sdv.end())
      os << ", " << st->first << ':' << st->second;
  }
  os << '}';
  return os;
}

static double factorial(const int64_t n, const double accumulator = 1) {
  if(n < 2)
    return accumulator;
  return factorial(n - 1, accumulator * n);
}

static double combination(const int64_t n, const int64_t r) {
  return factorial(n) / (factorial(n - r) * factorial(r));
}

static double permutation(const int64_t n, const int64_t r) {
  return factorial(n) / factorial(n - r);
}

static void bw_gen_states(SD_Values &sdv, State_Description sd, const double &multiplier, const int64_t num_blocks, int64_t largest_stack_size, int64_t num_of_largest) {
  sd += largest_stack_size;

  if(largest_stack_size == 1) {
    sdv += make_pair(sd, multiplier);
    return;
  }

  const double stack_0 = multiplier * permutation(num_blocks, largest_stack_size) / num_of_largest++;
  const int64_t remaining_blocks = num_blocks - largest_stack_size;

  if(!remaining_blocks) {
    sdv += make_pair(sd, stack_0);
    return;
  }

  if(remaining_blocks < largest_stack_size) {
    largest_stack_size = remaining_blocks;
    num_of_largest = 1;
  }

  do {
    bw_gen_states(sdv, sd, stack_0, remaining_blocks, largest_stack_size--, num_of_largest);
    num_of_largest = 1;
  } while(largest_stack_size != 0);
}

static void bw_gen_states(SD_Values &sdv, const int64_t num_blocks) {
  for(int64_t largest_stack_size = num_blocks; largest_stack_size != 0; --largest_stack_size)
    bw_gen_states(sdv, State_Description(), 1, num_blocks, largest_stack_size, 1);
}

static double bw_num_states(const int64_t num_blocks) {
  SD_Values sdv;

  bw_gen_states(sdv, num_blocks);

  return sum(sdv);
}

static void bw_gen_sap(SD_Values &sdv, State_Description sd, const int64_t num_blocks, int64_t largest_stack_size, const int64_t num_preceding_stacks) {
  sd += largest_stack_size;

  if(largest_stack_size == 1) {
    const int64_t num_stacks = num_preceding_stacks + num_blocks - 1;
    sdv += make_pair(sd, num_stacks * num_stacks - num_blocks);
    return;
  }

  const int64_t remaining_blocks = num_blocks - largest_stack_size;

  if(!remaining_blocks) {
    sdv += make_pair(sd, num_preceding_stacks * num_preceding_stacks);
    return;
  }

  if(remaining_blocks < largest_stack_size)
    largest_stack_size = remaining_blocks;

  do {
    bw_gen_sap(sdv, sd, remaining_blocks, largest_stack_size--, num_preceding_stacks + 1);
  } while(largest_stack_size != 0);
}

static void bw_gen_sap(SD_Values &sdv, const int64_t num_blocks) {
  for(int64_t largest_stack_size = num_blocks; largest_stack_size != 0; --largest_stack_size)
    bw_gen_sap(sdv, State_Description(), num_blocks, largest_stack_size, 1);
}

static double bw_num_sap(const int64_t num_blocks) {
  SD_Values sdv;
  SD_Values sda;

  bw_gen_states(sdv, num_blocks);
  bw_gen_sap(sda, num_blocks);
  
//   cerr << sdv << endl;
//   cerr << sda << endl;
  
  return sdv * sda;
}

static tuple<int64_t, double, int64_t, double> bw_num_states_sap(const int64_t num_blocks) {
  SD_Values sdv;
  SD_Values sda;

  bw_gen_states(sdv, num_blocks);
  bw_gen_sap(sda, num_blocks);
  
//   cerr << sdv << endl;
//   cerr << sda << endl;
  
  return make_tuple(int64_t(sdv.size()), sum(sdv), sum(sda), sdv * sda);
}

class comma_numpunct : public numpunct<char>
{
  protected:
    virtual char do_thousands_sep() const
    {
        return ',';
    }

    virtual string do_grouping() const
    {
        return "\03";
    }
};

int main(int argc, char **argv) {
  // this creates a new locale based on the current application default
  // (which is either the one given on startup, but can be overriden with
  // locale::global) - then extends it with an extra facet that 
  // controls numeric output.
  locale comma_locale(locale(), new comma_numpunct());

  // tell cout to use our new locale.
  cerr.imbue(comma_locale);
  cout.imbue(comma_locale);

  cerr.precision(numeric_limits<double>::max_digits10);
  cout.precision(numeric_limits<double>::max_digits10);

  for(int i = 1; i < argc; ++i) {
    const int64_t num_blocks = atoi(argv[i]);
    auto num_states_sap = bw_num_states_sap(num_blocks);

    cout << "Blocks World with " << num_blocks << " blocks has " << get<0>(num_states_sap) << " NI states, " << get<1>(num_states_sap) << " states, " << get<2>(num_states_sap) << " NI SAPs and " << get<3>(num_states_sap) << " SAPs." << endl;
  }
}