STAT_Combinatorial.hpp

Go to the documentation of this file.

#ifndef STAT_Combinatorial_HPP
#define STAT_Combinatorial_HPP

#include <limits>
#include <stdexcept>

template<class T>
const tFlag STAT_Combinatorial<T>::COMBINATION=0;

template<class T>
const tFlag STAT_Combinatorial<T>::PERMUTATION=1;

template<class T>
const tFlag STAT_Combinatorial<T>::CIRCULAR_PERMUTATION=2;




// see http://howardhinnant.github.io/combinations.html
// for efficient code

template<class T>
STAT_Combinatorial<T>::STAT_Combinatorial():CORE_Object() {
  setType("STAT_Combinatorial");
  mVector=null;
  mSize=0;
  mCount=0;
  mFunction=COMBINATION;
  mFirstIndex=0;
  mLastIndex=0;
  mDistanceIndex=0;
  mCombinatorialList=null;
  mIsDisplay=false;
}

 
template<class T>
STAT_Combinatorial<T>::~STAT_Combinatorial() {
    if (mVector!=null) delete[] mVector;
}

template<class T>
tBoolean  STAT_Combinatorial<T>::oneComputationDone() {

  
    switch(mFunction) {
    case COMBINATION:
        return oneCombinationDone();
        break;
    case PERMUTATION:
        return permute(mFirstIndex,mLastIndex,mDistanceIndex);
        break;
    case CIRCULAR_PERMUTATION:
        if (mDistanceIndex<=1) return onePermutationDone();
        return permute(mFirstIndex+1,mLastIndex,mDistanceIndex-1);
        break;
    }

  

    return false;

}
template<class T>
void  STAT_Combinatorial<T>::display() {
    cout << mCount<<": ";
    size_t r=0;
    T* iter=mFirstIndex;

    // print the 1st element
    cout << *iter;
    r++;
    // print all the element until last Index
    for (iter++; iter !=mLastIndex; iter++) {
        cout << ", " << *iter;
        r++;
    }

    if (r<mSize) {
        cout << " | ";
        cout << *iter;
        for (iter++; iter !=mEnd; iter++)
        {
            cout << ", " << *iter;
        }
    }
    cout << "\n";

    
    
}

template <class T>
tBoolean STAT_Combinatorial<T>::permute_(T* first1,T* last1,const int& d1) {
    using std::swap;
    switch (d1) {
    case 0:
    case 1:
        return onePermutationDone();
    case 2:
        if (onePermutationDone()) return true;
        
        swap(*first1, *(first1+1));
        return onePermutationDone();
    case 3:
        {
            if (onePermutationDone()) return true;
            
            T *f2 = first1+1;
            T *f3 = f2+1;
            swap(*f2, *f3);
            if (onePermutationDone()) return true;
            
            swap(*first1, *f3);
            swap(*f2, *f3);
            if (onePermutationDone()) return true;
            
            swap(*f2, *f3);
            if (onePermutationDone()) return true;
            
            swap(*first1, *f2);
            swap(*f2, *f3);
            if (onePermutationDone()) return true;
            
            swap(*f2, *f3);
            return onePermutationDone();
        }
    }
    T *fp1 = first1+1;
    T *p;
    for (p = fp1; p != last1; p++) {
        
        if (permute_(fp1, last1, d1-1)) return true;
        std::reverse(fp1, last1);
        swap(*first1, *p);
    }
    return permute_(fp1, last1, d1-1);
}



