/*********************************************************************
   Returns the binomial coefficient (n k),  0 <= k <= n, as a floating
   point number.
   C.A. Bertulani        May/15/2000
*********************************************************************/

typedef double Number;

#include <math.h>
#include<iostream.h>

int main(){
	Number bico(int , int );
	int j, n, k;
	j=1;
	cout << "\n Enter n and k. To stop, enter k<0\n";
	cin >> n >> k;
    while(k>=0){
		if(j>10){ cout << "\n My patience is over. Stop, please!\n";
			break;
		}
		if(j!=1){
			cout << "\n Enter n and k. To stop, enter k<0.\n";
			cin >> n >> k;
		}
		cout << "\n Here is (n   k): " << bico(n,k);
		j=j+1;
	}
	return 0;
}
/*********************************************************************
   Returns the binomial coefficient (n k),  0 <= k <= n, as a floating
   point number.
*********************************************************************/
Number bico(int n, int k)
{
	Number factln(int n);

	return floor(0.5+exp(factln(n)-factln(k)-factln(n-k)));
	/* the floor function clean up roundoff error for smaller values of n and k. */
}
/********************************************************************
   Returns the value of ln[Gamma(xx)] for xx > 0
********************************************************************/

Number gamma_ln(Number xx)
{
	Number x,y,tmp,ser;
	static Number cof[6]={76.18009172947146,-86.50532032941677,
		24.01409824083091,-1.231739572450155,
		0.1208650973866179e-2,-0.5395239384953e-5};
	int j;

	y=x=xx;
	tmp=x+5.5;
	tmp -= (x+0.5)*log(tmp);
	ser=1.000000000190015;
	for (j=0;j<=5;j++) ser += cof[j]/++y;
	return -tmp+log(2.5066282746310005*ser/x);
}
/********************************************************************
   Returns ln(n!)
********************************************************************/
Number factln(int n)
{
	Number gamma_ln(Number xx);
	static Number a[101];  /* a static array is automatically initialized to zero */

	if (n < 0) cerr << "Negative factorial in routine factl\n";
	if (n <= 1) return 0.0;
	if (n <= 100) return a[n] ? a[n] : (a[n]=gamma_ln(n+1.0));   /* In range of table */
	else return gamma_ln(n+1.0);  /* Out of range of table */
}
/********************************************************************/