/*
  Copyright (C) 2021 Alessandro Languasco
*/ 

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

#ifdef DBL_DECIMAL_DIG
  #define OP_DBL_Digs (DBL_DECIMAL_DIG)
#else  
  #ifdef DECIMAL_DIG
    #define OP_DBL_Digs (DECIMAL_DIG)
  #else  
    #define OP_DBL_Digs (DBL_DIG + 3)
  #endif
#endif

#define REAL 0
#define IMAG 1


 
int main() {
  ptrdiff_t a;  
  long double c; 
  long double d; 
  ptrdiff_t i;
  ptrdiff_t j;
  ptrdiff_t k;
  ptrdiff_t m;
  ptrdiff_t q;   
  ptrdiff_t primroot; 
  ptrdiff_t NUM_POINTS;
  ptrdiff_t nhalf;
  ptrdiff_t nquarter; 
  long double Minimumoddresult;
  long double Minimumevenresult; 
  long double Maximumoddresult;
  long double Maximumevenresult;  
  long double pi;   
  long double Eulergamma;   
  Eulergamma=0.57721566490153286060651209008240243104;
  pi= M_PI;
  long double logpi = logl(pi);
  
  printf("********** A. LANGUASCO  *********** \n");
  printf("Computation of the maxmin_{chi neq chi_0} |L(1,chi)| with the logGamma function. \n");
  printf("[Long Double PRECISION - GURU64 - OUTPLACE transforms]\n");
  printf("Acquiring q from file: q = ");
  FILE *myFile;
  myFile = fopen("primroot.res", "r");
  
  if (myFile == NULL){
    printf("File primroot.res doesn't exist\n");
    exit (0);
  }
  fscanf(myFile, "%tu", &q);
  printf("%td\n",q);
  
 /*getting the value of the primitive root g mod q used in the precomputation */       
  fscanf(myFile, "%tu", &primroot);
  
  fclose(myFile);
  
  NUM_POINTS=q-1;
  nhalf=NUM_POINTS/2;  
  nquarter=nhalf/2;  

  /*BEGIN GURU 64 PLAN PREPARATION */
  printf("%s", "Start FFT plan_lngamma_even (even) preparation\n"); 
  clock_t startplan_lngamma_even = clock();  /* start clock */
   
  /***********  LOGGAMMA PLAN PREPARATION ***********/
    /***** for an OUTPLACE TRANSFORM *****/
  long double *lngammavalues;
  lngammavalues = (long double*) fftwl_malloc(sizeof(long double) * nhalf); 

  fftwl_complex *lngammatrasf;
  lngammatrasf = (fftwl_complex*) fftwl_malloc(sizeof(fftwl_complex) * (nquarter+1));
  
  /***** for an INPLACE TRANSFORM *****/
 /* for an INPLACE TRANSFORM; lngammavalues should be of the same dim of the output ; 2*(nquarter+1) */
  
  /* long double *lngammavalues; */
  /*lngammavalues = (long double*) fftwl_malloc(sizeof(long double) * nhalf+2); */
  /* the FFT of a real sequence with length = nhalf has nquarter + 1 complex positions = nhalf+2 long double positions */
  /* fftwl_complex *lngammatrasf;  */
  /* lngammatrasf = (fftwl_complex*)lngammavalues;  */


  int rank=1;
  fftw_iodim64 *dims=malloc(rank*sizeof(fftw_iodim64));
  dims[0].n= nhalf;  /* length of the trasform */
  dims[0].is=1;
  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_lngamma_even = fftwl_plan_guru64_dft_r2c(rank, /*rank=1 means monodimensional vector*/
              dims, /* NUM_POINTS */
              howmany_rank,
              howmany_dims,
              lngammavalues,  /*input*/
              lngammatrasf,  /*output; for outplace transforms*/
              FFTW_ESTIMATE); /*flags*/
              
  if (plan_lngamma_even ==NULL){fprintf(stderr,"loggamma plan creation failed\n");exit(1);}
  clock_t endplan_lngamma_even = clock(); /* stop clock */
  printf("%s", "End FFT plan-lngamma preparation\n"); 
  

  
  /*********** acquire the values to generate sequence g^k%q ***********/
  
  printf("%s", "Generating log(gamma(a/q)) [(even) decimated in frequency]"); 
  
  
  /******** initialisation for log gamma (even)  **********/
  
 
  clock_t startinitseq_lngamma = clock();  /* start clock */

  a = 1; 
  for (k = 0; k <nhalf ; k++){  
     d = (long double) a/q ;
     lngammavalues[k]= logpi - logl(sinl(pi*d));  /* values of log(gamma(a_k/q) via the reflection formula*/
     a = (a*primroot)%q;   
  } 
   clock_t endinitseq_lngamma = clock();  /* stop clock */
  printf("%s\n", ". Done"); 
  
  printf("%s", "Start first FFT transform (logGamma even)"); 

  clock_t start = clock();  /* start clock */
  /* transform loggamma */
  fftwl_execute(plan_lngamma_even);

  /* lngammatrasf contains the DFT of loggamma; due to the simmetries the output has half elements than the input */
  /* lngammatrasf[i][REAL]=lngammatrasf[nhalf-i][REAL];  lngammatrasf[i][nhalf]=-lngammatrasf[nhalf-i][IMAG] */
    /*   SO THEY HAVE THE SAME ABSOLUTE VALUES  */
  
  clock_t end = clock(); /* stop clock */
  printf("%s\n", " -- End first FFT transform."); 
  
  /* destroy the FFTW plans */  
  fftwl_destroy_plan(plan_lngamma_even);  /* DO NOT free lngammatrasf; we need it later */
 
  
  /* Computing maxmin_{chi\ne chi_0} \vert L(1,chi) \vert =: M(q) */
  /* sum starts from i=1 to avoid the trivial character contribution */
  
  printf("%s", "Computing the maxmin over even nontrivial characters:"); 
  
  clock_t startcompeven = clock(); /* start clock */
  long double Minimum=LDBL_MAX;/* even characters */
  long double Maximum=0; 
  long double value=0;
  long int jminL = 0;
  long int jmaxL = 0;
          
  for (i=1;  i < nquarter; ++i) { /* computes |(a+ib)| */ 
         value = lngammatrasf[i][REAL] * lngammatrasf[i][REAL] +  lngammatrasf[i][IMAG]*  lngammatrasf[i][IMAG]; 
          if (value < Minimum) {Minimum=value; jminL = 2*i;}; 
          if (value > Maximum) {Maximum=value; jmaxL = 2*i;}; 
    /*printf("%Lf\n", M);*/
  }
  long double radq= sqrtl(q);
  Minimumevenresult=2*sqrtl(Minimum)/radq; /* min for even characters */
  Maximumevenresult=2*sqrtl(Maximum)/radq; /* max for even characters */

  printf("\n***\n");
  printf("indexes for minimal and maximal even characters |L (1,chi)| \n ");
  printf("prime = ");
  printf("%ld",q);
  printf("; primitive root = ");
  printf("%ld",primroot);
  printf("; chiminevenL index = ");
  printf("%ld",jminL);
  printf("; chimaxevenL index = ");
  printf("%ld",jmaxL);
  printf("\n***\n");  
  clock_t endcompeven = clock(); /* stop clock */
  printf("%s\n", "Done");

  /* freeing the fftw_malloc*/ 
  fftwl_free(lngammatrasf);  
  
    printf("%s", "Start FFT plan_ak preparation\n"); 
   clock_t startplan_ak = clock();  /* start clock */
  /********** ak PLAN PREPARATION ***********/
   /***** for an OUTPLACE TRANSFORM *****/
   
   fftwl_complex *akvalues;
  akvalues = (fftwl_complex*) fftwl_malloc(sizeof(fftwl_complex) * nhalf); 
  fftwl_complex *akvalues_trasf;
    /* the FFT of a complex sequence with length = nhalf has nhalf complex positions */
  akvalues_trasf = (fftwl_complex*) fftwl_malloc(sizeof(fftwl_complex) * nhalf); 
  
    /***** for an INPLACE TRANSFORM *****/
     /* This is a FORWARD transform; corresponds to the sum over a of conj(chi)(a) */
 
    /* for an INPLACE TRANSFORM;   */
    /*fftwl_complex *akvalues_trasf;*/
    /*akvalues_trasf = (fftwl_complex*) fftwl_malloc(sizeof(fftwl_complex) * nhalf);*/
  
  fftw_iodim64 *ak_dims=malloc(rank*sizeof(fftw_iodim64));
  ak_dims[0].n= nhalf;   /* length of the trasform */
  ak_dims[0].is=1;
  ak_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,  input OUTPLACE*/
              /* akvalues_trasf,   output OUTPLACE*/
              akvalues_trasf,  /*input INPLACE*/
              akvalues_trasf,  /*output INPLACE*/
              FFTW_FORWARD, /* -1 because it is the sign of FORWARD transform */
              FFTW_ESTIMATE); /*flags*/
              
  if (plan_ak ==NULL){fprintf(stderr,"ak plan creation failed\n");exit(1);}
    clock_t endplan_ak = clock(); 

  /*END GURU 64 PLAN PREPARATION */
      printf("%s", "End FFT plan-ak preparation\n"); 
    
    
      printf("%s", "Generating g^k%q [(odd) decimated in frequency]");   
   /******** initialisation for  a_k **********/
  
  /*  COMMENT: the speediness of the DIF for a_k can be improved
  by replacing the nhalf calls to sin and cos with the summation
  formulae for the sin and cos functions (so using just one call of
  sin and cos and about nhalf sums and products). This  works nicely in
  quadruple precision but in the long double precision it introduces
  some computation errors in the final value of the euler-kronecker constants.
  So we are using here the slower but more precise version of the DIF
  for the a_k sequence.
  */
  clock_t startinitseq_ak = clock();  /* start clock */
  long double sinvalue;
  long double cosvalue;
  long double twist = 0.0L;
  long double angle = (long double) pi/nhalf;
  a = 1; 
  for (k = 0; k <nhalf ; k++){  
	 c = (long double) (2*a-q);
     c = (long double) c/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 */
     /* sincosl(twist, &sinvalue, &cosvalue);  works in LINUX; not in C standard; slightly faster */
     cosvalue = cosl(twist);
     sinvalue = sinl(twist); 
     akvalues_trasf[k][REAL]= c * cosvalue;  
     akvalues_trasf[k][IMAG]= (-1.0L) * c * sinvalue; /* depends on the sign = -1 of the transform */ 
     a = (a*primroot)%q;  
     twist = twist + angle;
  } 
   clock_t endinitseq_ak = clock();  /* stop clock */
  printf("%s\n", ". Done"); 
  
    printf("%s", "Start second FFT transform (g^k%q)"); 
    clock_t start1 = clock();  /* start clock */
  /* transform ak */
  fftwl_execute(plan_ak);

  /*akvalues_trasf contains the DFT of ak; it's a complex to complex transform */
  
  clock_t end1 = clock(); /* stop clock */

  printf("%s\n", " -- End second FFT transform."); 
  
  /* destroy the ak plan */ 
  fftwl_destroy_plan(plan_ak); 
  

  printf("%s", "Computing the maxmin over odd characters:");
  
    clock_t startcompodd = clock(); /* start clock */
  
  Minimum=LDBL_MAX;/* odd characters */
  Maximum=0; 
  value=0;  
  jminL = 0;
  jmaxL = 0;
     
  for (i=0;  i < nhalf; ++i) { 
          value = akvalues_trasf[i][REAL] * akvalues_trasf[i][REAL] +  akvalues_trasf[i][IMAG]*akvalues_trasf[i][IMAG];
          if (value < Minimum) {Minimum=value; jminL = 2*i+1;};
          if (value > Maximum) {Maximum=value; jmaxL = 2*i+1;}; 
    /*printf("%Lf\n", M);*/
  }
  Minimumoddresult=pi*sqrtl(Minimum)/radq; /* min for odd characters */
  Maximumoddresult=pi*sqrtl(Maximum)/radq; /* max for odd characters */
 
  printf("\n***\n");
  printf("indexes for minimal and maximal odd characters |L (1,chi)| \n ");
  printf("prime = ");
  printf("%ld",q);
  printf("; primitive root = ");
  printf("%ld",primroot);  
  printf("; chiminoddL index = ");
  printf("%ld",jminL);
  printf("; chimaxoddL index = ");
  printf("%ld",jmaxL);
  printf("\n***\n");  
   
  clock_t endcompodd = clock(); /* stop clock */
  printf("%s\n", "Done");
  
 

  /* freeing the fftw_malloc*/ 
  fftwl_free(akvalues_trasf);  
  /* clean the fftw setting */ 
  fftwl_cleanup();
 
  
  /******************** output of the results *********************/ 
  printf("%s\n", "--------------");
  printf("%s\n\n", "*** RESULTS:");
   
   /* print min_q */
  printf("%s", "min_even(");
  printf("%td",q);
  printf("%s", ") = "); 
  /*("%Lf", Minimumevenresult);*/
  printf("%.*Lf\n", OP_DBL_Digs + 6, Minimumevenresult);

     
      /* print min_q */
  printf("%s", "min_odd(");
  printf("%td",q);
  printf("%s", ") = ");
  /*printf("%Lf", Minimumoddresult);*/
  printf("%.*Lf\n", OP_DBL_Digs + 6, Minimumoddresult); 
  
    /* print min_q */
  printf("%s", "min_{chi neq chi_0}|L(1,chi)|(");
  printf("%td",q);
  printf("%s", ") = "); 
  /*printf("%Lf", fminl(Minimumevenresult,Minimumoddresult));*/
  printf("%.*Lf\n", OP_DBL_Digs + 6, fminl(Minimumevenresult,Minimumoddresult));  
  
  printf("%s\n", "--------------");
    
   /* print max_q */
  printf("%s", "max_even(");
  printf("%td",q);
  printf("%s", ") = "); 
  /*printf("%Lf", Maximumevenresult);*/
  printf("%.*Lf\n", OP_DBL_Digs + 6, Maximumevenresult);
  
      /* print max_q */
  printf("%s", "max_odd(");
  printf("%td",q);
  printf("%s", ") = ");
  /*printf("%Lf", Maximumoddresult);*/
  printf("%.*Lf\n", OP_DBL_Digs + 6, Maximumoddresult);
  
    /* print max_q */
  printf("%s", "max_{chi neq chi_0}|L(1,chi)|(");
  printf("%td",q);
  printf("%s", ") = "); 
  /*printf("%Lf", fmaxl(Maximumevenresult,Maximumoddresult));*/
  printf("%.*Lf\n", OP_DBL_Digs + 6, fmaxl(Maximumevenresult,Maximumoddresult));

  printf("%s\n", "--------------");
  
 
  /********************* output of the computation times *********************/
  
  printf("%s\n", "*** TIMES:");
  double time_elapsed_in_sec = (end - start+ end1 - start1)/(double)CLOCKS_PER_SEC;   

  /* output of the FFT total computation time */
  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", "FFT - total computation time: ");
  printf("%f", min);
  printf("%s", " min. ");
  printf("%f", sec);
  printf("%s", " sec. ");
  printf("%f", msec);
  printf("%s\n", " msec. "); 
  
  /* output of the PLAN FFT total computation time */
  time_elapsed_in_sec = (endplan_lngamma_even - startplan_lngamma_even +endplan_ak - startplan_ak)/(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", "PLAN FFT - total computation time: ");
  printf("%f", min);
  printf("%s", " min. ");
  printf("%f", sec);
  printf("%s", " sec. ");
  printf("%f", msec);
  printf("%s\n", " msec. ");
  
  /* output of the total computation time */
  time_elapsed_in_sec = (endcompeven - startcompeven + endcompodd - startcompodd)/(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 maxmin-rq-(EK(q)-EK(q)+) - total computation time: ");
  printf("%f", min);
  printf("%s", " min. ");
  printf("%f", sec);
  printf("%s", " sec. ");
  printf("%f", msec);
  printf("%s\n", " msec. "); 
  
  
    /* output initialisation time */
  time_elapsed_in_sec = (endinitseq_lngamma - startinitseq_lngamma + endinitseq_ak-startinitseq_ak)/(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", "Initialisation - 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 = (endinitseq_lngamma - startinitseq_lngamma  
  + endinitseq_ak-startinitseq_ak 
  + end - start+ end1 - start1 
  + endplan_lngamma_even - startplan_lngamma_even  
  + endplan_ak - startplan_ak
  + endcompeven - startcompeven  
  + endcompodd - startcompodd)/(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\n", "TOTAL 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:
 
  
 /usr/bin/gcc  -o  MaxLMinL-fftwl-inplace.exe  MaxLMinL-fftwl-inplace.c -lfftw3l -lm  (on the optiplex linux box at the dept) 

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