template<class T>
tBoolean STAT_Combinatorial<T>::permute(T* first1,T* last1,const int& d1) {
    
    
    using std::swap;
    switch (d1) {
    case 0:
    case 1:
        return onePermutationDone();
    case 2:
    {
        if (onePermutationDone()) return true;
        T* i = first1+1;
        swap(*first1, *i);
        if (onePermutationDone()) return true;
        swap(*first1, *i);
    }
    break;
    case 3:
    {
        if (onePermutationDone()) return true;
        T *f2 = first1+1;
        T *f3 = f2+1;
        swap(*f2, *f3);
        if (onePermutationDone()) return true;
        
        swap(*first1, *f3);
        swap(*f2, *f3);
        if (onePermutationDone()) return true;
       
        swap(*f2, *f3);
        if (onePermutationDone()) return true;
        
        swap(*first1, *f2);
        swap(*f2, *f3);
        if (onePermutationDone()) return true;
      
        swap(*f2, *f3);
        if (onePermutationDone()) return true;
       
        swap(*first1, *f3);
    }
    break;
    default:
        T* fp1 = first1+1;
        T* p;
        for (p = fp1; p != last1; p++)
        {
            if (permute_(fp1, last1, d1-1)) return true;
            std::reverse(fp1, last1);
            swap(*first1, *p);
        }
        if (permute_(fp1, last1, d1-1)) return true;
        std::reverse(first1, last1);
        break;
    }
    return false;
}



template<class T>
int  STAT_Combinatorial<T>::computeCombinations(const int& p) {
    ASSERT_IN(p>0);
    ASSERT_IN(p<=(int)mSize);
    mFirstIndex=&mVector[0];
    mLastIndex=&mVector[p];
    mDistanceIndex=p;
    mCount=0;
    mFunction=COMBINATION;
    combine_discontinuous(mFirstIndex, mLastIndex, mDistanceIndex,
                          mLastIndex, mEnd, mSize-p);
    return mCount;
}



template<class T>
int  STAT_Combinatorial<T>::computePermutations(const int& p) {
    ASSERT_IN(p>0);
    ASSERT_IN(p<=(int)mSize);
    mFirstIndex=&mVector[0];
    mLastIndex=&mVector[p];
    mDistanceIndex=p;
    mFunction=PERMUTATION;
    mCount=0;
    combine_discontinuous(mFirstIndex, mLastIndex, mDistanceIndex,
                          mLastIndex, mEnd, mSize-p);
    return mCount;
}
template<class T>
int  STAT_Combinatorial<T>::computeCircularPermutations(const int& p) {
    ASSERT_IN(p>0);
    ASSERT_IN(p<=(int)mSize);
    mFirstIndex=&mVector[0];
    mLastIndex=&mVector[p];
    mDistanceIndex=p;
    mFunction=CIRCULAR_PERMUTATION;
    mCount=0;
    combine_discontinuous(mFirstIndex, mLastIndex, mDistanceIndex,
                          mLastIndex, mEnd, mSize-p);
    return mCount;
}


template<class T>
void STAT_Combinatorial<T>::rotate_discontinuous(T* first1, T* last1,const int&  d1,
                                                 T* first2, T* last2,const int&  d2)
{
    using std::swap;
    if (d1 <= d2)
        std::rotate(first2, std::swap_ranges(first1, last1, first2), last2);
    else
    {
        T* i1 = last1;
        while (first2 != last2)
            swap(*--i1, *--last2);
        std::rotate(first1, i1, last1);
    }
}


