/***************************************************************************
*  Copyright (c) 2017
*  Joel Berman & Peter Koehler
*
*  Project:       FDn
*  Filename:      a132581.cpp
*  Purpose:       Caclulation Algorithm for A132581
*  Author:        Peter Koehler
*
* - HoC -
*   Date  | Author | Description
* --------+--------+--------------------------------
* 05.07.17 pk       Submission to OEIS
***************************************************************************/

//! \file fdseqownmpi.cpp
//!
//! This program uses a straightforward implementation of the "Splitting result" of
//! J.Berman and P.Koehler, Cardinalities of finite distributive lattices, 
//! Mitteilungen aus dem Mathematischen Seminar Giessen, 121 (1976), 103 - 124.
//! which is implemented in the function f below.
//!
//! The function f is called for the numbers 2^i-1, i=0,...,n, where n is the last
//! index of the sequence to be computed. The numbers 2^i-1 correspond to the posets
//! of the sequence. If i is of the form 2^k for some k, then the poset corresponding
//! to 2^i-1 is the Boolean algebra with i elements. The computed values are printed
//! to Standard Output.
//!
//! As f is called recursively, for speedup the computed values are stored in a map,
//! where they can be looked up efficiently if used in subsequent steps.
//! The number of entries of the map is also printed to Standard Output after the last
//! step.
//!
//! The Multi-precision integer arithmetic is implemented in a separate class
//! BigInt - it would also be possible to use the GNU multiprecion package GMP.

#define _CRT_SECURE_NO_WARNINGS

#include <iostream>
#include <map>
#include <vector>

#include <cstdio>
#include <cstdlib>
#include <cstring>

using namespace std;

//! \class BigInt
//!
//! The class BigInt provides a multiprecision
//! integer arithmetic. Oppeosed to e.g. GNU
//! GMP it uses a (compile time) fixed bitsize, in order
//! to make the operations simple and fast.

class BigInt {
  static const int Dim = 16;                        // 256 Bit - BigEndian
public:
  BigInt(unsigned int n = 0) { memset(nums, 0, sizeof(nums)); nums[Dim - 1] = n & 0xffff; nums[Dim - 2] = (n >> 16) & 0xffff; }
  BigInt(const BigInt & other) { memcpy(nums, other.nums, sizeof(nums)); }

  BigInt operator+ (const BigInt & other);
  BigInt & operator+= (const BigInt & other);
  bool operator!= (int val);
  BigInt operator<< (unsigned int count);
  BigInt & operator>>= (unsigned int count);
  BigInt operator& (const BigInt & other);
  BigInt & operator&= (const BigInt & other);
  BigInt operator ~();
  BigInt operator/ (unsigned short value);
  unsigned short operator% (unsigned short value);
  BigInt & andComp(const BigInt & other);

private:
  unsigned short nums[Dim];

  static void div(const BigInt & first, unsigned short second, BigInt &quotient, unsigned short & remainder);
  static void print(std::ostream & stream, BigInt arg);
  friend bool operator< (const BigInt & first, const BigInt & second);
  friend std::ostream & operator<<(std::ostream &, BigInt);
};

BigInt BigInt::operator+ (const BigInt & other)
{
  BigInt result;
  unsigned short carry = 0;
  for (int n = Dim - 1; n >= 0; n--) {
    unsigned int sum = carry + nums[n] + other.nums[n];
    result.nums[n] = sum & 0xffff;
    carry = sum >= 0x10000 ? 1 : 0;
  }
  return result;
}

BigInt & BigInt::operator+= (const BigInt & other)
{
  unsigned short carry = 0;
  for (int n = Dim - 1; n >= 0; n--) {
    unsigned int sum = carry + nums[n] + other.nums[n];
    nums[n] = sum & 0xffff;
    carry = sum >= 0x10000 ? 1 : 0;
  }
  return *this;
}

bool  BigInt::operator!= (int val)
{
  for (int n = 0; n < Dim - 2; n++) {
    if (nums[n] != 0) return true;
  }
  return (nums[Dim - 1] != (val & 0xffff)) || (nums[Dim - 2] != ((val >> 16) & 0xffff));
}

BigInt BigInt::operator<< (unsigned int count)
{
  BigInt result;
  int offset = count / 16;
  unsigned int shift = count % 16;
  for (int n = 0; n < Dim; n++) {
    unsigned shiftValue = nums[n] << shift;
    if (n - offset >= 0) {
      result.nums[n - offset] |= shiftValue & 0xffff;
      if (n - offset > 0) result.nums[n - offset - 1] |= (shiftValue >> 16) & 0xffff;
    }
  }
  return result;
}

