/**********************************************************************/
/*                                                                    */
/* Program: fac.c                                                     */
/*                                                                    */
/* This program calculates large n ! (factorials) exactly.            */
/*                                                                    */
/*                                                                    */
/* Author: Stefan Spaennare, Lund, Sweden                             */
/*                                                                    */
/* E-mail: stefan@spaennare.se                                        */
/*                                                                    */
/* First version:         1993                                        */
/* Latest update:         2004-06-23                                  */
/*                                                                    */
/*                                                                    */
/* Benchmark                                                          */
/* =========                                                          */
/*                                                                    */
/* The test was run on a Intel Celeron (Tualatin core) computer       */
/* with 256 kbyte cache memory and 512 Mbyte primary memory.          */
/* The operating system was Red Hat 7.2 Linux. Compiler tested        */
/* was GNU GCC 2.96 (gcc) for Linux. Under Windows XP the compiler    */
/* GNU DJGPP 3.0.4 (gcc, djgpp) for DOS was used.                     */
/*                                                                    */
/* Numbers tested:                                                    */
/* ---------------                                                    */
/*                                                                    */
/* na =    10000 =  10^4 = 10*10^3                                    */
/*                                                                    */
/* Results:                                                           */
/* --------                                                           */
/*                                                                    */
/* n:  Computer:         CPU-time (s):  Compiler options:             */
/* ================================================================   */
/*                                                                    */
/* na  Celeron 1200 MHz       2.90      gcc -O3                       */
/* na  Celeron 1200 MHz       2.69      gcc -O2 (djgpp)               */
/*                                                                    */
/*                                                                    */
/* Old benchmark                                                      */
/* =============                                                      */
/*                                                                    */
/* n=10000                                                            */
/*                                                                    */
/* Computer:                CPU-time (s)           Compilation:       */
/*                                                                    */
/* Pentium 180 MHz              8.54                  gcc -O3         */
/*                                                                    */
/**********************************************************************/

/********************************************************************************/
/*                                                                              */
/* Notice                                                                       */
/* ======                                                                       */
/*                                                                              */
/* I make no warranties that this program is (1) free of errors, (2) consistent */
/* with any standard merchantability, or (3) meeting the requirements of a      */
/* particular application. This software shall not, partly or as a whole,       */
/* participate in a process, whose outcome can result in injury to a person or  */
/* loss of property. It is solely designed for analytical work. Permission to   */
/* use, copy and distribute is hereby granted without fee, providing that the   */
/* header above including this notice appears in all copies.                    */
/*                                                                              */
/*                                                            Stefan Spaennare  */
/*                                                                              */
/********************************************************************************/


#include <math.h>
#include <stdio.h>
#include <unistd.h>
#include <string.h>
#include <stdlib.h>
#include <time.h>


double timetic=(double)(CLOCKS_PER_SEC);

char fname[80]="fac.dat";  /* Output filename */

double pi=3.14159265358979323;


void fppp(t,outfile)

int t; 
FILE *outfile;

{
  if (t > 9999999) {
     fprintf(outfile,"%8d",t);
  }
  else if (t > 999999) {
     fprintf(outfile,"0%7d",t);
  } 
  else if (t > 99999) {
     fprintf(outfile,"00%6d",t);
  }
  else if (t > 9999) {
     fprintf(outfile,"000%5d",t);
  } 
  else if (t > 999) {
     fprintf(outfile,"0000%4d",t);
  }
  else if (t > 99) {
     fprintf(outfile,"00000%3d",t);
  } 
  else if (t > 9) {
     fprintf(outfile,"000000%2d",t);
  }
  else {
     fprintf(outfile,"0000000%1d",t);
  } /* if */

} /* fppp */



 
int main(argc,argv) 
int argc;
char *argv[];

{

   FILE *outfile;

   int i,k,n;
   int iel,m;

   double el,r,q2;
   double c2,cl,jd;
   double bas,basi;

   clock_t time3,timediff3;

   int *a; 

   if (argc != 2) {
      printf("Usage: %s ni\n",argv[0]);
      exit(0);
   } /* if */

   n=atoi(argv[1]);


   bas=100000000.0;
   basi=0.00000001;

   c2=1/(8.0*log(10.0));
   cl=1/log(10.0);

   el=0.0;
   
   for (i=1; i<=n; i++) {
      el=el+log((double)(i));
   } /* for i */

   el=el*cl;

   printf("\n");

   printf("%d ! = %f %d\n\n",n,exp((el-(double)((int)(el)))/cl),(int)(el));

   printf("Calculating ...\n");
 
   a = (int *)calloc((int)(el)/8+10,sizeof(int));

   time3=clock(); 

   a[1]=1;
   el=1.0;

   for (i=1; i<=n; i++) {
      jd=(double)(i);
      el=el+c2*log(jd);
      iel=(int)(el+0.5);
      r=0.0;
      for (k=1; k<=iel; k++) {
         q2=(double)(a[k])*jd+r+0.5;
         r=(double)((int)(basi*q2)); 
         a[k]=(int)(q2-bas*r);
      } /* for k */
    } /* for i */

   timediff3=clock()-time3; 
   printf("\nTime %8d !:      %10.2f\n",n,(double)(timediff3)/timetic);

   printf("\n");

   printf("Calculation finished, output to %s.\n\n",fname);

   outfile=fopen(fname,"w");


   m=0;

   iel=(int)(el+0.5);

   if (a[iel] > 0) {
      fprintf(outfile,"%8d",a[iel]);
      iel--;
      fppp(a[iel],outfile);
      m=m+2;
   } 
   else {
      iel--;
      fprintf(outfile,"%8d",a[iel]);
      m++;
   } /* if */

   iel--;

   for (i=iel; i>=1; i--) {
      fppp(a[i],outfile);
      m++;
      if ((m % 9)==0) {
         fprintf(outfile,"\n");
      } /* if */
   } /* for i */

   fprintf(outfile,"\n");

   printf("\n");


} /* end fac */