template<class T>
bool STAT_Combinatorial<T>::combine_discontinuous(T* first1, T* last1,const int&  d1,
                                                  T* first2, T* last2,const int&  d2,
                                                  const int& d)
{
    
    using std::swap;
    if (d1 == 0 || d2 == 0)
        return oneComputationDone();
    if (d1 == 1)
    {
        T *i2=first2;
        for ( i2= first2; i2 != last2; i2++)
        {
            if (oneComputationDone()) return true;
            swap(*first1, *i2);
        }
    }
    else
    {
        T *f1p = first1+1;
        T *i2 = first2;
        for (int d22 = d2; i2 != last2; ++i2, d22--)
        {
            if (combine_discontinuous(f1p, last1, d1-1, i2, last2, d22,d+1)) return true;
            swap(*first1, *i2);
        }
    }
    if (oneComputationDone()) return true;
    if (d != 0)
        rotate_discontinuous(first1, last1, d1, first2+1, last2, d2-1);
    else
        rotate_discontinuous(first1, last1, d1, first2, last2, d2);
    return false;
}
template <class T>
template <class UINT>
UINT STAT_Combinatorial<T>::getCircularPermutationsNumber(const UINT& p, const UINT& n) {
    UINT d1=p;
    UINT d2=n-p;
    // return d1 > 0 ? (d1+d2)!/(d1*d2!) : 1
    if (d1 == 0)
        return 1;
    UINT r;
    if (d1 <= d2)
    {
        try
        {
            r = getCombinationsNumber(d1, d1+d2);
        }
        catch (const std::overflow_error&)
        {
            throw std::overflow_error("overflow in STAT_Combinatorial::getCircularPermutationsNumber()");
        }
        for (--d1; d1 > 1; --d1)
        {
            if (r > std::numeric_limits<UINT>::max()/d1)
                throw std::overflow_error("overflow in STAT_Combinatorial::getCircularPermutationsNumber()");
            r *= d1;
        }
    }
    else
    {   // functionally equivalent but faster algorithm
        if (d1 > std::numeric_limits<UINT>::max() - d2)
            throw std::overflow_error("overflow in STAT_Combinatorial::getCircularPermutationsNumber()");
        UINT s = d1 + d2;
        r = 1;
        for (; s > d1; --s)
        {
            if (r > std::numeric_limits<UINT>::max()/s)
                throw std::overflow_error("overflow in STAT_Combinatorial::getCircularPermutationsNumber()");
            r *= s;
        }
        for (--s; s > d2; --s)
        {
            if (r > std::numeric_limits<UINT>::max()/s)
                throw std::overflow_error("overflow in STAT_Combinatorial::getCircularPermutationsNumber()");
            r *= s;
        }
    }
    return r;
}
template <class T>
template <class UINT>
UINT STAT_Combinatorial<T>::getPermutationsNumber(const UINT& p, const UINT& n) {
    UINT d1=p;
    UINT d2=n-p;
    // return (d1+d2)!/d2!
    if (d1 > std::numeric_limits<UINT>::max() - d2)
        throw std::overflow_error("overflow in STAT_Combinatorial::getPermutationsNumber");
    UINT s = d1 + d2;
    UINT r = 1;
    for (; s > d2; --s)
    {
        if (r > std::numeric_limits<UINT>::max() / s)
            throw std::overflow_error("overflow in STAT_Combinatorial::getPermutationsNumber");
        r *= s;
    }
    return r;
}

template<class T>
template<class UINT>
UINT STAT_Combinatorial<T>::getCombinationsNumber(const UINT& p,const UINT& n) {
    UINT d1=p;
    UINT d2=n-p;

    // return (d1+d2)!/(d1!*d2!)
    if (d2 < d1)
        std::swap(d1, d2);
    if (d1 == 0)
        return 1;
    if (d1 > std::numeric_limits<UINT>::max() - d2)
        throw std::overflow_error("overflow in STAT_Combinatorial::getCombinationsNumber");
    UINT s = d1 + d2;
    UINT r = s;
    --s;
    UINT k,g,t;
    for (k = 2; k <= d1; ++k, --s) {
        // r = r * n / k, known to not not have truncation error
        g = greatestCommonDivisor(r, k);
        r /= g;
        t = s / (k / g);
        if (r > std::numeric_limits<UINT>::max() / t)
            throw std::overflow_error("overflow in STAT_Combinatorial::getCombinationsNumber");
        r *= t;
    }
    return r;

}
template<class T>
template<class UINT>
UINT STAT_Combinatorial<T>::greatestCommonDivisor(const UINT& p,
                                                 const UINT& n) {
    
    UINT x=p,y=n,t=0;
    while (y != 0) {
        t = x % y;
        x = y;
        y = t;
    }
    return x;
}

#endif