MTToolBox  0.2.10
AlgorithmEquidistribution.hpp
[詳解]
1 #ifndef MTTOOLBOX_ALGORITHM_EQUIDISTRIBUTION_HPP
2 #define MTTOOLBOX_ALGORITHM_EQUIDISTRIBUTION_HPP
3 
35 #if __cplusplus >= 201103L
36 #include <memory>
37 #else
38 #if defined(__clang__)
39 #if __has_include(<tr1/memory>)
40 #define MTTOOLBOX_USE_TR1
41 #include <tr1/memory>
42 #else
43 #include <memory>
44 #endif // __has_indlude
45 #else // not clang
46 #define MTTOOLBOX_USE_TR1
47 #include <tr1/memory>
48 #endif // clang
49 #endif // cplusplus version
50 #include <stdexcept>
52 #include <MTToolBox/period.hpp>
53 #include <MTToolBox/util.hpp>
54 
55 namespace MTToolBox {
56 #if defined(MTTOOLBOX_USE_TR1)
57  using std::tr1::shared_ptr;
58 #else
59  using std::shared_ptr;
60 #endif
61 
87  template<typename U, typename V = U>
89  public:
90 
102 
116  linear_generator_vector<U, V>(const ECGenerator& generator) {
117  shared_ptr<ECGenerator> r(generator.clone());
118  rand = r;
119  //rand->seed(1);
120  count = 0;
121  zero = false;
122  setZero(next);
123  }
124 
147  linear_generator_vector<U, V>(const ECGenerator& generator,
148  int bit_pos) {
149  shared_ptr<ECGenerator> r(generator.clone());
150  rand = r;
151  rand->setZero();
152  count = 0;
153  zero = false;
154  next = getOne<U>() << (bit_size<U>() - bit_pos - 1);
155 #if defined(DEBUG)
156  cout << "DEBUG:" << dec << bit_pos << ":"
157  << hex << next << endl;
158 #endif
159  }
160 
161  void add(const linear_generator_vector<U, V>& src);
162  void next_state(int bit_len);
163  void debug_print();
164 
174  shared_ptr<ECGenerator> rand;
186  int count;
187 
197  bool zero;
198 
211  U next;
212  };
213 
237  template<typename U, typename V = U>
239 
250 
251  /*
252  *\japanese
253  * 均等分布次元計算可能な疑似乱数生成器
254  *\endjapanese
255  *\english
256  * Pseudo random number generator which can calculate dimension
257  * of equi-distribution.
258  *\endenglish
259  */
261 
262  public:
263 
296  AlgorithmEquidistribution(const ECGenerator& rand, int bit_length,
297  int mexp = 0) {
298  bit_len = bit_length;
299  size = bit_len + 1;
300  basis = new linear_vec * [size];
301  if (mexp == 0) {
302  stateBitSize = rand.bitSize();
303  } else {
304  stateBitSize = mexp;
305  }
306  for (int i = 0; i < bit_len; i++) {
307  basis[i] = new linear_vec(rand, i);
308  }
309  basis[bit_len] = new linear_vec(rand);
310  basis[bit_len]->next_state(bit_len);
311  }
312 
322  for (int i = 0; i < size; i++) {
323  delete basis[i];
324  }
325  delete[] basis;
326  }
327 
328  int get_all_equidist(int veq[]);
329  int get_equidist(int *sum_equidist);
330  private:
331  int get_equidist_main(int bit_len);
332  void adjust(int new_len);
333 
350  linear_vec **basis;
351 
360  int bit_len;
361 
371  int stateBitSize;
372 
382  int size;
383  };
384 
398  template<typename U, typename V>
399  void AlgorithmEquidistribution<U, V>::adjust(int new_len) {
400  using namespace std;
401  U tmp;
402  setZero(tmp);
403  U mask = ~tmp << (bit_size<U>() - new_len);
404  for (int i = 0; i < size; i++) {
405  basis[i]->next = basis[i]->next & mask;
406  if (isZero(basis[i]->next)) {
407  basis[i]->next_state(new_len);
408  }
409  }
410  }
411 
412 #if defined(DEBUG)
413 
421  template<typename U, typename V>
423  using namespace std;
424 
425  cout << "debug ====" << endl;
426  cout << "count = " << dec << count << endl;
427  cout << "zero = " << zero << endl;
428  cout << "next = " << hex << next << endl;
429  cout << "debug ====" << endl;
430  }
431 #else
432  template<typename U, typename V>
434  }
435 #endif
436 
469  template<typename U, typename V>
471  using namespace std;
472 
473  int sum = 0;
474 
475  veq[bit_len - 1] = get_equidist_main(bit_len);
476 #if defined(DEBUG)
477  for (int i = 0; i < size; i++) {
478  basis[i]->debug_print();
479  }
480 #endif
481  sum += stateBitSize / bit_len - veq[bit_len - 1];
482  bit_len--;
483  for (; bit_len >= 1; bit_len--) {
484  adjust(bit_len);
485  veq[bit_len - 1] = get_equidist_main(bit_len);
486  sum += stateBitSize / bit_len - veq[bit_len - 1];
487  }
488  return sum;
489  }
490 
523  template<typename U, typename V>
525  using namespace std;
526 
527  int veq = get_equidist_main(bit_len);
528  int sum = 0;
529  bit_len--;
530  for (; bit_len >= 1; bit_len--) {
531  adjust(bit_len);
532  sum += stateBitSize / bit_len - get_equidist_main(bit_len);
533  }
534  *sum_equidist = sum;
535  return veq;
536  }
537 
549  template<typename U, typename V>
551  const linear_generator_vector<U, V>& src) {
552  using namespace std;
553 
554  rand->add(*src.rand);
555  next ^= src.next;
556  }
557 
575  template<typename U, typename V>
577  using namespace std;
578 
579  if (zero) {
580  return;
581  }
582  int zero_count = 0;
583  next = rand->generate(bit_len);
584  count++;
585  while (isZero(next)) {
586  zero_count++;
587  if (zero_count > rand->bitSize() * 2) {
588  zero = true;
589  if (rand->isZero()) {
590  zero = true;
591  }
592  break;
593  }
594  next = rand->generate(bit_len);
595  count++;
596  }
597  }
598 
615  template<typename U, typename V>
617  using namespace std;
618  using namespace NTL;
619  int bit_len = v;
620  int pivot_index;
621  int old_pivot = 0;
622 
623  pivot_index = calc_1pos(basis[bit_len]->next);
624  while (!basis[bit_len]->zero) {
625 #if defined(DEBUG)
626  if (pivot_index == -1) {
627  cout << "pivot_index = " << dec << pivot_index << endl;
628  cout << "zero = " << basis[bit_len]->zero << endl;
629  cout << "next = " << hex << basis[bit_len]->next << endl;
630  throw new std::logic_error("pivot error 0");
631  }
632  if (pivot_index >= bit_len) {
633  cout << "pivot_index = " << dec << pivot_index << endl;
634  cout << "bit_len = " << bit_len << endl;
635  cout << "next = " << hex << basis[bit_len]->next << endl;
636  throw new std::logic_error("pivot error 0.1");
637  }
638  if (pivot_index != calc_1pos(basis[pivot_index]->next)) {
639  cerr << "pivot error 1" << endl;
640  cerr << "pivot_index:" << dec << pivot_index << endl;
641  cerr << "calc_1pos:" << dec
642  << calc_1pos(basis[pivot_index]->next) << endl;
643  cerr << "next:" << hex << basis[pivot_index]->next << endl;
644  throw new std::logic_error("pivot error 1");
645  }
646 #endif
647  // アルゴリズムとして、全部のcount を平均的に大きくしたい。
648  // 従って count の小さい方を変化させたい
649  if (basis[bit_len]->count > basis[pivot_index]->count) {
650  swap(basis[bit_len], basis[pivot_index]);
651  }
652 #if defined(DEBUG)
653  cout << "before add bit_len next = " << hex
654  << basis[bit_len]->next << " count = " << dec
655  << basis[bit_len]->count << endl;
656  cout << "before add pivot next = " << hex
657  << basis[pivot_index]->next << " count = " << dec
658  << basis[pivot_index]->count << endl;
659 #endif
660  basis[bit_len]->add(*basis[pivot_index]);
661 #if defined(DEBUG)
662  cout << "after add bit_len next = " << hex
663  << basis[bit_len]->next << " count = " << dec
664  << basis[bit_len]->count << endl;
665 #endif
666  // add の結果 next の最後の1 は必ず 0 になる。
667  // 全部0なら次の状態に進める。(内部でcount が大きくなる)
668  if (isZero(basis[bit_len]->next)) {
669  basis[bit_len]->next_state(bit_len);
670  pivot_index = calc_1pos(basis[bit_len]->next);
671 #if defined(DEBUG)
672  cout << "zero" << endl;
673  cout << "pivot_index = " << dec << pivot_index << endl;
674  if (pivot_index >= bit_len) {
675  cout << "pivot_index = " << dec << pivot_index << endl;
676  cout << "bit_len = " << bit_len << endl;
677  cout << "next = " << hex << basis[bit_len]->next << endl;
678  throw new std::logic_error("pivot error 1.1");
679  }
680  if (basis[bit_len]->zero) {
681  cout << "loop exit condition" << endl;
682  } else if (pivot_index == -1) {
683  cerr << "pivot error 1.1" << endl;
684  throw new std::logic_error("pivot error 1.2");
685  }
686 #endif
687  // 全部0でなければ、pivot_index は小さくなる。
688  // pivot_index が 0 になれば最上位bit のみ1なので
689  // 次の add で全部0になる。
690  } else {
691  old_pivot = pivot_index;
692  pivot_index = calc_1pos(basis[bit_len]->next);
693  if (pivot_index >= bit_len) {
694  cout << "pivot_index = " << dec << pivot_index << endl;
695  cout << "bit_len = " << bit_len << endl;
696  cout << "next = " << hex << basis[bit_len]->next << endl;
697  throw new std::logic_error("pivot error 2");
698  }
699  if (old_pivot <= pivot_index) {
700  cerr << "pivot error 2" << endl;
701  cerr << "old_pivot = " << dec << old_pivot << endl;
702  cerr << "pivot_index = " << dec << pivot_index << endl;
703  throw new std::logic_error("pivot error 2.1");
704  }
705  }
706  }
707 
708  // 計算終了したので最長のベクトルを求める。(長いとはcountが少ないこと)
709  int min_count = basis[0]->count;
710  for (int i = 1; i < bit_len; i++) {
711  if (min_count > basis[i]->count) {
712  min_count = basis[i]->count;
713  }
714  }
715 #if defined(DEBUG)
716  cout << "min_count = " << min_count << endl;
717  cout << "stateBitSize = " << stateBitSize << endl;
718  cout << "bit_len = " << bit_len << endl;
719  cout << "theoretical bound " << (stateBitSize / bit_len) << endl;
720 #endif
721  if (min_count > stateBitSize / bit_len) {
722  cerr << basis[0]->rand->getParamString() << endl;
723  for(int i = 0; i < size; i++) {
724  basis[i]->debug_print();
725  }
726  cerr << "min_count = " << min_count << endl;
727  cerr << "stateBitSize = " << stateBitSize << endl;
728  cerr << "bit_len = " << bit_len << endl;
729  cerr << "over theoretical bound " << (stateBitSize / bit_len)
730  << endl;
731  throw new std::logic_error("over theoretical bound");
732  }
733  return min_count;
734  }
735 }
736 #if defined(MTTOOLBOX_USE_TR1)
737 #undef MTTOOLBOX_USE_TR1
738 #endif
739 #endif // MTTOOLBOX_ALGORITHM_EQUIDISTRIBUTION_HPP
AlgorithmEquidistribution(const ECGenerator &rand, int bit_length, int mexp=0)
コンストラクタ
Definition: AlgorithmEquidistribution.hpp:296
GF(2)ベクトルとしてのGF(2)疑似乱数生成器
Definition: AlgorithmEquidistribution.hpp:88
int count
next_state() が呼ばれた回数 これは多項式としてみた場合の次数に関係がある。
Definition: AlgorithmEquidistribution.hpp:186
virtual int bitSize() const =0
内部状態空間のビットサイズを返す。
void debug_print()
Definition: AlgorithmEquidistribution.hpp:433
~AlgorithmEquidistribution()
デストラクタ
Definition: AlgorithmEquidistribution.hpp:321
void add(const linear_generator_vector< U, V > &src)
ベクトルの加法
Definition: AlgorithmEquidistribution.hpp:550
このクラスはGF(2)線形疑似乱数生成器の均等分布次元を計算するためのクラスである。
Definition: EquidistributionCalculatable.hpp:50
void next_state(int bit_len)
疑似乱数生成器の状態遷移
Definition: AlgorithmEquidistribution.hpp:576
疑似乱数生成器の均等分布次元を計算する
Definition: AlgorithmEquidistribution.hpp:238
bool zero
ゼロベクトルであるかどうかを示す。
Definition: AlgorithmEquidistribution.hpp:197
ユーティリティ関数群
shared_ptr< ECGenerator > rand
GF(2)線形疑似乱数生成器
Definition: AlgorithmEquidistribution.hpp:174
int get_all_equidist(int veq[])
vビット精度の均等分布次元を計算する。 v = bit_len から 1までの均等分布次元を計算して、veq[] に入れる...
Definition: AlgorithmEquidistribution.hpp:470
このクラスはGF(2)線形疑似乱数生成器の均等分布次元を計算するためのクラスである。
U next
疑似乱数生成器の最新の出力(の上位vビット)または 多項式の最高次の係数のなすベクトル ...
Definition: AlgorithmEquidistribution.hpp:211
void setZero(U &x)
ゼロをセットする SIMD型は、そのSIMD型のファイルでこのテンプレートを特殊化する。
Definition: util.hpp:586
static int calc_1pos(uint16_t x)
入力をビット列とみなして最上位の1の位置を0とした最も右側の(下位の)1の位置を返す。 ...
Definition: util.hpp:279
最小多項式の計算と原始性の判定
int get_equidist(int *sum_equidist)
vビット精度の均等分布次元を計算する。
Definition: AlgorithmEquidistribution.hpp:524
EquidistributionCalculatable< U, V > ECGenerator
均等分布次元計算可能な疑似乱数生成器
Definition: AlgorithmEquidistribution.hpp:101
bool isZero(U x)
ゼロかどうか、判定する。 SIMD型は、そのSIMD型のファイルでこのテンプレートを特殊化する。 ...
Definition: util.hpp:691
MTToolBox の名前空間