Paper Implementation : Sequential Monte Carlo Methods for Statistical

Analysis of Tables

Implemented this in C and Julia, lol it took me 2 years to understand the paper — skill issue on my end

I implement the paper Sequential Monte Carlo Methods for Statistical Analysis of Tables by Yuguo Chen, Persi Diaconis, Susan P. Holmes, and Jun S.Liu

These are part of my startup notes —

• Test using contingency table in Diaconis’ 1995 paper RECTANGULAR ARRAYS WITH FIXED MARGINS

5 | 2 | 3 | 10 50| 7| 5 |62 3 | 6 | 4 |13 5 | 3 | 3 | 11

2 | 7 |30| 39

Update — I changed to using libgmp, my accuracy issues were a consequence of integer overflow.

3rd Feb 2025 Update

I’m obsessed with this paper. I saw that only the row sums and column sums are needed. I changed the code (3rd entry with libgmp) to use less memory.

C implementation — rand() is less accurate than Julia’s mersenne twister

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <assert.h>
#include <string.h> // For memcpy
#define MAX_ChenDiaconis(a, b) ((a) > (b) ? (a) : (b))
#define MIN_ChenDiaconis(a, b) ((a) < (b) ? (a) : (b))
#define MATRIX_INDEX(x, y, cols) ((x) * (cols) + (y))
int *GenerateRowSums(int height, int width, int *matrix)
{
 int *result = calloc(height, sizeof(int));
 for(int i = 0; i < height; i++)
 {
  int sum = 0;
  for(int j = 0; j < width; j++)
  {
   sum += matrix[(i*width) + j];
  }
  result[i]=sum;
 }
 return result;
}

int *GenerateColumnSums(int height, int width, int *matrix)
{
 int *result = calloc(width, sizeof(int));
 for(int j = 0; j < width; j++)
 {
  int sum = 0;
  for(int i = 0; i < height; i++)
  {
   sum += matrix[(i*width) + j];
  }
  result[j] = sum;
 }
 return result;
}

void PrintIntArray_ChenDiaconis(int length, int *array)
{
 for(int i = 0; i < length; i++)
 {
  printf("%3d, ", array[i]);
 }
 printf("\n");
}

int RandRange(int lowerBound, int upperBound)
{
    return lowerBound + rand() % (upperBound - lowerBound + 1);
}

int MakeDecision(int currentColumnSum, int currentRowSum, int currentTotalSum, int *currentQ)
{
 int decision = 0;
 int lowerBound = MAX_ChenDiaconis(0, currentColumnSum + currentRowSum - currentTotalSum);
 int upperBound = MIN_ChenDiaconis(currentColumnSum, currentRowSum);
 decision = RandRange(lowerBound, upperBound);
 *currentQ = upperBound - lowerBound + 1;
 return decision;
}

void PrintTable_ChenDiaconis(int numberOfRows, int numberOfColumns, int *table)
{
 for(int i = 0; i < numberOfRows; i++)
 {
  for(int j = 0; j < numberOfColumns; j++)
  {
   printf("%3d, ", table[MATRIX_INDEX(i,j,numberOfColumns)]);
  }
  printf("\n");
 }
}

