プログラムのエラーの原因がわかりません

Question

###前提・実現したいこと
次のような問題のプログラムで、コンパイルはできるのですが実行するとエラーが出てしまいます。
（cygwinにgccを入れてやっています）

コインを100回投げる試行を行うことを考える。
メトロポリスの方法を使って、二項分布に従うマルコフ過程から10000個分のサンプル列を作成せよ。
ただし、確率変数kとすると、kの次にはk-1かk+1が1/2の確率でそれぞれ選ばれ、0と100はつながっているものとする。

おそらく除算をしているところで不具合が起きているのではないかと考えているのですが・・・。
お手数ですがよろしくお願いいたします。

###発生している問題・エラーメッセージ

Floating point exception (コアダンプ)

###該当のソースコード

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

#define K 0

uint32_t xor128(void) {             //乱数の生成(0以上2^32未満)
  static uint32_t x = 123456789;
  static uint32_t y = 362436069;
  static uint32_t z = 521288629;
  static uint32_t w = 88675123; 
  uint32_t t;
 
  t = x ^ (x << 11);
  x = y; y = z; z = w;
  return w = (w ^ (w >> 19)) ^ (t ^ (t >> 8)); 
}

int fac(int n) {              //階乗の計算
	int result = 1;
	int k;

	for (k = 1; k <= n; k++) {
		result *= k;
	}
	return result;
}

int cmb(int n, int k){        //combinationの計算
	int c;
	c = fac(n) / (fac(k) * fac(n - k));
	return c;
}

int main(void){
	int P[101];
	int W[101] = {0};
	int i, k, l;
	long a;
	double r;

	P[0] = 1;
	P[100] = 1;

	for (i = 1; i <= 99; i++){
		P[i] = cmb(100,i);
	}

	W[K]++;
	k = K;

	for (i = 1; i < 10000; i++){
		a = xor128();
		if (a < pow(2, 31)){
			l = k - 1;
			if (l == -1){
				l = 100;
			}
		}
		else{
			l = k + 1;
			if (l == 101){
				l = 0;
			}
		}

		a = xor128();
		r = (double)P[l] / (double)P[k];
		r = fmin(r, 1);
		if (a < pow(2, 32)*r){
			k = l;
		}

		W[k]++;
	}

	for (i = 0; i <= 100;i++){
		printf("%d\n",W[i]);
	}

  return 0;
}

Answer 1

途中の計算結果が大きすぎて、intで表現できていないのではないでしょうか?

こちらのサイトによると、100C50は100891344545564193334812497256だそうです。
また、combinationの計算方法も、先に分子と分母を求めて除算結果を返しておられるようですが、ここで必要になる100!はさらに大きく93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000です。
こちらのサイトで紹介されているような、漸化式による方法がいいかもしれません。

Comments

[jackie687456]

無事解決しました。ありがとうございました。

Answer 2

web関係以外は全くわからないのですが、ちらっとint型は扱える数が限られていると聞いたことがあります。
Bongoさんと同じく、intが使えないものと思います。