/*

COMBIN.C:  A computer program for verifying numerically the identities
developed in the paper "Combinatorial identities associated with CFTP",
by G.O. Roberts and J.S. Rosenthal, available from

	http://probability.ca/jeff/research.html


Copyright (c) 2001 by Jeffrey S. Rosenthal.

Licensed for general copying, distribution and modification according to
the GNU General Public License (http://www.gnu.org/copyleft/gpl.html).


SAVE AS: combin.c

COMPILE WITH:  cc combin.c -lm

RUN WITH:  a.out

MODIFY VALUES BY:  Changing the line below having comment "For now".

DESIRED RESULT:  If 0 <= n <= N and 0 <= i <= N, then the two output values,
                 "Totalval" and "1/(N+1)", should be identical.

*/



#include <stdio.h>
#include <stdlib.h>
#include <math.h>

main()

{

    double Q(), Adiff(), B(), C(), R();
    double fact();
    int N,n,i;
    double totalval;
    int tmpadd;
    double tmpval;
    int j, l, m, k;

    N=40; n=10; i=23;  /* For now -- these values may be modified!! */

    totalval = 0.0;
    for (j=0; j <= N; j++) {  /* Value of top process, from N. */
	for (l=0; l <= imin(n,N); l++) { /* Num holds of top. */
	    for (m=0; m <= imin(n,N-l); m++) { /* Num holds of bottom. */
		k = j-N+imin(N,l+m);  /* Value of bottom process, from 0. */
		tmpadd = 0;
		if (k > j)
		    printf("ERROR, k >= j HERE.\n");
		else if (k < 0)
		    /* printf("NOTE: k < 0 HERE. (j=%d,l=%d,m=%d,k=%d)\n",
		    					j,l,m,k) */ ;
		else {
                    if (j==i)
		        tmpadd = tmpadd + l;
		    if (k==i)
			tmpadd = tmpadd + m;
		    if ((k<=i) && (i<=j))
		    tmpadd = tmpadd + 1;
		}
		tmpval = tmpadd * R(N,n,l,m,k) / (N+1);
		totalval = totalval + tmpval;
                if (tmpval<0.0) {
		    printf("NEG! l=%d, m=%d, k=%d, tmpval=%f, totalval=%f\n", 
		    				l,m,k,tmpval, totalval);
		}
	    }
	}
    }

    printf("N=%d, n=%d, i=%d:  Totalval = %.16f;  1/(N+1)=%.16f \n", 
					N,n,i,totalval, 1.0/(N+1));

}

int imin(int a, int b)
{
    if (a < b)
	return(a);
    else
	return(b);
}

int imax(int a, int b)
{
    if (a > b)
	return(a);
    else
	return(b);
}

/* Q(n,x) = Prob[ssrw run for time n, ends up at x] */
double Q(int n, int x)
{
    double fact();
    if ( (-n <= x) && (x <= n) && ( 2*((x+n)/2) == x+n ) )
	return( pow(0.5,n) * fact(n) / fact((n+x)/2) / fact((n-x)/2) );
    else
	return(0);
}

double fact(int n)
{
    if (n<=1)
	return(1.0);
    else
	return(n*fact(n-1));
}

double Adiff(int n, int a, int b, int j)
{
    double resultval;
    int i;

    if ( (a*b>0) || ((b-a)*(j-b)>0) ) {
	printf("ERROR, NOT SIGN-ALTERNATING: n=%d,a=%d,b=%d,j=%d.\n", n,a,b,j);
    } else if ( (a==0) && (b==0) ) {
	return(0.0);
    } else {
	resultval =  0.0;
	for (i=0; i<= n / iabs(b-a) / 2; i++) {
	    resultval = resultval 
			+ Q(n, iabs(a) + (2*i+1)*iabs(b-a) + iabs(j-b))
			- Q(n, iabs(a) + (2*i+2)*iabs(b-a) + iabs(j-a));
	}
	if (resultval < 0) {
	    printf("ERROR, resultval<0 (==%f); n=%d, a=%d, b=%d, j=%d.\n", 
				resultval,n,a,b,j);
	    return(0.0);
	} else {
	    return( resultval );
	}
    }
}

/* B(n,a,b,j) = prob[ max >= a, min <= b, end at j]. */
double B(int n, int a, int b, int j)
{
    if ( (j>imax(a,b)) || (j<imin(a,b)) )
	return( 0.0 );
    else
	return( Adiff(n,a,b,j) + Adiff(n,b,a,j) );
}

/* C(n,a,b,j) = prob[ max == a, min == b, end at j]. */
double C(int n, int a, int b, int j)
{
    if ( (a<0) || (b>0) )
	printf("ERROR, HAVE a<0 AND/OR b>0 HERE.\n");
    else if ( (a==0) && (b==0) ) {
        if ( (n==0) && (j==0) )
	    return( 1.0 );
	else
	    return( 0.0 );
    } else
	return( B(n,a,b,j) - B(n,a+1,b,j) - B(n,a,b-1,j) + B(n,a+1,b-1,j) );
}

double R(int N, int n, int l, int m, int k)
{
    if ( (l+m<=n) && (k >= 0) && (k <= N) && (l >= 0) && (m >= 0) )
	return( C(n, l, -m, k-m) );
    else
	return( 0.0 );
}

int iabs(int a)
{
    if (a < 0)
	return(-a);
    else
	return(a);
}