BigInt & BigInt::operator>>= (unsigned int count)
{
  BigInt temp;
  unsigned int offset = count / 16;
  unsigned int shift = count % 16;
  for (int n = 0; n < Dim; n++) {
    unsigned shiftValue = (((unsigned)nums[n]) << 16) >> shift;
    if (n + offset < Dim) {
      temp.nums[n + offset] |= (shiftValue >> 16) & 0xffff;
      if (n + offset < Dim - 1) temp.nums[n + offset + 1] |= (shiftValue & 0xffff);
    }
  }
  *this = temp;
  return *this;
}

BigInt BigInt::operator& (const BigInt & other)
{
  BigInt result;
  for (int n = 0; n < Dim; n++) result.nums[n] = nums[n] & other.nums[n];
  return result;
}

BigInt & BigInt::operator&= (const BigInt & other)
{
  for (int n = 0; n < Dim; n++) nums[n] &= other.nums[n];
  return *this;
}

BigInt BigInt::operator ~()
{
  BigInt result;
  for (int n = 0; n < Dim; n++) result.nums[n] = (~nums[n]) & 0xffff;
  return result;
}


bool operator< (const BigInt & first, const BigInt & second)
{
  int index = 0;
  while (first.nums[index] == second.nums[index] && index < BigInt::Dim) index++;
  return index != BigInt::Dim && first.nums[index] < second.nums[index];
}

BigInt BigInt::operator/ (unsigned short value)
{
  unsigned short remainder;
  BigInt result;
  div(*this, value, result, remainder);
  return result;
}

unsigned short  BigInt::operator% (unsigned short value)
{
  unsigned short remainder;
  BigInt result;
  div(*this, value, result, remainder);
  return remainder;
}

//! \brief And-With-Complement (x &= ~y)
BigInt & BigInt::andComp(const BigInt & other)
{
  for (int n = 0; n < Dim; n++) nums[n] &= ~other.nums[n];
  return *this;
}

void BigInt::div(const BigInt & first, const unsigned short second, BigInt & quotient, unsigned short & remainder)
{
  unsigned int carry = 0;
  for (int n = 0; n < Dim; n++) {
    unsigned int temp = 0x10000 * carry + first.nums[n];
    quotient.nums[n] = temp / second;
    carry = temp % second;
  }
  remainder = carry & 0xffff;
}

void BigInt::print(ostream & stream, BigInt arg)
{
  if (arg < 10) stream << arg.nums[Dim - 1];
  else {
    unsigned short remainder;
    div(arg, 10, arg, remainder);
    stream << arg << remainder;
  }
}
std::ostream & operator<<(std::ostream & stream, BigInt arg)
{
  BigInt::print(stream, arg);
  return stream;
}

typedef std::map<BigInt, BigInt> ValueMap;
static ValueMap valueMap;
static vector<BigInt> powers;

static unsigned int bitsize;

//! \brief Computation Function
//! \param m: Poset expressed as number
//! \result: number of semiideals of the poset
BigInt f ( BigInt m )
{
  register unsigned x, t ;
  BigInt u, s, p;

  // If the value is known, take it from the map
  ValueMap::iterator iter = valueMap.find(m);
  if ( iter != valueMap.end() ) {
    return iter->second;
  }

  u = m;
  // Determine log2(m)
  for (x = 0; u != 0; x++) {
    u >>= 1;
  }
  x -= 1;
    
  // Size of Poset P-{x}
  BigInt p1 = m;
  p1.andComp(powers[x]);
  s = f(p1);

  // Size of Poset P-cone(x)
  p = m;
  for ( t = 0; t <= x; t++ ) {
    if ( (t & x) == t ) {
      p.andComp(powers[t]);
    }
  }

  s += f(p);

  // Store the value
  valueMap[m] = s;

  return s; 
}

//! \brief Main Function
//! \param argc: Number of command line arguments
//! \param args: Argument strings
int main (int argc, char *argv[] )
{
  int x = -1 ;

  if ( ( argc == 1 )                             ||
       ( *argv[1] == '?' )                       ||
       ( 1 != sscanf ( argv [1], "%d", &x ) )    ||
       ( x < 0 )
       ) cerr << "The program calculates the first terms of\n"
                "the integer sequence A132581 of OEIS" << endl;
  else {
    bitsize = x;
    BigInt n(1);
    int i;

    valueMap[0] = BigInt(1);
    valueMap[1] = BigInt(2);
    for(int n=0; n<x; n++) {
      powers.push_back(BigInt(1) << n);
    }
    cout << "# First " << x << " terms of A132581 - computed by Peter Koehler (05.07.2017)" << endl;
    for (i = 1; i <= x; i = i + 1, n = n + n + BigInt(1)) {
      BigInt cur = f(n);
      std::cout <<  i << " " << cur << std::endl;
    }
    cout << "# Map usage: " << valueMap.size() << " entries" << endl;;
  }
  return 0;
}