int GenerateTable(int numberOfTables, int numberOfRows, int numberOfColumns, int tableIndex, int totalSum, int *rowSums, int *columnSums, int **tableHolder)
{
 assert(tableIndex > -1);
 assert(tableIndex < numberOfTables-1);
 long currentQ = 1; //maybe change to mpz_t
 int currentColumnSum = 0;
 int currentRowSum = 0;
 int currentTotalSum = 0;
 int currentTableValue = 0;
 int decision = 0;
 int qHolder = 0;
 for(int i = 0; i < numberOfColumns; i++)
 {
  currentColumnSum = columnSums[i];
  currentTotalSum = totalSum;
  for(int j = 0; j < numberOfRows; j++)
  {
   currentTableValue = tableHolder[tableIndex][MATRIX_INDEX(j,i,numberOfColumns)];
   currentRowSum = rowSums[j];
   decision = MakeDecision(currentColumnSum, currentRowSum, currentTotalSum, &qHolder);
   currentQ *= qHolder;
   tableHolder[tableIndex+1][MATRIX_INDEX(j,i,numberOfColumns)] = decision; 


   currentColumnSum -= decision;
   currentTotalSum -= currentRowSum;
   rowSums[j] -= decision;
   totalSum -= decision;

  }
 }
 return currentQ;
}
int main()
{
 srand(55455);
 int numberOfTables = 1000000;
 int numberOfRows = 5;
 int numberOfColumns = 3;
 long qT = 0;
 int **tableHolder = malloc(numberOfTables * sizeof(int*));for(int i = 0; i < numberOfTables; i++){tableHolder[i] = calloc(numberOfRows * numberOfColumns, sizeof(int));}

 tableHolder[0][MATRIX_INDEX(0,0,numberOfColumns)] = 5; tableHolder[0][MATRIX_INDEX(0,1,numberOfColumns)] = 2;tableHolder[0][MATRIX_INDEX(0,2,numberOfColumns)] = 3;
 tableHolder[0][MATRIX_INDEX(1,0,numberOfColumns)] = 50;tableHolder[0][MATRIX_INDEX(1,1,numberOfColumns)] = 7;tableHolder[0][MATRIX_INDEX(1,2,numberOfColumns)] = 5;
 tableHolder[0][MATRIX_INDEX(2,0,numberOfColumns)] = 3; tableHolder[0][MATRIX_INDEX(2,1,numberOfColumns)] = 6;tableHolder[0][MATRIX_INDEX(2,2,numberOfColumns)] = 4;
 tableHolder[0][MATRIX_INDEX(3,0,numberOfColumns)] = 5; tableHolder[0][MATRIX_INDEX(3,1,numberOfColumns)] = 3;tableHolder[0][MATRIX_INDEX(3,2,numberOfColumns)] = 3;
 tableHolder[0][MATRIX_INDEX(4,0,numberOfColumns)] = 2; tableHolder[0][MATRIX_INDEX(4,1,numberOfColumns)] = 7;tableHolder[0][MATRIX_INDEX(4,2,numberOfColumns)] = 30;
 int *rowSums = GenerateRowSums(numberOfRows, numberOfColumns, tableHolder[0]);
 int *columnSums = GenerateColumnSums(numberOfRows, numberOfColumns, tableHolder[0]);

 int *rowSumsCopy = GenerateRowSums(numberOfRows, numberOfColumns, tableHolder[0]);
 int *columnSumsCopy = GenerateColumnSums(numberOfRows, numberOfColumns, tableHolder[0]);

 int totalSum = 0; for(int i = 0; i < numberOfRows; i++){totalSum += rowSums[i];} 
 //PrintIntArray_ChenDiaconis(numberOfRows, rowSums);PrintIntArray_ChenDiaconis(numberOfColumns, columnSums);printf("%3d\n", totalSum);


 for(int tableIndex = 0; tableIndex < numberOfTables-1; tableIndex++)
 {
  memcpy(rowSumsCopy, rowSums, numberOfRows * sizeof(int));
  memcpy(columnSumsCopy, columnSums, numberOfColumns * sizeof(int));
  qT += GenerateTable(numberOfTables, numberOfRows, numberOfColumns, tableIndex, totalSum, rowSumsCopy, columnSumsCopy, tableHolder);
  //PrintTable_ChenDiaconis(numberOfRows, numberOfColumns, tableHolder[tableIndex+1]);printf("\n");
 }
 printf("Approximate Count : %3ld : %d\n", qT / numberOfTables, 238243776);


 free(rowSums);free(columnSums);free(rowSumsCopy);free(columnSumsCopy);
 for(int i = 0; i < numberOfTables; i++){free(tableHolder[i]);}free(tableHolder);
 return 0;
}

Julia — accurate because of Mersenne Twister

using Random
using HypothesisTests
using Dates

myGlobalProgramRNG = MersenneTwister(123) 

function FindRowSums(table, tableIndex)
 for i in 1:numberOfRows-1
  rowSum = 0
  for j in 1:numberOfColumns-1
   rowSum += table[i, j ,tableIndex]
  end
  table[i, numberOfColumns ,tableIndex] = rowSum
 end
end

function FindColumnSums(table, tableIndex)
 for i in 1:numberOfColumns
  columnSum = 0
  for j in 1:numberOfRows-1
   columnSum += table[j,i,tableIndex]
  end
  table[numberOfRows,i,tableIndex] = columnSum
 end
end

function MakeDecision(c1,r1,M)
 #lowerBound = max(0, c1+r1-M)
 lowerBound = max(0, c1+r1-M) 
 upperBound = min(r1,c1) 
 #println("Lower: ", lowerBound, " Upper: ", upperBound)
 result = rand(myGlobalProgramRNG,lowerBound:upperBound)
 currentQ = ((upperBound - lowerBound) + 1)
 return (result, currentQ)
