Source Code

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

#define MAX_POP_SIZE 10000

int popB[MAX_POP_SIZE], popR[MAX_POP_SIZE];

void GetPara(int *popSizePtr, double *freqPtr);
void InitPop(int *popPtr, int popSize, double initFreq);
void CreateNextGenPanmictic(int *curPop, int *nextPop, int popSize);
int RandInt(int max);
int MyRound(double x);

int main(void) {
  int gen, *currentPtr, *nextPtr;
  int popSize, numGen = 100000;
  double initFreq;
  
  GetPara(&popSize, &initFreq);  /* get the values from the user */
  
  if (popSize > MAX_POP_SIZE) {  /* Make sure we have big enough arrays */
    fprintf(stderr, "pop size: %d is bigger than the max: %d\n",
	    popSize, MAX_POP_SIZE);
    exit(EXIT_FAILURE);
  }
  srand(time(NULL));
  InitPop(popB, popSize, initFreq);  /* randomly assign alleles */

  for (gen = 0; gen < numGen; gen++) {
    if (gen % 2 == 0) {
      currentPtr = popB;
      nextPtr = popR;
    } else {
      currentPtr = popR;
      nextPtr = popB;
    }
    
    CreateNextGenPanmictic(currentPtr, nextPtr, popSize);
  }
  return 0;
}

/* Initialize the parameters */
void GetPara(int *popSizePtr, double *freqPtr) {
  fprintf(stderr, "Population Size?: ");
  scanf("%d", popSizePtr);
  fprintf(stderr, "Init freq. of allele 1?: ");
  scanf("%lf", freqPtr);
  return;
}

/*
 * Assume that there are two alleles: 0 and 1
 * Initialize the population, so that the allele frequency of allele 1
 * is initFreq.
 */
void InitPop(int *popPtr, int popSize, double initFreq) {
  int i, numAlleles, cntr, index;
  numAlleles =  MyRound(popSize *  initFreq);
  
  for (i = 0; i < popSize; i++) {
    popPtr[i] = 0;
  }
  
  cntr = 0;
  while (cntr < numAlleles) {
    index = RandInt(popSize - 1);
    if (popPtr[index] == 0) {
      popPtr[index] = 1;
      cntr++;
    }
  }
  return;
}

void CreateNextGenPanmictic(int *curPop, int *nextPop, int popSize) {
  int i, parent;
  
  for (i = 0; i < popSize; i ++) {
    parent = RandInt(popSize - 1);
    nextPop[i] = curPop[parent];
  }
  
  return;
}

/* 
 * Returns a uniform random integer between 0 and max (ends inclusive) 
 * It uses rand(), so you need to srand() in main()
 */
int RandInt(int max) {
  int r = rand();
  return r % (max+1);
}

/* 
 * Round the floating points to the nearest integer 
 * For example, 9.5 <= x < 10.5 become 10
 * This function can cause integer over-flow when x is large
 */
int MyRound(double x) {
  int intPart;
  
  if (x > 0) {
    return (int) (x + 0.5);
  } else {
    return (int) (x - 0.5);
  }
}