// Dackral Phillips
// COMP-6700
// Assignment 5A
// 7/23/2K1

// In this assignment, we are given a value for t, the number of degrees of freedom, and the number
// of tails of a given t distribution and are supposed to produce a probability from the information.
// This probability is calculated from an integration using Simpson's rule to calculate it.

// Without further delay, let's dive in and code:  Necessary Libraries

#include <stdlib.h>
#include <iostream.h>
#include <strings.h>
#include <iomanip.h>
#include <math.h>
#define PI 3.14159

   double integrate(double lower, double upper, int n, double t);
   double t-gamma(int n);
   double gamma(double n);
   double t-function(double t, int n);
   double simpsons(double lower, double upper, int n, double t);

   int main()
   {
      double t;
      int n;
      int tails;
      double result;
      char ans = 'y';
      int stop = 0;
   int test = 1;
   
      while (!(stop))
      {
         system("clear");
         cout << endl << endl 
            << "Please enter a probability (t): ";
         cin >> t;
         cout << endl << endl
            << "Please enter the degrees of freedom (n): ";
         cin >> n;
         cout << endl << endl
            << "Please enter the number of tails (1 or 2): ";
         cin >> tails;
      
         result = integrate(0, t, n);
      
         if (tails == 1)
         {
            result = result + .5;
         }
         else
         {
            result = result * 2;
         }
      
      cout << endl << endl 
           << "Test\t\t t \t\t\tDegrees of Freedom\tTails\t\tProbability" << endl
           << "____\t\t___\t\t\t__________________\t_____\t\t___________" << endl
           <<  test << "\t\t" << t << "\t\t\t" << n << "\t\t" << tails << "\t\t" << result; 
      
         cout << endl << endl 
            << "Would you like to calculate another probability? [Y,y,N,n]: ";
         cin >> ans;
      
         if ((ans == 'Y') || (ans == 'y'))
         {
            stop = 0;
         test = test + 1;
         }
         else
         {
            stop = 1;
         }
      }
      return 0;
   }

   double gamma(double n)
   {
      double answer;
   
      if (n == 1)
      {
         answer = 1;
      }
      else if (n == .5)
      {
         answer = sqrt(PI);
      }
      else
      {
         answer = (n - 1) * gamma(n - 1);
      }
   
      return answer;
   }

   double t-gamma(int n)
   {
      double answer;
   
      answer = (gamma(((n + 1)/2)) / (sqrt(n * PI) * gamma(n / 2)));
      return answer;
   }

   double t-function(double t, int n)
   {
      double answer;
   
      answer = pow((1 + ((t * t) / n)), (((n + 1) / 2) * -1));
      return answer;
   }

   double simpsons(double lower, double upper, int n, double t)
   {
      double answer = 1;
      double oldAnswer = 0;
      double e = 0.01;
      double w = (upper - lower) / n;
      int flag = 0;
      int i;
   
   
      while (abs(answer - oldAnswer) <= e)
      {
         oldAnswer = answer; 
         answer = t-function(lower, n);
      
         for (i = (lower + w); i < upper; i = (i + w))
         {
            if (flag == 0)
            {
               answer = answer + (4 * t-function(i, n));
               flag = 1;
            }
            else if (flag == 1)
            {
               answer = answer + (2 * t-function(i, n));
               flag = 0;
            }
            times = times + 1;
         }
         answer = answer + t-function(upper, n);
         answer = answer * (w / 3);
         w = w * 2;
      }
   
      return answer;
   }

   double integrate(double lower, double upper, int n, double t)
   {
      double answer;
   
      answer = t-gamma(n) * simpsons(lower, upper, n, t);
      return answer;
   }