end

function GenerateTable(table, tableIndex)
 currentQ  = 1
 for i in 1:numberOfColumns-1
  c1 = table[numberOfRows,i,tableIndex] #Column sum
  M  = table[numberOfRows,numberOfColumns,tableIndex]
  for j in 1:numberOfRows-1
   a = table[j,i,tableIndex] 
   r1 = table[j,numberOfColumns,tableIndex]
   #println(c1," ", r1," ",M)
    decisionResult = MakeDecision(c1,r1,M)
   a11 = decisionResult[1]
   currentQ *= decisionResult[2]
   c1 -= a11
   M  -= r1 
   table[j,i,tableIndex+1] = a11
   table[j,numberOfColumns,tableIndex] -= a11
   table[numberOfRows,numberOfColumns,tableIndex] -= a11
  end
 end
 return currentQ
end

function FindMarginalSums(table, tableIndex)
 FindRowSums(table,tableIndex)
 FindColumnSums(table,tableIndex)
end

function GenerateTables(table, numberOfTables)
    t1 = now()
 qT = 0
 for i in 1 : numberOfTables-1
  FindMarginalSums(table,i)
  qt = GenerateTable(table,i) #returns an integer
  FindMarginalSums(table,i)
  qT += qt
  qTArray[i] = qt
 end
 FindMarginalSums(table,numberOfTables)
    t2 =now() + Hour(1)
 println("Approximate Total: ", qT/numberOfTables)
    println("Time", t2 - t1)
end

function Indicator(originalChiStatistic, chiSquareStatistic)
 if(chiSquareStatistic <= originalChiStatistic)
  return 1
 else
  return 0
 end
end
function CalculateMu(table, expectedTable, numberOfTables)
 #degreesOfFreedom = (numberOfRows - 2) * (numberOfColumns - 2)
 originalChiStatistic = 0;

 muNumerator = 0
 muDenominator = 0
 z = 0
 for i in 1 : numberOfTables-1
  chiSquareStatistic = 0
  for j in 1 : numberOfRows-1
   for k in 1 : numberOfColumns-1
    expectedValue = convert(Float64, table[j, numberOfColumns,i]) * convert(Float64, table[numberOfRows, k,i]) / convert(Float64, table[numberOfRows, numberOfColumns,i]) 
    expectedValue = (convert(Float64, table[j, k,i]) - expectedValue)^2 / expectedValue
    chiSquareStatistic += expectedValue
   end
  end
  if i == 1
   originalChiStatistic = chiSquareStatistic
   println(originalChiStatistic)
  end

  if i == 2
   println(chiSquareStatistic)
  end

  fT = Indicator(originalChiStatistic, chiSquareStatistic)
  muNumerator += convert(Float64, fT * (1 / qTArray[i]))
  muDenominator += convert(Float64, (1 / qTArray[i]))
  #z += fT
  #println(chiSquareStatistic, " ", fT," ", qTArray[i])
 end
 mu = muNumerator / muDenominator

 println("Mu value: ", mu, " ") 
end

numberOfTables = 100000
numberOfRows = 6 #6
numberOfColumns = 4 #4
table = zeros(Int64, (numberOfRows,numberOfColumns,numberOfTables))
expectedTable = zeros(Float64, (numberOfRows-1,numberOfColumns-1,numberOfTables))
qTArray = zeros(Int64,numberOfTables, 1)

table[1,1,1] =  5; table[1,2,1] = 2; table[1,3,1] = 3;
table[2,1,1] = 50; table[2,2,1] = 7; table[2,3,1] = 5;
table[3,1,1] =  3; table[3,2,1] = 6; table[3,3,1] = 4;
table[4,1,1] =  5; table[4,2,1] = 3; table[4,3,1] = 3;
table[5,1,1] =  2; table[5,2,1] = 7; table[5,3,1] =30;
#Perhaps try
#table[1,1,1] =  5; table[1,2,1] = 2; table[1,3,1] = 3;
#table[2,1,1] = 50; table[2,2,1] = 7; table[2,3,1] = 5;

GenerateTables(table, numberOfTables)
CalculateMu(table, expectedTable, numberOfTables)

libgmp implementation in C, super accurate

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gmp.h>
#include <assert.h>
#include <string.h> // For memcpy
#define MAX_ChenDiaconis(a, b) ((a) > (b) ? (a) : (b))
#define MIN_ChenDiaconis(a, b) ((a) < (b) ? (a) : (b))
#define MATRIX_INDEX(x, y, cols) ((x) * (cols) + (y))

