gomu/ext/permutations/permutations.hpp

/**
 * This file is part of Gomu.
 *
 *  Copyright 2016 by Jean Fromentin <jean.fromentin@math.cnrs.fr>
 *
 * Gomu is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * Gomu is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with Gomu. If not, see <http://www.gnu.org/licenses/>. 
 */


#ifndef PERMUTATIONS_HPP
#define PERMUTATIONS_HPP

#include <cassert>
#include <cstdint>
#include <iostream>
#include "smmintrin.h"
#include "coxeter.hpp"

using namespace std;

typedef int8_t int8;
typedef uint32_t uint32;
typedef uint64_t uint64;
typedef int64_t int64;

static const __m128i id128=_mm_setr_epi8(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15);
static const __m128i shmask=_mm_setr_epi8(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,15);

void swap(int8& a,int8& b);
uint64 fact(uint64 n);

template<CoxeterType T> class Permutation{
public:
  uint64 n;
  static Permutation Ptemp;
  union{
    __m128i m128;
    int8 win[16];
  };
  Permutation(uint64=0);
  void setInverse(Permutation&) const;
  const Permutation&  getInverse() const;
  uint64 descentInteger() const;
  int cmp(const Permutation&) const;
  string display() const;
  template<CoxeterType U> friend ostream& operator<<(ostream&,const Permutation<U>&);
};


template<CoxeterType T> class PermutationEnumerator;

  
template<> class PermutationEnumerator<CoxeterA>{
protected:
  Permutation<CoxeterA> sigma;
  uint64 n;
public:
  PermutationEnumerator(uint64);
  void reset();
  bool next();
  size_t size() const;
  const Permutation<CoxeterA>& get() const;
  int cmp(const PermutationEnumerator<CoxeterA>& E) const;
  string display() const;
  friend ostream& operator<<(ostream&,const PermutationEnumerator<CoxeterA>&);
};

template<> class PermutationEnumerator<CoxeterB>:public PermutationEnumerator<CoxeterA>{
protected:
  Permutation<CoxeterB> sigmaB;
  uint64 sign;
public:
  PermutationEnumerator(uint64);
  void reset();
  bool next();
  size_t size() const;
  const Permutation<CoxeterB>& get() const;
  int cmp(const PermutationEnumerator<CoxeterB>& E) const;
  string display() const;
  friend ostream& operator<<(ostream&,const PermutationEnumerator<CoxeterB>&);
};

template<> class PermutationEnumerator<CoxeterD>:public PermutationEnumerator<CoxeterA>{
protected:
  Permutation<CoxeterD> sigmaD;
  uint64 sign;
public:
  PermutationEnumerator(uint64);
  void reset();
  bool next();
  size_t size() const;
  const Permutation<CoxeterD>& get() const;
  int cmp(const PermutationEnumerator<CoxeterD>& E) const;
  string display() const;
  friend ostream& operator<<(ostream&,const PermutationEnumerator<CoxeterD>&);
};


typedef Permutation<CoxeterA> PermutationA;
typedef Permutation<CoxeterB> PermutationB;
typedef Permutation<CoxeterD> PermutationD;

inline void
swap(int8& a,int8& b){
  int8 t=a;
  a=b;
  b=t;
}

inline uint64
fact(uint64 n){
  if(n<2) return 1;
  return n*fact(n-1);
}

//***************
//* Permutation *
//***************

template<CoxeterType T> Permutation<T> Permutation<T>::Ptemp;

template<CoxeterType T> inline
Permutation<T>::Permutation(uint64 r):n(r),m128(id128){
}

template<CoxeterType T> void
Permutation<T>::setInverse(Permutation<T>& inv) const{
  inv.n=n;
  for(uint64 i=1;i<=n;++i){
    int8 v=win[i];
    if(v>0) inv.win[v]=i;
    else inv.win[-v]=-i;
  }
}

template<CoxeterType T> const Permutation<T>&
Permutation<T>::getInverse() const{
  setInverse(Ptemp);
  return Ptemp;
}

template<> inline uint64
Permutation<CoxeterA>::descentInteger() const{
  __m128i temp=_mm_shuffle_epi8(m128,shmask);
  temp=_mm_cmpgt_epi8(m128,temp);
  return _mm_movemask_epi8(temp);
}

template<> inline uint64
Permutation<CoxeterB>::descentInteger() const{
  __m128i temp=_mm_shuffle_epi8(m128,shmask);
  temp=_mm_cmpgt_epi8(m128,temp);
  return _mm_movemask_epi8(temp);
}

template<> inline uint64
Permutation<CoxeterD>::descentInteger() const{
  __m128i temp=_mm_shuffle_epi8(m128,shmask);
  temp=_mm_cmpgt_epi8(m128,temp);
  uint64 res=_mm_movemask_epi8(temp);
  if(win[1]+win[2]<0) res|=1L;
  else res&=(~1L);
  return res;
}

template<CoxeterType T> int
Permutation<T>::cmp(const Permutation& P) const{
  if(n!=P.n){
    if(n<P.n) return 1;
    return -1;
  }
  for(uint64 i=1;i<=n;++i){
    if(win[i]!=P.win[i]){
      if(win[i]<P.win[i]) return 1;
      return -1;
    }
  }
  return 0;
}


template<CoxeterType T> string
Permutation<T>::display() const{
  string str;
  if(n==0) return "()";
  str="(\033[34m"+to_string((int64)win[1])+"\033[0m";
  for(uint64 i=2;i<=n;++i){
    str+=",\033[34m"+to_string((int64)win[i])+"\033[0m";
  }
  str+=')';
  return str;
}

template<CoxeterType T> inline ostream&
operator<<(ostream& os,const Permutation<T>& sigma){
  return os<<sigma.display();
}

//
// PermutationEnumerator
//

inline size_t
PermutationEnumerator<CoxeterA>::size() const{
  return fact(n);
}

inline size_t
PermutationEnumerator<CoxeterB>::size() const{
  return (1<<n)*fact(n);
}

inline size_t
PermutationEnumerator<CoxeterD>::size() const{
  if(n==0) return 1;
  return (1<<(n-1))*fact(n);
}

inline string
PermutationEnumerator<CoxeterA>::display() const{
  return "Enumerator for permutations of type A and rank <= "+to_string(n);
}

inline string
PermutationEnumerator<CoxeterB>::display() const{
  return "Enumerator for permutations of type B and rank <= "+to_string(n);
}

inline string
PermutationEnumerator<CoxeterD>::display() const{
  return "Enumerator for permutations of type D and rank <= "+to_string(n);
}

template<CoxeterType T> inline ostream&
operator<<(ostream& os,const PermutationEnumerator<T>& E){
  return os<<E.display();
}


#endif