#include <iostream.h>
#include <stdlib.h>
#include <conio.h>
#include <complex>

// nCr*r! subroutine
long double combination(int num, int rum)
{
int temp;
long double fact=1;
	temp = num-rum;
    while (num > temp)
        fact *= num--;
        return(fact);
}

// factorial subroutine
long double factorial(int num)
{
    if (num == 1)
        return(1);
    else
        return(num * factorial(num - 1));
}

void main()
{
long double resn;
long double muts;
long double libsiz;
long double combins;
long double facto;
long double maxlib;
float coverage;
cout<<"copyright 2003 by Paul A Dalby, University College London"<<endl;
cout<<""<<endl;
  cout<<"enter total number of residues in target gene"<<endl;
  cin>>resn;
  cout<<"enter maximum number of mutations per gene"<<endl;
  cin>>muts;
    //calculate n fac/n-r fac
    facto = (combination(resn, muts));
	    //calculate combinations
        combins = facto/factorial(muts);
          //input choice
	      float choice;
          cout<<"enter allowed number of possible variants at each position"<<endl;
          cout<<"Eg. 1 for DNA transitions or transversions only (1 base mutation/gene)"<<endl;
          cout<<"Eg. 2 for DNA transitions or transversions only (> 1 base mutation/gene)"<<endl;
          cout<<"Eg. 3 for only other DNA bases (1 base mutation/gene)"<<endl;
	      cout<<"Eg. 4 for all DNA bases (> 1 base mutation/gene)"<<endl;
	      cout<<"Eg. 5.7 for only other amino acids in EPPCR (1 residue mutation/gene)"<<endl;
	      cout<<"Eg. 6.7 for all amino acids in EPPCR (>1 residue mutation/gene)"<<endl;
          cin>>choice;
	// calculate library size
	maxlib = combins*pow(choice,muts);
	cout<<"maximum library size is "<<maxlib<<endl;
	cout<<"(based on "<<choice<<"average number of variants possible per site)"<<endl;
       // calculate library size required for any given % coverage
	   cout<<"enter the % coverage of maximum library required (1-100%)"<<endl;
       cin>>coverage;
	     libsiz = -1*maxlib*log(1-(coverage/100));
         cout<<"required number of clones for "<<coverage<<"% coverage is "<<libsiz<<endl;
	cout<<"Note - this uses a Poisson distribution and assumes an UNBIASED library"<<endl;
cout<<"press any key to continue"<<endl;
getch();
}