int *GenerateRowSums(int height, int width, int *matrix)
{
 int *result = calloc(height, sizeof(int));
 for(int i = 0; i < height; i++)
 {
  int sum = 0;
  for(int j = 0; j < width; j++)
  {
   sum += matrix[(i*width) + j];
  }
  result[i]=sum;
 }
 return result;
}

int *GenerateColumnSums(int height, int width, int *matrix)
{
 int *result = calloc(width, sizeof(int));
 for(int j = 0; j < width; j++)
 {
  int sum = 0;
  for(int i = 0; i < height; i++)
  {
   sum += matrix[(i*width) + j];
  }
  result[j] = sum;
 }
 return result;
}

void PrintIntArray_ChenDiaconis(int length, int *array)
{
 for(int i = 0; i < length; i++)
 {
  printf("%3d, ", array[i]);
 }
 printf("\n");
}

int RandRange(int lowerBound, int upperBound)
{
    return lowerBound + rand() % (upperBound - lowerBound + 1);
}

int MakeDecision(int currentColumnSum, int currentRowSum, int currentTotalSum, int *currentQ)
{
 int decision = 0;
 int lowerBound = MAX_ChenDiaconis(0, currentColumnSum + currentRowSum - currentTotalSum);
 int upperBound = MIN_ChenDiaconis(currentColumnSum, currentRowSum);
 decision = RandRange(lowerBound, upperBound);
 *currentQ = upperBound - lowerBound + 1;
 return decision;
}

void PrintTable_ChenDiaconis(int numberOfRows, int numberOfColumns, int *table)
{
 for(int i = 0; i < numberOfRows; i++)
 {
  for(int j = 0; j < numberOfColumns; j++)
  {
   printf("%3d, ", table[MATRIX_INDEX(i,j,numberOfColumns)]);
  }
  printf("\n");
 }
}

void GenerateTable(int numberOfTables, int numberOfRows, int numberOfColumns, int tableIndex, int totalSum, int *rowSums, int *columnSums, int **tableHolder, mpz_t qT)
{
 assert(tableIndex > -1);
 assert(tableIndex < numberOfTables-1);
 mpz_t currentQ; mpz_init(currentQ);mpz_set_ui(currentQ, 1);
 int currentColumnSum = 0;
 int currentRowSum = 0;
 int currentTotalSum = 0;
 int currentTableValue = 0;
 int decision = 0;
 int qHolder = 0;
 for(int i = 0; i < numberOfColumns; i++)
 {
  currentColumnSum = columnSums[i];
  currentTotalSum = totalSum;
  for(int j = 0; j < numberOfRows; j++)
  {
   currentTableValue = tableHolder[tableIndex][MATRIX_INDEX(j,i,numberOfColumns)];
   currentRowSum = rowSums[j];
   decision = MakeDecision(currentColumnSum, currentRowSum, currentTotalSum, &qHolder);
   mpz_mul_ui(currentQ, currentQ, (unsigned long int) qHolder);
   tableHolder[tableIndex+1][MATRIX_INDEX(j,i,numberOfColumns)] = decision; 
   currentColumnSum -= decision;
   currentTotalSum -= currentRowSum;
   rowSums[j] -= decision;
   totalSum -= decision;
  }
 }
 mpz_add(qT, qT, currentQ);
 mpz_clear(currentQ);
}

