/*
  Copyright (C) 2019 Alessandro Languasco
*/

#include <time.h>
#include <fftw3.h> 
#include <stdio.h>
#include <stdlib.h>
#include <math.h> 
 
#include <errno.h> 

#define REAL 0
#define IMAG 1

/*   for checking results */
void prntc(fftwl_complex* result, unsigned long long int  NUM_POINTS) { 
  int i;
  for (i = 0; i < NUM_POINTS; ++i) {
  printf("%Lf\n", result[i][REAL]);
  printf("%Lf\n", result[i][IMAG]);
  }
}
 

 
int main(int argc, char *argv[]) {
  ptrdiff_t  a;  
  long double  c; 
  ptrdiff_t  i;
  ptrdiff_t  j;
  ptrdiff_t  k;
  ptrdiff_t  q;  
  ptrdiff_t  qcheck; 
  ptrdiff_t  primroot; 
  long double pi;
  long double Kummer;  
  long double realcorrection;
  ptrdiff_t NUM_POINTS;
  ptrdiff_t nhalf;

  pi=3.1415926535897932384626433832795028841971693993751;
  
  printf("********** A. LANGUASCO  *********** \n");
  printf("Computation of the Kummer ratio [LONG DOUBLE PRECISION - GURU64]\n");
  printf("with generalised Bernoulli numbers. [DECIMATED in FREQUENCY]\n");
  
   if (argc != 2) {
  fprintf(stderr, "USAGE: %s <prime_number_q>\n", argv[0]);
  exit(1);
  }
  
  q = strtol(argv[1], NULL,10);  
   
  FILE *myFile;
  myFile = fopen("primroot.res", "r");
  
  if (myFile == NULL){
  printf("\n ******  ERROR: File primroot.res doesn't exist\n");
  exit (0);
  }
  
  /*checking the value of q in input = the one for which we got the primitive root */
  fscanf(myFile, "%tu", &qcheck);
  if (qcheck != q) {printf("\n%s", "******* ERROR: primroot of = ");
      printf("%td\n",qcheck); 
      exit (0);
      };
      
  /*getting the value of the primitive root g mod q used in the precomputation */
  fscanf(myFile, "%tu", &primroot);
  fclose(myFile);
 
  
  /* fftw initialisations */
  
  NUM_POINTS = (ptrdiff_t) q-1;   
  nhalf=NUM_POINTS/2;  
  

  long double res; /* to store the result of the final computation */
 
   
  clock_t startplan = clock();  /* start clock plan generation */
  
  /* checking if there are fftw wisdoms available */
   if (!fftwl_import_wisdom_from_filename("particular-wisdomFFTWL_ESTIMATE.txt") && !fftwl_import_wisdom_from_filename("generic-wisdomFFTWL_ESTIMATE.txt"))
  {  fprintf(stderr, "*** NOT available particular wisdom FFTWL_ESTIMATE file \n");
  fprintf(stderr, "*** NOT available generic wisdom FFTWL_ESTIMATE file \n");
  }
  if (!fftwl_import_wisdom_from_filename("particular-wisdomFFTWL_ESTIMATE.txt") && fftwl_import_wisdom_from_filename("generic-wisdomFFTWL_ESTIMATE.txt"))
  {  fprintf(stderr, "*** NOT available particular wisdom FFTWL_ESTIMATE file \n");
  fprintf(stderr, "*** Using a generic wisdom FFTWL_ESTIMATE file \n");
  }
  
  if (fftwl_import_wisdom_from_filename("particular-wisdomFFTWL_ESTIMATE.txt"))
  {  fprintf(stderr, "*** Using a particular wisdom FFTWL_ESTIMATE file \n");
  }
  
  /* guru64 plan preparation */ /* This is a BACKWARD transform; corresponds to the sum over a of  (chi)(a) */
 
  
 
  fftwl_complex *akvalues_trasf;
  akvalues_trasf = (fftwl_complex*) fftwl_malloc(sizeof(fftwl_complex) * nhalf); /*for INPLACE transforms*/ 
  
  int rank=1;
  fftw_iodim64 *ak_dims=malloc(rank*sizeof(fftw_iodim64));
	ak_dims[0].n= nhalf;  
	ak_dims[0].is=1;
  ak_dims[0].os=1; 
  int howmany_rank=1;
  fftw_iodim64 *howmany_dims=malloc(howmany_rank*sizeof(fftw_iodim64)); 
  howmany_dims[0].n=1;
  howmany_dims[0].is=1;
  howmany_dims[0].os=1; 
  fftwl_plan plan_ak = fftwl_plan_guru64_dft(rank, /*rank=1 means monodimensional vector*/
									  ak_dims, /* nhalf */
									  howmany_rank,
									  howmany_dims,
									  akvalues_trasf,  /*input for INPLACE transforms*/ 
									  akvalues_trasf,   /*output; for INPLACE transforms*/ 
									  FFTW_BACKWARD, /* 1 because it is the sign of BACKWARD transform */
									  FFTW_ESTIMATE); /*flags*/
									  
  if (plan_ak ==NULL){fprintf(stderr," ak plan creation failed\n");exit(1);}
  
  /* saving a generated wisdom for future runs */
  if (!fftwl_import_wisdom_from_filename("particular-wisdomFFTWL_ESTIMATE.txt"))
  {  fprintf(stderr, "*** Saving on file a particular wisdom \n");
  fftwl_export_wisdom_to_filename("particular-wisdomFFTWL_ESTIMATE.txt");
  }
  clock_t endplan = clock(); 
  /* stop clock plan generation */

  /*END GURU 64 PLAN PREPARATION */ 
   
  
  printf("%s", "Building the sequence g^k%q [decimated in frequency]: ");  
  
  clock_t startffttime = clock();  /* start clock */
  long double twist;
  long double angle; 
  twist = 0.0L;
  angle=(long double) pi/nhalf; 
  a=1;
  c = (long double) (a-NUM_POINTS); /*2*1-q*/
  akvalues_trasf[0][REAL]= c; /*g^0- g^((q-1)/2 */
  akvalues_trasf[0][IMAG]= 0.0L;   

  for (k = 1; k <nhalf ; k++){   
     a = (a*primroot)%q;   
     c = (long double) (2*a-q); /*g^k -R g^(k+m) \equiv g^k  -R (q-1)g^k \equiv g^k -R (-g^k) \equiv g^k -R (q -g^k) = 2*g^k -R q */
     twist = twist + angle;
     akvalues_trasf[k][REAL]= c * cosl(twist);  
     akvalues_trasf[k][IMAG]= c * sinl(twist); 
     /* depends on the sign = +1 of the transform */
  } 
  
   printf("%s\n\n", "Done.");   
   

  /* compute the Fast Fourier transform of the sequence g^k%q */
  /*IT'S A BACKWARD TRANSFORM, i.e. computes sum [a * chi(a)] */
  fftwl_execute(plan_ak);
  
   clock_t endffttime = clock(); /* stop clock */ 
   
/* prntc( akvalues_trasf, nhalf);*/
   
  
  /* destroy the FFTW plan */ 
  fftwl_destroy_plan(plan_ak); 
  
  clock_t startfinalcomp = clock();  /* start clock */
  /* instead of computing the product of generalised Bernoulli numbers (normalised) over odd characters
  we compute its complex logarithm in res */
  res=0;   
  realcorrection = (logl(pi)-(1.5)*logl(q)) * nhalf; 
  /*imagcorrection = pi * nhalf*/;
   
  for (i=0;  i < nhalf; ++i) {  
    /*z=clogl(result[j][REAL] - I * result[j][IMAG]);*/ 
    /*  sum [a * chi(a)]; it's a backward transform */
    res = res + (0.5)* logl(akvalues_trasf[i][REAL]*akvalues_trasf[i][REAL]+akvalues_trasf[i][IMAG]*akvalues_trasf[i][IMAG]); /* creall(z)*/ 
    }
 
 
  /* getting the product_(\chi odd) */
  Kummer= expl(res + realcorrection);
  
  clock_t endfinalcomp = clock(); /* stop clock */ 
  
  /* freeing the fftw_malloc*/ 
  fftwl_free(akvalues_trasf);   
  fftwl_cleanup();
  
  /* output of the result */
  printf("%s\n", "*** RESULT:");
  printf("%s", "Kummer_ratio(");
  printf("%td",q);
  printf("%s", ") = ");
  printf("%Lf\n", Kummer);
  printf("\n");  
  
    
    
  printf("%s\n", "*** TIMES:");
  
   
  /* output of the plan preparation time */
  double time_elapsed_in_sec = (endplan - startplan)/(double)CLOCKS_PER_SEC; 
  
  double min=floor(time_elapsed_in_sec/60);
  double sec=floor(time_elapsed_in_sec-min*60);
  double msec=(time_elapsed_in_sec-min*60-sec)*1000;
  printf("%s", "Plan and wisdom preparation time: ");
  printf("%f", min);
  printf("%s", " min. ");
  printf("%f", sec);
  printf("%s", " sec. ");
  printf("%f", msec);
  printf("%s\n", " msec. ");  
    
   /* output of the fft computation time */
  time_elapsed_in_sec = (endffttime - startffttime)/(double)CLOCKS_PER_SEC; 
  
  min=floor(time_elapsed_in_sec/60);
  sec=floor(time_elapsed_in_sec-min*60);
  msec=(time_elapsed_in_sec-min*60-sec)*1000;
  printf("%s", "FFT - computation time (I/O excluded): ");
  printf("%f", min);
  printf("%s", " min. ");
  printf("%f", sec);
  printf("%s", " sec. ");
  printf("%f", msec);
  printf("%s\n", " msec. "); 
  
  /* output of the final computation time */
  time_elapsed_in_sec = (endfinalcomp - startfinalcomp)/(double)CLOCKS_PER_SEC; 
  
  min=floor(time_elapsed_in_sec/60);
  sec=floor(time_elapsed_in_sec-min*60);
  msec=(time_elapsed_in_sec-min*60-sec)*1000;
  printf("%s", "Final comp Kq - time: ");
  printf("%f", min);
  printf("%s", " min. ");
  printf("%f", sec);
  printf("%s", " sec. ");
  printf("%f", msec);
  printf("%s\n", " msec. "); 
  
    /* total time */
  time_elapsed_in_sec = (endplan - startplan + endffttime - startffttime + endfinalcomp - startfinalcomp)/(double)CLOCKS_PER_SEC; 
  
  min=floor(time_elapsed_in_sec/60);
  sec=floor(time_elapsed_in_sec-min*60);
  msec=(time_elapsed_in_sec-min*60-sec)*1000;
  printf("%s", "Total elapsed time: ");
  printf("%f", min);
  printf("%s", " min. ");
  printf("%f", sec);
  printf("%s", " sec. ");
  printf("%f", msec);
  printf("%s\n", " msec. "); 
  printf("********** END PROGRAM  *********** \n");
  return 0;
}
 
/******  to compile
 Compilation instructions vary by platform. On the Mac and other UNIX-like systems, you should just be able to type:
 
clang -o Kummer-bernoulli-fftwl.exe Kummer-bernoulli-fftwl.c  -lfftw3l -lm   
 
gcc -o Kummer-bernoulli-fftwl.exe Kummer-bernoulli-fftwl.c  -lfftw3l -lm   (on the optiplex linux box at the dept)
 
 

 ******************/