monoid.hpp 8.9 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317
  1. #ifndef MONOID_HPP
  2. #define MONOID_HPP
  3. #include <cstdint>
  4. #include <fstream>
  5. // We cant have those as C++ constant because of the #define used below in
  6. // remove_generator manual loop unrolling. I don't know if the unrolling is
  7. // doable using template metaprogamming. I didn't manage to figure out how.
  8. #ifndef MAX_GENUS
  9. #error "Please define the MAX_GENUS macro"
  10. #endif
  11. #define SIZE_BOUND (3*(MAX_GENUS-1))
  12. #define NBLOCKS ((SIZE_BOUND+15) >> 4)
  13. #define SIZE (NBLOCKS << 4)
  14. // const uint_fast8_t MAX_GENUS = 40;
  15. // const uint_fast8_t SIZE_BOUND = (3*(MAX_GENUS-1));
  16. // const uint_fast8_t NBLOCKS = ((SIZE_BOUND+15) >> 4);
  17. // const uint_fast8_t SIZE = (NBLOCKS << 4);
  18. #include <x86intrin.h>
  19. typedef uint8_t epi8 __attribute__ ((vector_size (16)));
  20. typedef uint8_t dec_numbers[SIZE] __attribute__ ((aligned (16)));
  21. typedef epi8 dec_blocks[NBLOCKS];
  22. typedef uint_fast64_t ind_t; // The type used for array indexes
  23. using namespace std;
  24. struct monoid
  25. {
  26. union {
  27. dec_numbers decs;
  28. dec_blocks blocks;
  29. };
  30. // Dont use char as they have to be promoted to 64 bits to do pointer arithmetic.
  31. ind_t conductor, min, genus,left_primitive,left,e,wilf;
  32. };
  33. void init_full_N(monoid &);
  34. void init_ordinary(monoid& O,int m);
  35. void print_monoid(const monoid &);
  36. void print_monoid_gen(const monoid&);
  37. void print_epi8(epi8);
  38. void output(const monoid& m,fstream& f);
  39. inline void copy_blocks( dec_blocks &__restrict__ dst,
  40. const dec_blocks &__restrict__ src);
  41. inline void remove_generator(monoid &__restrict__ dst,
  42. const monoid &__restrict__ src,
  43. ind_t gen,ind_t pos);
  44. inline monoid remove_generator(const monoid &src, ind_t gen,ind_t pos);
  45. // typedef enum { ALL, CHILDREN } generator_type;
  46. class ALL {};
  47. class CHILDREN {};
  48. // template <generator_type T> class generator_iter
  49. template <class T> class generator_iter
  50. {
  51. private:
  52. const monoid &m;
  53. unsigned int mask; // movemask_epi8 returns a 32 bits values
  54. ind_t iblock, gen, bound;
  55. public:
  56. generator_iter(const monoid &mon);
  57. bool move_next();
  58. uint8_t count();
  59. inline ind_t get_gen() const {return gen; };
  60. };
  61. ///////////////////////////////////////////////////////////////////////////
  62. ///////////////////////////////////////////////////////////////////////////
  63. //////////////// Implementation part here for inlining ////////////////
  64. ///////////////////////////////////////////////////////////////////////////
  65. ///////////////////////////////////////////////////////////////////////////
  66. /*
  67. Note: for some reason the code using gcc vector variables is 2-3% faster than
  68. the code using gcc intrinsic instructions.
  69. Here are the results of two runs:
  70. data_mmx = [9.757, 9.774, 9.757, 9.761, 9.811, 9.819, 9.765, 9.888, 9.774, 9.771]
  71. data_epi8 = [9.592, 9.535, 9.657, 9.468, 9.460, 9.679, 9.458, 9.461, 9.629, 9.474]
  72. */
  73. extern inline epi8 load_unaligned_epi8(const uint8_t *t)
  74. { return (epi8) _mm_loadu_si128((__m128i *) (t)); };
  75. extern inline epi8 shuffle_epi8(epi8 b, epi8 sh) // Require SSE 3
  76. { return (epi8) _mm_shuffle_epi8((__m128i) b, (__m128i) sh);}
  77. extern inline int movemask_epi8(epi8 b)
  78. { return _mm_movemask_epi8((__m128i) b);};
  79. const epi8 zero = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
  80. const epi8 block1 = zero + 1;
  81. const uint8_t m1 = 255;
  82. const epi8 shift16[16] =
  83. { { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9,10,11,12,13,14,15},
  84. {m1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9,10,11,12,13,14},
  85. {m1,m1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9,10,11,12,13},
  86. {m1,m1,m1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9,10,11,12},
  87. {m1,m1,m1,m1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9,10,11},
  88. {m1,m1,m1,m1,m1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9,10},
  89. {m1,m1,m1,m1,m1,m1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9},
  90. {m1,m1,m1,m1,m1,m1,m1, 0, 1, 2, 3, 4, 5, 6, 7, 8},
  91. {m1,m1,m1,m1,m1,m1,m1,m1, 0, 1, 2, 3, 4, 5, 6, 7},
  92. {m1,m1,m1,m1,m1,m1,m1,m1,m1, 0, 1, 2, 3, 4, 5, 6},
  93. {m1,m1,m1,m1,m1,m1,m1,m1,m1,m1, 0, 1, 2, 3, 4, 5},
  94. {m1,m1,m1,m1,m1,m1,m1,m1,m1,m1,m1, 0, 1, 2, 3, 4},
  95. {m1,m1,m1,m1,m1,m1,m1,m1,m1,m1,m1,m1, 0, 1, 2, 3},
  96. {m1,m1,m1,m1,m1,m1,m1,m1,m1,m1,m1,m1,m1, 0, 1, 2},
  97. {m1,m1,m1,m1,m1,m1,m1,m1,m1,m1,m1,m1,m1,m1, 0, 1},
  98. {m1,m1,m1,m1,m1,m1,m1,m1,m1,m1,m1,m1,m1,m1,m1, 0} };
  99. const epi8 mask16[16] =
  100. { {m1,m1,m1,m1,m1,m1,m1,m1,m1,m1,m1,m1,m1,m1,m1,m1},
  101. { 0,m1,m1,m1,m1,m1,m1,m1,m1,m1,m1,m1,m1,m1,m1,m1},
  102. { 0, 0,m1,m1,m1,m1,m1,m1,m1,m1,m1,m1,m1,m1,m1,m1},
  103. { 0, 0, 0,m1,m1,m1,m1,m1,m1,m1,m1,m1,m1,m1,m1,m1},
  104. { 0, 0, 0, 0,m1,m1,m1,m1,m1,m1,m1,m1,m1,m1,m1,m1},
  105. { 0, 0, 0, 0, 0,m1,m1,m1,m1,m1,m1,m1,m1,m1,m1,m1},
  106. { 0, 0, 0, 0, 0, 0,m1,m1,m1,m1,m1,m1,m1,m1,m1,m1},
  107. { 0, 0, 0, 0, 0, 0, 0,m1,m1,m1,m1,m1,m1,m1,m1,m1},
  108. { 0, 0, 0, 0, 0, 0, 0, 0,m1,m1,m1,m1,m1,m1,m1,m1},
  109. { 0, 0, 0, 0, 0, 0, 0, 0, 0,m1,m1,m1,m1,m1,m1,m1},
  110. { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,m1,m1,m1,m1,m1,m1},
  111. { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,m1,m1,m1,m1,m1},
  112. { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,m1,m1,m1,m1},
  113. { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,m1,m1,m1},
  114. { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,m1,m1},
  115. { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,m1} };
  116. template<> inline generator_iter<ALL>::generator_iter(const monoid &mon)
  117. : m(mon), bound((mon.conductor+mon.min+15) >> 4)
  118. {
  119. epi8 block;
  120. iblock = 0;
  121. block = m.blocks[0];
  122. mask = movemask_epi8(block == block1);
  123. mask &= 0xFFFE; // 0 is not a generator
  124. gen = - 1;
  125. };
  126. template<> inline generator_iter<CHILDREN>::generator_iter(const monoid &mon)
  127. : m(mon), bound((mon.conductor+mon.min+15) >> 4)
  128. {
  129. epi8 block;
  130. iblock = m.conductor >> 4;
  131. block = m.blocks[iblock] & mask16[m.conductor & 0xF];
  132. mask = movemask_epi8(block == block1);
  133. gen = (iblock << 4) - 1;
  134. };
  135. template <class T> inline uint8_t generator_iter<T>::count()
  136. {
  137. uint8_t nbr = _mm_popcnt_u32(mask); // popcnt returns a 8 bits value
  138. for (ind_t ib = iblock+1; ib < bound; ib++)
  139. nbr += _mm_popcnt_u32(movemask_epi8(m.blocks[ib] == block1));
  140. return nbr;
  141. };
  142. template <class T> inline bool generator_iter<T>::move_next()
  143. {
  144. while (!mask)
  145. {
  146. iblock++;
  147. if (iblock > bound) return false;
  148. gen = (iblock << 4) - 1;
  149. mask = movemask_epi8(m.blocks[iblock] == block1);
  150. }
  151. unsigned char shift = __bsfd (mask) + 1; // Bit Scan Forward
  152. gen += shift;
  153. mask >>= shift;
  154. return true;
  155. };
  156. inline __attribute__((always_inline))
  157. void copy_blocks(dec_blocks &dst, dec_blocks const &src)
  158. {
  159. for (ind_t i=0; i<NBLOCKS; i++) dst[i] = src[i];
  160. }
  161. #include <cassert>
  162. inline __attribute__((always_inline))
  163. void remove_generator(monoid &__restrict__ dst,
  164. const monoid &__restrict__ src,
  165. ind_t gen,
  166. ind_t pos)
  167. {
  168. ind_t start_block, decal;
  169. epi8 block;
  170. assert(src.decs[gen] == 1);
  171. ind_t t=gen+1;
  172. dst.conductor = t;
  173. dst.genus = src.genus + 1;
  174. dst.min=src.min;
  175. int delta=(src.decs[gen+src.min]==2)?0:-1;
  176. dst.e=src.e+delta;
  177. assert(dst.e==((gen==src.min)?src.e+1:((src.decs[gen+src.min]==2)?src.e:src.e-1)));
  178. ind_t k=gen-src.conductor;
  179. assert(dst.conductor==src.conductor+k+1);
  180. dst.left=src.left+k;
  181. dst.left_primitive=src.left_primitive+pos;
  182. //dst.wilf=dst.e*dst.left-dst.conductor;//src.wilf+delta*(src.left+k)-k-1;
  183. dst.wilf=src.wilf+delta*(src.left+k)+(src.e-1)*k-1;
  184. copy_blocks(dst.blocks, src.blocks);
  185. start_block = gen >> 4;
  186. decal = gen & 0xF;
  187. // Shift block by decal uchar
  188. block = shuffle_epi8(src.blocks[0], shift16[decal]);
  189. dst.blocks[start_block] -= ((block != zero) & block1);
  190. #if NBLOCKS >= 5
  191. #define CASE_UNROLL(i_loop) \
  192. case i_loop : \
  193. dst.blocks[i_loop+1] -= (load_unaligned_epi8(srcblock) != zero) & block1; \
  194. srcblock += sizeof(epi8);
  195. {
  196. const uint8_t *srcblock = src.decs + sizeof(epi8) - decal;
  197. switch(start_block)
  198. {
  199. CASE_UNROLL(0);
  200. #if NBLOCKS > 2
  201. CASE_UNROLL(1);
  202. #endif
  203. #if NBLOCKS > 3
  204. CASE_UNROLL(2);
  205. #endif
  206. #if NBLOCKS > 4
  207. CASE_UNROLL(3);
  208. #endif
  209. #if NBLOCKS > 5
  210. CASE_UNROLL(4);
  211. #endif
  212. #if NBLOCKS > 6
  213. CASE_UNROLL(5);
  214. #endif
  215. #if NBLOCKS > 7
  216. CASE_UNROLL(6);
  217. #endif
  218. #if NBLOCKS > 8
  219. CASE_UNROLL(7);
  220. #endif
  221. #if NBLOCKS > 9
  222. CASE_UNROLL(8);
  223. #endif
  224. #if NBLOCKS > 10
  225. CASE_UNROLL(9);
  226. #endif
  227. #if NBLOCKS > 11
  228. CASE_UNROLL(10);
  229. #endif
  230. #if NBLOCKS > 12
  231. CASE_UNROLL(11);
  232. #endif
  233. #if NBLOCKS > 13
  234. CASE_UNROLL(12);
  235. #endif
  236. #if NBLOCKS > 14
  237. CASE_UNROLL(13);
  238. #endif
  239. #if NBLOCKS > 15
  240. CASE_UNROLL(14);
  241. #endif
  242. #if NBLOCKS > 16
  243. CASE_UNROLL(15);
  244. #endif
  245. #if NBLOCKS > 17
  246. CASE_UNROLL(16);
  247. #endif
  248. #if NBLOCKS > 18
  249. CASE_UNROLL(17);
  250. #endif
  251. #if NBLOCKS > 19
  252. #error "Too many blocks"
  253. #endif
  254. }
  255. }
  256. #else
  257. #warning "Loop not unrolled"
  258. for (auto i=start_block+1; i<NBLOCKS; i++)
  259. {
  260. // The following won't work due to some alignment problem (bug in GCC 4.7.1?)
  261. // block = *((epi8*)(src.decs + ((i-start_block)<<4) - decal));
  262. block = load_unaligned_epi8(src.decs + ((i-start_block)<<4) - decal);
  263. dst.blocks[i] -= ((block != zero) & block1);
  264. }
  265. #endif
  266. assert(dst.decs[dst.conductor-1] == 0);
  267. }
  268. inline monoid remove_generator(const monoid &src, ind_t gen,ind_t pos)
  269. {
  270. monoid dst;
  271. remove_generator(dst, src, gen,pos);
  272. return dst;
  273. }
  274. #endif