int TestDiaconis()
{
 srand(55455);
 int numberOfTables = 10000;
 int numberOfRows = 5;
 int numberOfColumns = 3;
 mpz_t qT; mpz_init(qT);mpz_set_ui(qT, 0);
 int **tableHolder = malloc(numberOfTables * sizeof(int*));for(int i = 0; i < numberOfTables; i++){tableHolder[i] = calloc(numberOfRows * numberOfColumns, sizeof(int));}

 tableHolder[0][MATRIX_INDEX(0,0,numberOfColumns)] = 5; tableHolder[0][MATRIX_INDEX(0,1,numberOfColumns)] = 2;tableHolder[0][MATRIX_INDEX(0,2,numberOfColumns)] = 3;
 tableHolder[0][MATRIX_INDEX(1,0,numberOfColumns)] = 50;tableHolder[0][MATRIX_INDEX(1,1,numberOfColumns)] = 7;tableHolder[0][MATRIX_INDEX(1,2,numberOfColumns)] = 5;
 tableHolder[0][MATRIX_INDEX(2,0,numberOfColumns)] = 3; tableHolder[0][MATRIX_INDEX(2,1,numberOfColumns)] = 6;tableHolder[0][MATRIX_INDEX(2,2,numberOfColumns)] = 4;
 tableHolder[0][MATRIX_INDEX(3,0,numberOfColumns)] = 5; tableHolder[0][MATRIX_INDEX(3,1,numberOfColumns)] = 3;tableHolder[0][MATRIX_INDEX(3,2,numberOfColumns)] = 3;
 tableHolder[0][MATRIX_INDEX(4,0,numberOfColumns)] = 2; tableHolder[0][MATRIX_INDEX(4,1,numberOfColumns)] = 7;tableHolder[0][MATRIX_INDEX(4,2,numberOfColumns)] = 30;
 int *rowSums = GenerateRowSums(numberOfRows, numberOfColumns, tableHolder[0]);
 int *columnSums = GenerateColumnSums(numberOfRows, numberOfColumns, tableHolder[0]);

 int *rowSumsCopy = GenerateRowSums(numberOfRows, numberOfColumns, tableHolder[0]);
 int *columnSumsCopy = GenerateColumnSums(numberOfRows, numberOfColumns, tableHolder[0]);

 int totalSum = 0; for(int i = 0; i < numberOfRows; i++){totalSum += rowSums[i];} 
 //PrintIntArray_ChenDiaconis(numberOfRows, rowSums);PrintIntArray_ChenDiaconis(numberOfColumns, columnSums);printf("%3d\n", totalSum);


 for(int tableIndex = 0; tableIndex < numberOfTables-1; tableIndex++)
 {
  memcpy(rowSumsCopy, rowSums, numberOfRows * sizeof(int));
  memcpy(columnSumsCopy, columnSums, numberOfColumns * sizeof(int));
  GenerateTable(numberOfTables, numberOfRows, numberOfColumns, tableIndex, totalSum, rowSumsCopy, columnSumsCopy, tableHolder, qT);
  //PrintTable_ChenDiaconis(numberOfRows, numberOfColumns, tableHolder[tableIndex+1]);printf("\n");
 }
 mpz_div_ui(qT, qT, (unsigned long int) numberOfTables);
 gmp_printf("Approximate Count : %3Zd : %d\n", qT , 238243776);


 mpz_clear(qT);free(rowSums);free(columnSums);free(rowSumsCopy);free(columnSumsCopy);
 for(int i = 0; i < numberOfTables; i++){free(tableHolder[i]);}free(tableHolder);
 return 0;
}

3rd Feb 2025 update

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gmp.h>
#include <assert.h>
#include <string.h> // For memcpy
#define MAX_ChenDiaconis(a, b) ((a) > (b) ? (a) : (b))
#define MIN_ChenDiaconis(a, b) ((a) < (b) ? (a) : (b))
#define MATRIX_INDEX(x, y, cols) ((x) * (cols) + (y))

int *GenerateRowSums(int height, int width, int *matrix)
{
 int *result = calloc(height, sizeof(int));
 for(int i = 0; i < height; i++)
 {
  int sum = 0;
  for(int j = 0; j < width; j++)
  {
   sum += matrix[(i*width) + j];
  }
  result[i]=sum;
 }
 return result;
}

int *GenerateColumnSums(int height, int width, int *matrix)
{
 int *result = calloc(width, sizeof(int));
 for(int j = 0; j < width; j++)
 {
  int sum = 0;
  for(int i = 0; i < height; i++)
  {
   sum += matrix[(i*width) + j];
  }
  result[j] = sum;
 }
 return result;
}

void PrintIntArray_ChenDiaconis(int length, int *array)
{
 for(int i = 0; i < length; i++)
 {
  printf("%3d, ", array[i]);
 }
 printf("\n");
}

int RandRange(int lowerBound, int upperBound)
{
 return lowerBound + rand() % (upperBound - lowerBound + 1);
}

int MakeDecision(int currentColumnSum, int currentRowSum, int currentTotalSum, int *currentQ)
{
 int decision = 0;
 int lowerBound = MAX_ChenDiaconis(0, currentColumnSum + currentRowSum - currentTotalSum);
 int upperBound = MIN_ChenDiaconis(currentColumnSum, currentRowSum);
 decision = RandRange(lowerBound, upperBound);
 *currentQ = upperBound - lowerBound + 1;
 return decision;
}

