一様乱数の周期を求めるプログラムについて。

[いちこ]

一様乱数の周期を求めるプログラムを作成しています。
線形合同法を用い、初期値を1として、次に1が出てきたらループを終え、周期を出力したいです。
（漸化式Xn+1 = (a*Xn+c)modM a=1103515245、c=12345、M=2^31 を用います）

プログラムを作成し、実行してみたのですが、ループが最大周期であるはずのMを超えても止まりません。
どこを直せばよいかわからず困っています。
詳しい方がおられましたら、ご教示お願いいたします。

以下に、ソースコード、処理の説明、実行環境を示します。

コード:

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
int main(void){
   double a=1103515245;
   double c=12345;
   double Xn=1; //初期値
   double M=1;
   double n=0;
   int i;

  //最大周期 2^31 を求める
   for(i=31;i>0;i--)
   {
      M=M*2;
   }

   while(1){

      Xn=fmod(a*Xn+c,M); //乱数の計算
      n++;

      if(fmod(n,10000000)==0)
         printf("%.0lf\t Xn=%10.0f\n",n,Xn); //10000000ごとに出力
 
      if(Xn==1) break; //1周期終えたら、ループを抜ける

   }

   // 結果出力
   printf("%.0lf\t Xn=%10.0f\n",n,Xn);
   printf("周期　%.0f\n",n);
   return EXIT_SUCCESS;

}

実行環境
OS Windows7
コンパイラ gcc3.4.4

[h2so5]

doubleだと計算誤差が累積するのでlong longを利用したほうがよいと思います。

コード:

#include<stdio.h>
#include<stdlib.h>
#include<math.h>

int main(void){
  long long a=1103515245;
  long long c=12345;
  long long Xn=1; //初期値
  long long M=pow(2, 31);
  long long n=0;

  while(1){

    Xn=(a*Xn+c) % M; //乱数の計算
    n++;

    if(n % 10000000 == 0)
      printf("%lld\t Xn=%lld\n",n,Xn); //10000000ごとに出力

    if(Xn==1) break; //1周期終えたら、ループを抜ける

  }

  // 結果出力
  printf("%lld\t Xn=%lld\n",n,Xn);
  printf("周期　%lld\n",n);
  return EXIT_SUCCESS;
}

[いちこ]

>> h2so5さん
返信ありがとうございます。
確かに、整数値なのにdoubleはいやだなと思っていました。参考に致します…！

h2so5さんの書いて下さったプログラムを、上と同じ環境で実行してみたのですが、やはり周期Mでループを抜けてくれません。(実行速度はかなり速くなりました！)
なにか良い解決法はないでしょうか。。引き続きよろしくおねがいしますm(_ _)m

[h2so5]

こちらの環境(Apple LLVM version 5.1)ではちゃんと停止します。
そちらの環境はCygwinでしょうか。

[いちこ]

>>h2so5さん
返信ありがとうございます。
以下の2つの環境で実行してみました。
Windowsでは、Cygwinを使っています。
どちらも、停止しませんでした。
Windows7/gcc3.4.4/Cygwin
Ubuntu/gcc4.7.3

ちなみに、コンパイルは
gcc -o sample sample.c -lm
としています。

[nullptr]

h2so5さんに依頼されて動作確認をしました、
Windows 8.1
Visual C++2013 Professional
にて周期の終了を確認。
結構時間がかかったのでしばらく待ったら止まる、とかではないですか？

結果出力です

[いちこ]

>> 新々月さん
返信ありがとうございます。
動作確認もありがとうございます…！
もう一度しばらく待ってみたところ、ちゃんと止まりました。

どうやら、だめだったのはdouble型を使っていたからだったようです。
(以前は明らかにMをどれだけ過ぎても一向に止まらなかったので…)

ちゃんと確認せずに続けてしまってすみません。
h2so5さん、新々月さん、お手数をおかけしてすみませんでした。
ほんとうにありがとうございますm(_ _)m