/*********************************************************************/
/*    Jeff Balsley                                                   */
/*    gcc -lm probability2.c -o probability2                         */
/*                                                                   */
/*    This is basically the binomial distribution as taken from      */
/*    "Data reduction and Error analysis for the physical sciences"  */
/*    by Philip R. Bevington and D. Keith Robinson                   */
/*                                                                   */
/*********************************************************************/

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

double prob(int x, int num_events, double p);
double fac(int n);
double choose(double n, double k);

int main(void)
{
    double p;           /* The probability of the event */
    int num_events;  /* The total number of events */
    int x;           /* The number of occurances we want to see */

    printf ("What is the total number of times the event will be tried >");
    scanf ("%d", &num_events);
    printf ("What is the probability of the event >");
    scanf ("%lf", &p);
    printf ("What is the number of occurances you want to see >");
    scanf ("%d", &x);

    printf ("The probability of the event with prob %f occuring %d \n",
        p, x);
    printf ("times in %d tries is %f.\n", num_events,
        prob(x, num_events, p));

    return(0);
}

double prob(int x, int num_events, double p)
{
    double probability;

    probability = ( choose(num_events, x) ) * pow(p, x) *
        pow(1-p, num_events-x);

    /* debugging */

    return(probability);
}

double choose(double n, double k)
{
    int i;
    double product = 1;

    /* check a few things first */
    if (k > n ){
        return(0);
    }else if( k == n || k == 0 ){
        return(1.0);
    }

    if (k < n/2.0){
        for (i=n; i > (n-k); --i){
            product = product * i;
        }
        /*printf("1 product is %f\n", product);*/
        return( product/fac(k) );
    }
    if ( k >= n/2.0 ){
        for (i=n; i > k; --i){
            product = product * i;
        }
        /*printf("2 product is %f\n", product);*/
        return( product / fac( n-k ) );
    }
}

double fac(int n)
{
    int i;
    double product = 1;

    if (n > 1){
        for(i=n; i>1; --i){
            product = product * i;
        }
    }else if( n < 0){
        printf("Can't take factorial of a negative number!!\n");
        return(0);
    }
    /*printf("%d factorial is %f\n", n, product); */
    return(product);
}