void PrintTable_ChenDiaconis(int numberOfRows, int numberOfColumns, int *table)
{
 for(int i = 0; i < numberOfRows; i++)
 {
  for(int j = 0; j < numberOfColumns; j++)
  {
   printf("%3d, ", table[MATRIX_INDEX(i,j,numberOfColumns)]);
  }
  printf("\n");
 }
}

void GenerateTable(int numberOfRows, int numberOfColumns, int totalSum, int *rowSums, int *columnSums, int *tableHolder, mpz_t qT)
{
 mpz_t currentQ; mpz_init(currentQ);mpz_set_ui(currentQ, 1);
 int currentColumnSum = 0;
 int currentRowSum = 0;
 int currentTotalSum = 0;
 int currentTableValue = 0;
 int decision = 0;
 int qHolder = 0;
 for(int i = 0; i < numberOfColumns; i++)
 {
  currentColumnSum = columnSums[i];
  currentTotalSum = totalSum;
  for(int j = 0; j < numberOfRows; j++)
  {
   currentTableValue = tableHolder[MATRIX_INDEX(j,i,numberOfColumns)];
   currentRowSum = rowSums[j];
   decision = MakeDecision(currentColumnSum, currentRowSum, currentTotalSum, &qHolder);
   mpz_mul_ui(currentQ, currentQ, (unsigned long int) qHolder);
   tableHolder[MATRIX_INDEX(j,i,numberOfColumns)] = decision; 
   currentColumnSum -= decision;
   currentTotalSum  -= currentRowSum;
   rowSums[j] -= decision;
   totalSum -= decision;
  }
 }
 mpz_add(qT, qT, currentQ);
 mpz_clear(currentQ);
}

int TestDiaconis()
{
 srand(55455);
 int montecarloTrials = 10000;
 int numberOfRows = 5;
 int numberOfColumns = 3;
 mpz_t qT; mpz_init(qT);mpz_set_ui(qT, 0);
 int *tableHolder = calloc(numberOfRows *numberOfColumns, sizeof(int));

 tableHolder[MATRIX_INDEX(0,0,numberOfColumns)] = 5; tableHolder[MATRIX_INDEX(0,1,numberOfColumns)] = 2;tableHolder[MATRIX_INDEX(0,2,numberOfColumns)] = 3;
 tableHolder[MATRIX_INDEX(1,0,numberOfColumns)] = 50;tableHolder[MATRIX_INDEX(1,1,numberOfColumns)] = 7;tableHolder[MATRIX_INDEX(1,2,numberOfColumns)] = 5;
 tableHolder[MATRIX_INDEX(2,0,numberOfColumns)] = 3; tableHolder[MATRIX_INDEX(2,1,numberOfColumns)] = 6;tableHolder[MATRIX_INDEX(2,2,numberOfColumns)] = 4;
 tableHolder[MATRIX_INDEX(3,0,numberOfColumns)] = 5; tableHolder[MATRIX_INDEX(3,1,numberOfColumns)] = 3;tableHolder[MATRIX_INDEX(3,2,numberOfColumns)] = 3;
 tableHolder[MATRIX_INDEX(4,0,numberOfColumns)] = 2; tableHolder[MATRIX_INDEX(4,1,numberOfColumns)] = 7;tableHolder[MATRIX_INDEX(4,2,numberOfColumns)] = 30;
 int *rowSums = GenerateRowSums(numberOfRows, numberOfColumns, tableHolder);
 int *columnSums = GenerateColumnSums(numberOfRows, numberOfColumns, tableHolder);

 int *rowSumsCopy = GenerateRowSums(numberOfRows, numberOfColumns, tableHolder);
 int *columnSumsCopy = GenerateColumnSums(numberOfRows, numberOfColumns, tableHolder);

 int totalSum = 0; for(int i = 0; i < numberOfRows; i++){totalSum += rowSums[i];} 

 for(int i = 0; i < montecarloTrials; i++)
 {
  memcpy(rowSumsCopy, rowSums, numberOfRows * sizeof(int));
  memcpy(columnSumsCopy, columnSums, numberOfColumns * sizeof(int));
  GenerateTable(numberOfRows, numberOfColumns, totalSum, rowSumsCopy, columnSumsCopy, tableHolder, qT);
 }
 mpz_div_ui(qT, qT, (unsigned long int) montecarloTrials);
 gmp_printf("Approximate Count : %3Zd : %d\n", qT , 238243776);
 mpz_clear(qT);free(rowSums);free(columnSums);free(rowSumsCopy);free(columnSumsCopy);
 free(tableHolder);
 return 0;
}

int main()
{
 TestDiaconis();
 return 0;
}