疑似乱数

完全な乱数はアルゴリズムによっては生成できない。なぜなら、内部条件における状態量が有限であるからである。完全な乱数を得ようとするならば、量子力学に頼るしかない。
とはいえ、実用上「乱数」が必要なこともあるわけで、「アルゴリズムによってそれっぽい（質の高い）乱数を生成する」ことが要求されるのだが、いまのところ「質の高い乱数を生成するアルゴリズム」というのは数少ない。
広く知られているものとしては線形合同法があるが、あまり質がよくない。現在知られている質の良い乱数生成法としては「最長周期律法（M系列法）」程度であり、「内部状態を増やせばバレにくい」という感じのアルゴリズムである。このあたりの悩みはWeb検索されたい。

概要

M系列法は、QRコードと同じように群論（いわゆるガロア理論）に基いており、原始多項式を利用している。
読取と排他的論理和と代入だけで済むため、かなり高速である。

実装

C言語で書けば簡単なのだが、Javaのようなお上品な言語で読みやすく書くのはそれなりに面倒臭い。
原始多項式を用いた最長周期法を用いるときに、乱数の質は生成の種となる乱数テーブルの質によるが、このテーブルを用意するところから苦労がある。ビット列のどこかの桁がすべて 0 だったりすると、二進数で見たときにそこのビットは常に 0 となってうまくない。また、一般的には 0.0 から 1.0 までの値を double で返すのが普通だが、シャッフルなどに用いる場合は引数 n （非負の int）を渡して 0 から n-1 までの値がランダムに返ってくるというのが実用上はありがたい。そういった応用面まで考えると、実装は難しい。

参考文献

• 伏見正則「一様乱数の発生法について」（数学セミナー vol.7 No.10、（1979年10月号）、日本評論社）

コード

Java による、原始多項式を用いた最長周期法によるコードである。

/**
 * Luwis & payne
 * M 系列法（最長周期法）による
 * 疑似乱数の生成
 * x^pwp^ + x^pwq = 1
 */

public class RandLP {
    public static void main( String[] args ) {
      _main();
      return;
    }

    private static void _main() {

      System.out.println(Integer.MAX_VALUE);
      for (int ncnt = 0; ncnt < 100; ncnt+=1) {
       //LP(52);
       System.out.println(LP2(52));
      }
      /*
      for (int i = 0; i < 89; i+=1) {
      //  System.out.print(LP(32768)+ ", ");
       System.out.print((int)(LP(1) * Integer.MAX_VALUE) + ", ");

       if ((i % 5) == 0) System.out.println();
      }
      */
      return;
    }

  static final int IP  = 89;
  static final int PWQ = 38;
  static final int IQ  = IP - PWQ;
  static final int MAX  = Integer.MAX_VALUE;
  static int j = 0;

   static int m[] = {
       1592226944, 1984348928, 1592225792, 2015294208,  424637056,
       2077481472,  387402720, 2083704192,  624210816, 1722843392, //10
       1354712320, 1797627008,  444747744,  163130688, 1824820992,
        681686400,  907007552,  747478976, 1919608448, 2083701376, //20
       1824821504, 1592227968,  747476416, 1984354560, 1493592832,
       2083706368, 1919609856,  999998272, 1132738560,  681684288, //30
       1919607936, 1984353280,  132187704, 1147480320,  907002176,
       1654435200, 2083706752, 1592221696,  163132096, 2015290624, //40
       1919608832,  387402912,   63775904,  945283712, 1919610368,
        841217024,  945289024, 1824818816, 1465799808, 1632732160, //50
        227869760,  681688960, 1202196736,  907006336, 1906826624,
       2147478400, 1240478336, 1654438784,   63774792,  624207360, //60
       1400010240, 1309745920, 1919612288,  837735296, 1418951424,
        747473472, 1702736384, 1240478720,  444746848, 1078530432, //60
       1972103168,   70008336,  322661760, 2077480320, 1972105216,
       2015291904, 1760083584, 1702738816,  907004032,  653889920, //70
       1919610624, 1465797632, 1309751680, 1078532864, 1654435712,
           1523274624, 1760078464,  493044128, 1523277184, 1400006400, //80
          1972103168, 1202194688,  322661760, 2077480320, 1972105216,
           555259072, 1592223488, 1984352384, 1919613696
   };

    public static float LP(int range) {
       int ipx, iqx;
       j       = ( j + 1 )  % IP;
       int k   = ( j + IQ) % IP;
       m[ j ] ^= m[ k ];

//     return ((float)( m[j] * (double)range) / (float )MAX));
      return ((float) m[j]) / ((float)MAX);
    }

    public static long LP2(int range) {
       int ipx, iqx;
       j       = ( j + 1 )  % IP;
       int k   = ( j + IQ) % IP;
       m[ j ] ^= m[ k ];

//     return ((float)( m[j] * (double)range) / (float )MAX));
      return ((long) m[j]) * (long)range / ((long)MAX);
    }
}