/*
//  frac.cpp
//
//  match the expected value with a Bernoulli trial
//  (multinomial dist - really binomial)
//
//
//  theory
// the fraction of single hits =  (M|1) / sum_1_to_n (M|i)
//
// a few tricks for factorials and whatever
//
*/
// old style includes because this g++ is a little dated
//#include <exception>
#include <iostream.h>
#include <math.h>

// return the log of the binomial coefficient
double binomial( int m , int n)
{
	int o;
	int i;
	
	o = m-n;
// this should be correct, but my screwy g++ doesn't have the include
//	if( o < 0 ) throw range_error ;
	if( o < 0) return 0.;
	if( o == 0) return 0.;
	if( o < n) { i = n; n = o ; o = i;}

	double accum, x, lnx;
//  count down for the M!/(M-N)!
	accum = 0;
	for( i=m; i> n; i--)
	{
		x  =(double)i;
		lnx  = log(x);// this is clearer
		accum += lnx;
	}
	for( i=1; i<= o; i++)
	{
		x = (double)i;
		lnx = log(x);
		accum -= lnx;
	}
	return accum;
}




int main()
{
	double binomial( int, int);

	
	cout << "please enter the total number of contacts \n";
	int total ; cin >> total;
	cout << "please enter the total number of call signs\n";
	int calls ;
	cin >> calls;	
	cout << "please enter the total number of unique call signs\n";
	int unique;
	cin >> unique;
	double fract = (double) unique/ (double)calls;
	cout << "There are "<<fract<<" singletons\n";


	double bico[total];
	int i,j;
	
	for( i=0; i< total;i++)
		bico[i] = binomial(total,i+1);

	for( i=calls; i< calls* 500; i+=calls)
	{
		double accum, aden;
		accum = 0.;
		aden = 0.;
		accum = exp(bico[0] -log((float)i) + log( 1.-1./(float)i)*(total-1));
		
		for( j=0; j< total; j++)
		{
		aden += exp(bico[j] -log((float)i)*(j+1) + log( 1.-1./(float)i)*(total-j-1));
		}
		cout << i <<" accum " << accum << " aden "<< aden << " fraction "<< accum/aden <<"\n";
		if(accum/aden > fract) break;
	}//i
	int k;
	k = i;
	i = i-calls;	
	for( ; i< k; i++)
	{
		double accum, aden;
		accum = 0.;
		aden = 0.;
		accum = exp(bico[0] -log((float)i) + log( 1.-1./(float)i)*(total-1));
		
		for( j=0; j< total; j++)
		{
		aden += exp(bico[j] -log((float)i)*(j+1) + log( 1.-1./(float)i)*(total-j-1));
		}
		cout << i <<" accum " << accum << " aden "<< aden << " fraction "<< accum/aden <<"\n";
		if(accum/aden > fract) break;
	}//i
	
}