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

#define full_factorial 0	//set to 1 to also display full factorial along with truncated digits

//Calculates factorials up to 500! using the GNU Multiple Precision Arithmetic Library (GMP)
//and also the Stirling and Ramanujan approximations for the factorial function
//D.M. Gualtieri, November 3, 2020 (gualtieri@ieee.org)

//compile with gcc -o stirling_ramanujan stirling_ramanujan.c -lgmp -lm
//execute in terminal with ./stirling_ramanujan

void fact (mpz_t fact_n, int n)
{
  int i;
  mpz_t mpzi;

  mpz_init (mpzi);

  mpz_set_ui (fact_n, 1);
  if (n <= 1)
    return;
  for (i = 2; i <= n; i++)
    {
      mpz_set_ui (mpzi, i);
      mpz_mul (fact_n, fact_n, mpzi);
    }
}

void calc_stirling (mpf_t stirling, int n)
{
//Sirling approximation to factorial of n is sqrt(2*pi*n)*((n/e)^n)

  mpf_t two;
  mpf_t f;
  mpf_t e;
  mpf_init (two);
  mpf_init (f);
  mpf_init_set_str (e, "2.71828182845904523536028747135266250", 10);
  mpf_set_str (stirling, "3.141592653589793238462643383279502884", 10);
  mpf_set_str (two, "2.0", 10);
  mpf_mul (stirling, stirling, two);
  mpf_mul_ui (stirling, stirling, n);
  mpf_sqrt (stirling, stirling);
  mpf_set_ui (f, n);
  mpf_div (f, f, e);
  mpf_pow_ui (f, f, n);
  mpf_mul (stirling, stirling, f);
}

void calc_ramanujan (mpf_t ramanujan, int n)
{
//Ramanujan approximation to factorial of n is sqrt(pi)*((n/e)^n)*((8n^3)+(4n^2)+n+(1/30))^(1/6)

  double scratch;
  double e = 2.718281828459;
  mpf_t f;
  mpf_init (f);
  mpf_set_str (ramanujan, "3.141592653589793238462643383279502884", 10);
  mpf_sqrt (ramanujan, ramanujan);
  scratch = (double)n/e;
  mpf_set_d (f, scratch);
  mpf_pow_ui (f, f, n);
  mpf_mul (ramanujan, ramanujan, f);
  scratch = pow((8*pow(((double) (n)),3)+(4*pow(((double) (n)),2))+(double)n+(1/30)),0.166666667);
  mpf_set_d (f, scratch);
  mpf_mul (ramanujan, ramanujan, f);
}

int main (int argc, char *argv[])
{
  int i;
  int factorial;
  int stir;
  int ramanujan_fact;
  char result[2500] = "";
  mpz_t n;
  mpf_t stirling;
  mpf_t ramanujan;
  FILE *outdata;

  if ((outdata = fopen ("stirling_ramanujan.dat", "w")) == NULL)
    {
      printf ("\nOutput file cannot be opened.\n");
      exit (1);
    }

  //datafile header
  if (full_factorial == 1)
    {
      fprintf (outdata,
	       "n\tn!\tn!(truncated)\tstirling n!\tramanujan n!\tstirling error(%%)\tramanujan error(%%)\n");
    }
  else
    {
      fprintf (outdata,
	       "n\tn!\tstirling n!\tramanujan n!\tstirling error(%%)\tramanujan error(%%)\n");
    }

  //initialize variables
  mpz_init (n);
  mpf_init (stirling);
  mpf_init (ramanujan);

  mpf_set_str (stirling, "3.141592653589793238462643383279502884", 10);
  mpf_set_str (ramanujan, "3.141592653589793238462643383279502884", 10);

  for (i = 0; i <= 500; i++)
    {
      printf ("%d\t", i);
      fprintf (outdata, "%d\t", i);
      fact (n, i);
      mpz_get_str (result, 10, n);
      if (full_factorial == 1)
	{
	  printf ("%s\t", result);
	  fprintf (outdata, "%s\t", result);
	}
      result[9] = '\0';
      factorial = atoi (result);
      printf ("%s\t", result);
      fprintf (outdata, "%s\t", result);
      calc_stirling (stirling, i);
      gmp_snprintf (result, 10, "%.Ff", stirling);
      stir = atoi (result);
      printf ("%s\t", result);
      fprintf (outdata, "%s\t", result);
      calc_ramanujan (ramanujan, i);
      gmp_snprintf (result, 10, "%.Ff", ramanujan);
      ramanujan_fact = atoi (result);
      printf ("%s\t", result);
      fprintf (outdata, "%s\t", result);
      //calculate differences %
      printf ("%f\t", 100 * (double) (factorial - stir) / factorial);
      fprintf (outdata, "%f\t",
	       100 * (double) (factorial - stir) / factorial);
      printf ("%f\n", 100 * (double) (factorial - ramanujan_fact) / factorial);
      fprintf (outdata, "%f\n",
	       100 * (double) (factorial - ramanujan_fact) / factorial);
    }
  printf ("done\n");
  fclose (outdata);
  return 0;
}