#include "misc.h"
#include "sets.h"
#include "combin.h"

/*
** _allPermutations: recursive helper function of CombinAllPermutations,
** to list all permutations of the integers 0..n-1 inclusive.
** preArray is the array of
** things already filled up, so our job is to fill the rest of
** it with all possible permutations.  The top call is
**
**      _allPermutations(n, 0, [array of n ints], fcn)
** 
** where fcn is a function to call for each permutation.  The function
** will take the array and it's size as arguments.  (*fcn) should
** return 0 to continue listing the permutations, or non-zero to exit.
** The non-zero return value will be returned to the calling function.
*/

static Boolean _allPermutations
    (int n, int i, int *preArray, Boolean (*fcn)(int, int *))
{
    int j;
    if( n == i) /* output! */
	return fcn(n, preArray);

    for(j=0; j < n; j++)    /* put in slot i all j's not already appearing */
    {
	int k;
	for(k=0; k<i; k++)  /* see if this j's already been used */
	{
	    if(preArray[k] == j)
		break;
	}
	if(k == i)  /* this j hasn't appeared yet */
	{
	    int result;
	    preArray[i] = j;
	    if((result = _allPermutations(n, i+1, preArray, fcn)))
		return result;
	}
    }
    return false;
}


/*
** The externally-viewable function for calling all combinations.
*/
Boolean CombinAllPermutations(int n, Boolean (*fcn)(int, int *))
{
    int array[n];
    return _allPermutations(n, 0, array, fcn);
}



/*
** Helper function for CombinAllCombinations, similar to above.
*/
static Boolean _allCombinations
    (int n, int m, int i, int *array, Boolean (*fcn)(int, int *))
{
    int j;
    if(m == i)  /* output! */
	return fcn(m, array);

    for(j=array[i-1]+1; j < n; j++)
    {
	int result;
	array[i] = j;
	if((result = _allCombinations(n, m, i+1, array, fcn)))
	    return result;
    }
    return false;
}

/*
** Externally viewable function.
*/
Boolean CombinAllCombinations(int n, int m, Boolean (*fcn)(int, int *))
{
    int array[m], i;
    assert(n >= 0 && m >= 0);
    assert(m <= n);
    if(m == 0)
	return 0;
    for(i=0; i<n; i++)
    {
	int result;
	array[0] = i;
	if((result = _allCombinations(n, m, 1, array, fcn)))
	    return result;
    }
    return false;
}


/*
** The following routines are non-recursive.  You call them one at
** a time when you are ready to take the next combination, rather
** than the above functions which call *you* once for each combination.
*/
COMBIN *CombinZeroth(int N, int M, unsigned *A)
{
    COMBIN *c = (COMBIN*)Calloc(1,sizeof(COMBIN));
    int i;
    assert(M>=0 && N>=0);
    assert(M <= N);
    c->n = N;
    c->m = M;
    c->array = A;
    for(i=0; i < M; i++)
	A[i] = i;
    return c;
}

COMBIN *CombinIth(int N, int M, unsigned *A, unsigned long long I)
{
    COMBIN *c = CombinZeroth(N,M,A);
    unsigned long long currentIndex = 0;
    int element;
    
    assert(I < CombinChoose(N,M));

    /*
    ** We do a kind of weird binary search to find the Ith index.  We keep
    ** incrementing the first array element until we've gone too far, then do
    ** the same with the next element, and the next, and so on.  To increment
    ** the e'th element by 1 (everything starting at 0) requires skipping
    ** past Choose(n-1-e,m-1-e) combinations.
    */
    element = 0;
    while(element < M)
    {
	unsigned long long IofNextE =
	    currentIndex + CombinChoose(N-(A[element]+1), M-(element+1));
	if(IofNextE > I)
	    ++element;
	else
	{
	    int i;
	    for(i=element; i < M; i++)
		++A[i];
	    currentIndex = IofNextE;
	}
    }
    return c;
}


/*
** This should be done more cleverly
*/
Boolean CombinSkipN(COMBIN *c, int n)
{
    int i;
    for(i=0; i<n && CombinNext(c); i++)
	;
    /* if i==n, we didn't run out of combinations */
    return (i == n);
}

Boolean CombinNext(COMBIN *c)
{
    int i;
    
    if(c->m == 0)   /* special case */
	return false;

    i = c->m-1;
    do
    {
	if(++c->array[i] <= c->n - c->m + i)
	    break;
    } while(i--);

    if(i == -1)
	return false;
    else
    {
	int j;
	for(j=i+1; j < c->m; j++)
	    c->array[j] = c->array[j-1]+1;
	return true;
    }
}


Boolean CombinAssign(COMBIN *c, unsigned newComb[c->m])
{
    int i;
    for(i=0; i < c->m; i++)
    {
	if(i > 0) assert(newComb[i-1] < newComb[i]);
	assert(newComb[i] < c->n);
	c->array[i] = newComb[i];
    }
    return true;
}


/*
** The *best* way to do this would be to precompute the factorizations of
** all the numbers I'd ever use.  Then, the numers and denoms could be
** cancelled using addition, leaving only multiplications (no divisions).
** ie, overflow means overflow, and all representable answers would be
** computed properly.
**
** Note: all "impossible" values of n and m return 0.
*/
unsigned long long CombinChoose(int n, int m)
{
    unsigned long long result = 1;
    int low, high, i;
    SET *denoms;

    if(n < 0 || m < 0) return 0;
    if(n == 0) return m == 0 ? 1 : 0;
    if(m > n) return 0;

    if(m < n-m) { low = m; high = n-m; }
    else {low = n-m; high = m; }

    denoms = SetAlloc(low+1);

    for(i=2; i <= low; i++)
	SetAdd(denoms, i);

    for(i=n; i > high; i--)
    {
	int j;
	long long old = result;
	result *= i;
	if(result / i != old)	/* overflow */
	{
	    Warning("CombinChoose(%d,%d): overflow", n,m);
	    SetFree(denoms);
	    return 0;
	}
	for(j=2; j <= low; j++)
	    if(SetIn(denoms,j) && result % j == 0)
	    {
		result /= j;
		SetDelete(denoms, j);
	    }
    }

    assert(SetCardinality(denoms) == 0);
    SetFree(denoms);
    return result;
}


double CombinCumulativeBinomialCumulative(int n, int k, double x)
{
    if(x == 0.0 || x == 1.0)
	return 0.0;
{
    double sum = 0, x_j = IntPow(x,k), o_x = IntPow(1-x, n-k);
    int j;

    assert(k >= 0 && n >= k);

    for(j=k; j <=n; j++)
    {
	sum += CombinChoose(n,j) * x_j * o_x;
	x_j *= x;
	o_x /= 1-x;
    }
    return sum;
}
}


double CombinCumulativeBinomialDensity(int n, int k, double x)
{
    assert(false);
    return 0.;
}
