/* Copyright (C) 1994 by Thomas Glen Smith.	All Rights Reserved. */
/* binomp APL2 V1.0.0 **************************************************
* binomial.												 *
***********************************************************************/
#define INCLUDES 0
#include "includes.h"
void binomp(m,n,ret)
double *m,*n,*ret;
{
	Binomp; Factorl; Floor; Mod;
	extern int aplerr;
	double a,b,c,den,k,num,p,q,r,s;

	a = *m; /* Left. */
	b = *n; /* Right. */
	for(;;) { /* So I can use break. */
		if (floor(a) == a && floor(b) == b) { /* a, b are integers */
			if (a == b || a == 0e0) {
				*ret = 1e0;
				break;
			}
			if (b < 0e0 && a < b) { /* return (-1*b_a)X1!(|a+1) */
				p = -(b+1); /* New left = _right+1. */
				q = -(a+1); /* New right = -left+1. */
				binomp(&p,&q,&r);
				s = mod(b - a, 2e0) == 0e0 ? 1e0 : -1e0;
				*ret = s * r;
				break;
			}
			if (b >= 0e0)
				if (a > b || a < 0e0) {
					*ret = 0e0;
					break;
				}
				else;
			else if (a < 0e0 && a > b) {
				*ret = 0e0;
				break;
			}
			k = b - a;
			num = b;
			while((b -= 1e0) > k) num *= b;
			den = a;
			while((a -= 1e0) > 1e0) den *= a;
			*ret = num / den;
		} else { /* either a or b isn't an integer */
			if (b < 0e0 && floor(b) == b /* b is a negative integer */
				&& floor(a) != a /* *m is not an integer */) {
				aplerr = 40;
				*ret = 0e0;
				break;
			}
			*ret = (factorl(b) / (factorl(a) * factorl(b - a)));
		}
		break;
	} /* end for(;;) */
}