/*
  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 = 1; 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 psi values [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);
 
  
  /* fftwl 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 fftwl 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 *psivalues_trasf; 
  psivalues_trasf = (fftwl_complex*) fftwl_malloc(sizeof(fftwl_complex) * nhalf); /*for INPLACE transforms*/  
  
  int rank=1;
  fftw_iodim64 *psi_dims=malloc(rank*sizeof(fftw_iodim64));
  psi_dims[0].n= nhalf;  
  psi_dims[0].is=1;
  psi_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_psi = fftwl_plan_guru64_dft(rank, /*rank=1 means monodimensional vector*/
              psi_dims, /* nhalf */
              howmany_rank,
              howmany_dims,
              psivalues_trasf,  /*input for INPLACE transforms*/ 
              psivalues_trasf,  /*output; for INPLACE transforms*/ 
              FFTW_BACKWARD, /* 1 because it is the sign of BACKWARD transform */
              FFTW_ESTIMATE); /*flags*/
              
  if (plan_psi ==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) pi*a;
  c = (long double) c/q;
  psivalues_trasf[0][REAL] = -(1.0L)*pi/tanl(c); /* Decimated in frquency values of psi via the reflection formula */
  psivalues_trasf[0][IMAG] = 0.0L;  

  for (k = 1; k <nhalf ; k++){  
    a = (a*primroot)%q;
    c = (long double) pi*a;
    c = (long double) c/q; 
    c = -(1.0L)*pi/tanl(c); 
    twist = twist + angle;  
    /* printf("%Lf", twist);*/
    psivalues_trasf[k][REAL]= c * cosl(twist);  /* Decimated in frquency values of psi via the reflection formula */
    psivalues_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 psi(a/q) DIF */
  /*IT'S A BACKWARD TRANSFORM, i.e. computes sum L(1,chi] */
  fftwl_execute(plan_psi);
  
  
  clock_t endffttime = clock(); /* stop clock */ 
  
   
  /* destroy the FFTW plan */ 
  fftwl_destroy_plan(plan_psi);
  
  
  clock_t startfinalcomp = clock();  /* start clock */ 
 
  
  /* instead of computing the product of L(1.chi) over odd characters
  we compute its complex logarithm in res */
  res=0;  
  realcorrection = (-1.0)* nhalf * logl(q);  
  /*imagcorrection = pi * nhalf;*/
  
  for (i=0;  i < nhalf; ++i) {  
   /*  sum [a * chi(a)]; it's a backward transform */
   res = res + (0.5)* logl(psivalues_trasf[i][REAL]*psivalues_trasf[i][REAL]+psivalues_trasf[i][IMAG]*psivalues_trasf[i][IMAG]); /* creall(z)*/ 
   }
 
 
  /* getting the product_(\chi odd) */
  Kummer= expl(res + realcorrection);
  
  clock_t endfinalcomp = clock(); /* stop clock */ 
  
  /* freeing the fftwl_malloc*/ 
  fftwl_free(psivalues_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-psi-fftwl.exe Kummer-psi-fftwl.c -lfftw3l -lm 
         using long double precision (on mac ; fftwl installed using homebrew)
 
 gcc -o Kummer-psi-fftwl.exe Kummer-psi-fftwl.c -lfftw3l -lm
 ******************/