/****************************/
/* Fisher の直接確率 */
/* coded by Y.Suganuma */
/****************************/
#include <stdio.h>
#include <math.h>
/**********************/
/* n の階乗(対数計算) */
/**********************/
double kaijo(int n)
{
double s = 0.0;
for (int i1 = 2; i1 <= n; i1++)
s += log((double)i1);
return s;
}
/*******************************************/
/* 特定の状態における確率の分母(対数計算） */
/* n_row,n_col : 行と列数 */
/* d_sum : 総和の階乗(対数) */
/*******************************************/
double cal(int **mx, int n_row, int n_col, double d_sum)
{
double s = d_sum;
int i1, i2;
for (i1 = 0; i1 < n_row; i1++) {
for (i2 = 0; i2 < n_col; i2++)
s += kaijo(mx[i1][i2]);
}
return s;
}
/********************************************/
/* 様々な状態における確率を計算 */
/* matrix : 基本状態 */
/* mx : 状態 */
/* n_row,n_col : 行と列数 */
/* d_sum : 総和の階乗(対数) */
/* bunshi,bunbo : 基準行列の分子と分母 */
/* sum : 基準行列の総和 */
/* base : 基準確率 */
/* i : 値を設定する位置 */
/* i/n_col 行 i%n_col 列 */
/********************************************/
double set(int **matrix, int **mx1, int n_row, int n_col, double bunshi,
double bunbo, int sum, double d_sum, double base, int i)
{
double x, y, s = 0.0;
int i1, i2, k, k1, k2, n, n1, n2, sw = 0, ss = 0, **mx2;
// 周辺度数と一致しているかのチェック
for (i1 = 0; i1 < n_row; i1++)
mx1[i1][n_col] = 0;
for (i1 = 0; i1 < n_col; i1++)
mx1[n_row][i1] = 0;
for (i1 = 0; i1 < n_row; i1++) {
for (i2 = 0; i2 < n_col; i2++) {
mx1[i1][n_col] += mx1[i1][i2];
mx1[n_row][i2] += mx1[i1][i2];
ss += mx1[i1][i2];
}
}
for (i1 = 0; i1 < n_row && sw == 0; i1++) {
if (mx1[i1][n_col] != matrix[i1][n_col])
sw = 1;
}
if (sw == 0) {
for (i1 = 0; i1 < n_col && sw == 0; i1++) {
if (mx1[n_row][i1] != matrix[n_row][i1])
sw = 1;
}
}
// 一致している場合
if (sw == 0) {
y = cal(mx1, n_row, n_col, d_sum);
x = exp(bunshi - y);
if (x < base + 1.0e-10)
s = x;
// デバッグ用(計算した組み合わせの出力)
// for (i1 = 0; i1 < n_row; i1++) {
// for (i2 = 0; i2 < n_col; i2++)
// printf(" %d", mx1[i1][i2]);
// printf("\n");
// }
// printf("\n");
}
// 一致していない場合
else {
if (i < n_row*n_col-1 && (cal(mx1,n_row,n_col,d_sum) + kaijo(sum-ss)) > bunbo - 1.0e-10) {
i++;
k1 = i / n_col;
k2 = i % n_col;
n1 = 0;
n2 = 0;
for (i1 = 0; i1 < n_row; i1++)
n1 += mx1[i1][k2];
for (i1 = 0; i1 < n_col; i1++)
n2 += mx1[k1][i1];
n1 = matrix[n_row][k2] - n1;
n2 = matrix[k1][n_col] - n2;
n = (n1 < n2) ? n1 : n2;
k = 0;
mx2 = new int *[n_row+1];
for (i1 = 0; i1 < n_row+1; i1++)
mx2[i1] = new int [n_col+1];
while (k <= n) {
if (k2 < n_col-1 || k2 == n_col-1 && mx1[k1][n_col]+k == matrix[k1][n_col]) {
if (k1 < n_row-1 || k1 == n_row-1 && mx1[n_row][k2]+k == matrix[n_row][k2]) {
for (i1 = 0; i1 < n_row; i1++) {
for (i2 = 0; i2 < n_col; i2++)
mx2[i1][i2] = mx1[i1][i2];
}
mx2[k1][k2] = k;
s += set(matrix, mx2, n_row, n_col, bunshi, bunbo, sum, d_sum, base, i);
}
}
k++;
}
for (i1 = 0; i1 < n_row+1; i1++)
delete [] mx2[i1];
delete [] mx2;
}
}
return s;
}
/********/
/* main */
/********/
int main()
{
double bunshi, bunbo, d_sum, base, s;
int i1, i2, k, n_row, n_col, **mx, **matrix, sum;
// データの入力
// printf("行と列の数は？ ");
scanf("%d %d", &n_row, &n_col);
matrix = new int *[n_row+1];
for (i1 = 0; i1 < n_row+1; i1++)
matrix[i1] = new int [n_col+1];
for (i1 = 0; i1 < n_row; i1++) {
// printf("%d 行目のデータは(スペースで区切って）？ ", i1+1);
for (i2 = 0; i2 < n_col; i2++)
scanf("%d", &matrix[i1][i2]);
}
// 和の計算
for (i1 = 0; i1 < n_row; i1++)
matrix[i1][n_col] = 0;
for (i1 = 0; i1 < n_col; i1++)
matrix[n_row][i1] = 0;
sum = 0;
for (i1 = 0; i1 < n_row; i1++) {
for (i2 = 0; i2 < n_col; i2++) {
sum += matrix[i1][i2];
matrix[i1][n_col] += matrix[i1][i2];
matrix[n_row][i2] += matrix[i1][i2];
}
}
// 生起確率の計算
d_sum = kaijo(sum);
bunshi = 0.0;
for (i1 = 0; i1 < n_row; i1++)
bunshi += kaijo(matrix[i1][n_col]);
for (i1 = 0; i1 < n_col; i1++)
bunshi += kaijo(matrix[n_row][i1]);
bunbo = cal(matrix, n_row, n_col, d_sum);
base = exp(bunshi - bunbo);
printf("生起確率： %f\n", base);
// base より小さい確率の和
mx = new int *[n_row+1];
for (i1 = 0; i1 < n_row+1; i1++)
mx[i1] = new int [n_col+1];
k = 0;
s = 0.0;
while (k <= matrix[0][n_col] && k <= matrix[n_row][0]) {
for (i1 = 0; i1 < n_row; i1++) {
for (i2 = 0; i2 < n_col; i2++)
mx[i1][i2] = 0;
}
mx[0][0] = k;
s += set(matrix, mx, n_row, n_col, bunshi, bunbo, sum, d_sum, base, 0);
k++;
}
printf("確率の和： %f\n", s);
// 領域の開放
for (i1 = 0; i1 < n_row+1; i1++)
delete [] matrix[i1];
delete [] matrix;
for (i1 = 0; i1 < n_row+1; i1++)
delete [] mx[i1];
delete [] mx;
return 0;
}