/*
  Copyright (C) 2020 Alessandro Languasco
*/
#define _GNU_SOURCE
#include <time.h>
#include <fftw3.h> 
#include <stdio.h>
#include <stdlib.h>  
#include <float.h>
#include <math.h>
#include <quadmath.h>

#define REAL 0
#define IMAG 1


char bufRe[512];
char bufIm[512];

char bufPrint[512];
void printFloat128(__float128 x) {
  size_t n = quadmath_snprintf(bufPrint, sizeof(bufPrint), "%+-#*.20Qf", 92, x);
  if (n < sizeof(bufPrint))
  printf("%s\n", bufPrint);
  else
  printf("--error--: n = %ld\n", n);
} 
void getSvalues(__float128* Svalues, unsigned long long int dim) {
  FILE *myFile;
  myFile = fopen("precomp_S-DIF.res", "r");
  int i;
  int dummy;

  if (myFile == NULL){
   printf("File precomp_S-DIF.res doesn't exist\n");
   exit (0);
  }
  fscanf(myFile, "%d", &dummy);
  /* serve a 'consumare' il newline dopo l'intero dummy */
  fgets(bufRe, sizeof(bufRe), myFile);
  
  for (i = 0; i < dim; i++){
   fgets(bufRe, sizeof(bufRe), myFile);
   Svalues[i] = strtoflt128 (bufRe, NULL);
   /* printFloat128(Svalues[i]); */
  }
  /* printf("\ngetSValues() done.\n");*/

  /*for (i = 0; i < q-1; i++){
   printf("%Qg\n", Svalues[i]);
  }*/
  
  fclose(myFile);
  
}
  
void prntc(fftwq_complex* result, unsigned long long int NUM_POINTS) { /* for checking results */
  int i;
  int n;
  int width = 46;
  for (i = 0; i < NUM_POINTS; ++i) {
   /* printf("%llu\n", result[i][REAL]); */
   n = quadmath_snprintf(bufRe, sizeof(bufRe), "%+-#*.20Qf", width, result[i][REAL]);
   if ((size_t) n < sizeof bufRe)
  printf ("%s\n", bufRe);

   /* printf("%Qg\n", result[i][IMAG]); */
   n = quadmath_snprintf(bufIm, sizeof(bufIm), "%+-#*.20Qf", width, result[i][IMAG]);
   if ((size_t) n < sizeof bufIm)
  printf ("%s\n", bufIm);

  }
}
void prntr(__float128* result,unsigned long long int NUM_POINTS) { /* for checking results */
  int i;
  int n;
  int width = 46;
  for (i = 0; i < NUM_POINTS; ++i) {
   n = quadmath_snprintf(bufRe, sizeof(bufRe), "%+-#*.20Qf", width, result[i]);
   if ((size_t) n < sizeof bufRe)
  printf ("%s\n", bufRe);
  }
}


 
int main() {
  ptrdiff_t a;  
  __float128 c; 
  __float128 d; 
  ptrdiff_t i;
  ptrdiff_t j;
  ptrdiff_t k;
  ptrdiff_t m;
  ptrdiff_t q;  
  ptrdiff_t qcheck; 
  ptrdiff_t primroot; 
  ptrdiff_t NUM_POINTS;
  ptrdiff_t nhalf;
  ptrdiff_t nquarter; 
  __float128 oddresult;
  __float128 evenresult;
  __float128 res;
  __float128 Eulergamma;
  __float128 pi;  
  Eulergamma=0.57721566490153286060651209008240243104;
  pi=M_PIq;

 
  
  printf("********** A. LANGUASCO  *********** \n");
  printf("Computation of the min_{chi neq chi_0} |L'/L(1,chi)| with the S function. [Quadruple PRECISION - GURU64]\n");
  printf("Acquiring q from file: q = ");
  FILE *myFile;
  myFile = fopen("precomp_S-DIF.res", "r");
  
  if (myFile == NULL){
    printf("File precomp_S-DIF.res doesn't exist\n");
    exit (0);
  }
  fscanf(myFile, "%tu", &q);
  printf("%td\n",q);
  fclose(myFile);
  
  NUM_POINTS=q-1;
  nhalf=NUM_POINTS/2;  
  nquarter=nhalf/2;  

  /*BEGIN GURU 64 PLAN PREPARATION */
  printf("%s", "Start FFT plans preparation"); 
  clock_t startplan = clock();  /* start clock */
   
  /***********  LOGGAMMA PLAN PREPARATION ***********/
  /* __float128 *lngammavalues;
  lngammavalues = (__float128*) fftwq_malloc(sizeof(__float128) * NUM_POINTS); 

  fftwq_complex *lngammatrasf;
  lngammatrasf = (fftwq_complex*) fftwq_malloc(sizeof(fftwq_complex) * (nhalf+1));*/
  
  /***** for an INPLACE TRANSFORM *****/
  __float128 *lngammavalues;
  lngammavalues = (__float128*) fftwq_malloc(sizeof(__float128) * NUM_POINTS+2);  /* for an INPLACE TRANSFORM; lngammavalues should be of the same dim of the output ; 2*(nhalf+1) */
  /* the FFT of a real sequence with length = NUM_POINTS has nhalf + 1 complex positions = NUM_POINTS+2 __float128 positions*/
  fftwq_complex *lngammatrasf;
  lngammatrasf = (fftwq_complex*)lngammavalues;  

  int rank=1;
  fftw_iodim64 *dims=malloc(rank*sizeof(fftw_iodim64));
  dims[0].n= NUM_POINTS;  /* 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;  
  fftwq_plan plan_lngamma = fftwq_plan_guru64_dft_r2c(rank, /*rank=1 means monodimensional vector*/
              dims, /* NUM_POINTS */
              howmany_rank,
              howmany_dims,
              lngammavalues,  /*input*/
              lngammatrasf,  /*output IN PLACE*/  
              FFTW_ESTIMATE); /*flags*/
              
  if (plan_lngamma ==NULL){fprintf(stderr,"loggamma plan creation failed\n");exit(1);}
  
  /********** ak PLAN PREPARATION ***********/
  
    /***** for an INPLACE TRANSFORM *****/
  /*fftwq_complex *akvalues;
  akvalues = (fftwq_complex*) fftwq_malloc(sizeof(fftwq_complex) * nhalf); */
  fftwq_complex *akvalues_trasf;
   /* the FFT of a complex sequence with length = nhalf has nhalf complex positions*/
  akvalues_trasf = (fftwq_complex*) fftwq_malloc(sizeof(fftwq_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;  
  fftwq_plan plan_ak = fftwq_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_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);}
  
  /*********** S PLAN PREPARATION ***********/
  
    /***** for an INPLACE TRANSFORM *****/
  /*__float128 *Svalues;
  Svalues = (__float128*) fftwq_malloc(sizeof(__float128) * nhalf); 
  fftwq_complex *Strasf;
  Strasf = (fftwq_complex*) fftwq_malloc(sizeof(fftwq_complex) * (nquarter+1));*/
  
  __float128 *Svalues;
  Svalues = (__float128*) fftwq_malloc(sizeof(__float128) * nhalf+2);  /* for an INPLACE TRANSFORM; Svalues should be of the same dim of the output ; 2*(nquarter+1) */
   /* the FFT of a real sequence with length = nhalf has nquarter + 1 complex positions = nhalf+2 __float128 positions*/
  fftwq_complex *Strasf;
  Strasf = (fftwq_complex*)Svalues;  /* for an INPLACE TRANSFORM */
  
  fftw_iodim64 *S_dims=malloc(rank*sizeof(fftw_iodim64));
  S_dims[0].n= nhalf;   /* length of the trasform */
  S_dims[0].is=1;
  S_dims[0].os=1; 
  fftwq_plan plan_S = fftwq_plan_guru64_dft_r2c(rank, /*rank=1 means monodimensional vector*/
              S_dims, /* nhalf */
              howmany_rank,
              howmany_dims,
              Svalues,  /*input*/
              Strasf,  /*output*/ 
              FFTW_ESTIMATE); /*flags*/
              
  if (plan_S ==NULL){fprintf(stderr,"S plan creation failed\n");exit(1);}
  
  clock_t endplan = clock(); /* stop clock */
  printf("%s\n", ". Done"); 
  
  /*END GURU 64 PLAN PREPARATION */
  
  
  /*********** acquire the values to generate sequence g^k%q ***********/
  
  printf("%s", "Generating log(gamma(a/q)) and g^k%q [decimated in frequency]"); 
  myFile = fopen("primroot.res", "r");
  
  if (myFile == NULL){
    printf("\n ******  ERROR: File primroot.res doesn't exist\n");
    exit (0);
  }
  /*checking if it is the correct q */
  fscanf(myFile, "%tu", &qcheck);

  if (qcheck != q) {printf("\n%s", "******* ERROR: These is the generator  for q = ");
        printf("%td\n",qcheck); 
        exit (0);
        };
  /*getting the value of the primitive root g mod q used in the precomputation */       
  fscanf(myFile, "%tu", &primroot);
  
  /*printf("%lu",primroot);*/
  
  /******** initialisation for log gamma and 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 unfortunately 
  introduces some computation errors in the final value.
  So we are using here the slower but more precise version of the DIF
  for the a_k sequence.
  */ 
     clock_t startinitseq = clock();  /* start clock */
  __float128 sinvalue;
  __float128 cosvalue;
  __float128 logpi = logq(pi);
  __float128 twist = 0.0q;
  __float128 angle = (__float128) pi/nhalf;
  a = 1; 
  for (k = 0; k <nhalf ; k++){  
     c = (__float128) (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 */
     /* sincosq(twist, &sinvalue, &cosvalue);  works in LINUX; not in C standard; slightly faster */
     cosvalue = cosq(twist);
     sinvalue = sinq(twist); 
     akvalues_trasf[k][REAL]= c * cosvalue;  
     akvalues_trasf[k][IMAG]= (-1.0q) * c * sinvalue; /* depends on the sign = -1 of the transform */
     d = (__float128) a/q ;
     lngammavalues[k]=lgammaq(d);  /* values of log(gamma(a_k/q)*/
     lngammavalues[k+nhalf] = -lngammavalues[k] + logpi - logq(sinq(pi*d));/* values of log(gamma(a_(k+m)/q) via the reflection formula*/
     a = (a*primroot)%q;  
     twist = twist + angle;
  } 
       clock_t endinitseq = clock();  /* end clock */
  printf("%s\n", ". Done"); 
  
  printf("%s", "Start first FFT transform"); 

       clock_t start = clock();  /* start clock */
  /* transform loggamma */
  fftwq_execute(plan_lngamma);
  
  /* lngammatrasf contains the DFT of loggamma; due to the simmetries the output has half elements than the input */
  /* lngammatrasf[i][REAL]=lngammatrasf[NUM_POINTS-i][REAL];  lngammatrasf[i][IMAG]=-lngammatrasf[NUM_POINTS-i][IMAG] */
  
  clock_t end = clock(); /* stop clock */
  printf("%s\n", " -- End first FFT transform."); 
  
  /* destroy the FFTW plans */  
  fftwq_destroy_plan(plan_lngamma);  /* DO NOT free lngammatrasf; we need it later */
  
 
  printf("%s", "Start second FFT transform"); 
    clock_t start1 = clock();  /* start clock */
  /* transform ak */
  fftwq_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 */ 
  fftwq_destroy_plan(plan_ak);
  
  printf("%s", "Computing the min over odd characters"); 
  
  clock_t startcompodd = clock(); /* start clock */
  __float128 M=FLT128_MAX;/* odd characters */
  __float128 num=0;
  __float128 num1=0;
  __float128 den=0;
  __float128 quot=0;   
  __float128 corr= logq(2*pi)+Eulergamma;   
  for (i=1;  i <= nhalf; ++i) { /* computes |(a+ib)/(c+id)|= |a+ib|/|c+id| */
    j=2*i-1;
    m= i-1;
    k=NUM_POINTS-j;
    if ( j <= nhalf )
    {num = q*(lngammatrasf[j][REAL] * akvalues_trasf[m][REAL] + lngammatrasf[j][IMAG]*akvalues_trasf[m][IMAG]);
          den =akvalues_trasf[m][REAL]*akvalues_trasf[m][REAL]+akvalues_trasf[m][IMAG]*akvalues_trasf[m][IMAG];
          num1 = q*(lngammatrasf[j][IMAG] * akvalues_trasf[m][REAL] - lngammatrasf[j][REAL]*akvalues_trasf[m][IMAG]);
          quot = (corr+num/den)*(corr+num/den) + (num1/den)*(num1/den);
          if (quot< M) {M=quot;};
    }
   else
   {num = q*(lngammatrasf[k][REAL] * akvalues_trasf[m][REAL] - lngammatrasf[k][IMAG]*akvalues_trasf[m][IMAG]);
          den =akvalues_trasf[m][REAL]*akvalues_trasf[m][REAL]+akvalues_trasf[m][IMAG]*akvalues_trasf[m][IMAG];
          num1 = q*(-lngammatrasf[k][IMAG] * akvalues_trasf[m][REAL] - lngammatrasf[k][REAL]*akvalues_trasf[m][IMAG]);
          quot = (corr+num/den)*(corr+num/den) + (num1/den)*(num1/den);
          if (quot< M) {M=quot;};
    }
    /*printf("%Lf\n", M);*/
  }
  oddresult=sqrtq(M); /* min for odd characters */
  
  clock_t endcompodd = clock(); /* stop clock */
  printf("%s\n", ". Done");
  
  /* freeing the fftw_malloc*/ 
  fftwq_free(akvalues_trasf);
  
  /* Third TRANSFORM: S (decimated in frequency) */
  
  printf("%s", "Acquiring precomputed values for S(a/q) [decimated in frequency]");
  clock_t startreadonfile = clock(); /* start clock */  
  
  getSvalues(Svalues,nhalf);  /* Svalues has just nhalf meaningful entries; the remaining two are for the inplace transform*/

  clock_t endreadonfile = clock(); /* end clock */
  
  printf("%s\n", ". Done");
  
  printf("%s", "Start third FFT transform"); 
  clock_t start2 = clock();  /* start clock */
  
  fftwq_execute(plan_S);
  /* prnt(result,nhalf);  result contains the DFT of S */
  
  clock_t end2 = clock(); /* stop clock */
  printf("%s\n", " -- End third FFT transform."); 
  
   /* destroy the S plan */
  fftwq_destroy_plan(plan_S);
  
  /* Computing min_{chi\ne chi_0} \vert L'/L(1,\chi) \vert =: M(q) */
  /* sum starts from i=1 to avoid the trivial character contribution */
  
  /* Strasf contains the DFT of S; due to the simmetries the output has half elements than the input */
  /* Strasf[i][REAL]=result[nhalf-i][REAL];  Strasf[i][IMAG]=-result[nhalf-i][IMAG] */

  printf("%s", "Computing the min over even characters"); 
  clock_t startcompeven = clock(); /* start clock */
  
  M=FLT128_MAX;/* even characters */
  num=0;
  num1=0;
  den=0;
  quot=0;  
  for (i=1;  i <= nhalf-1; ++i) { 
        j=2*i;
        k=NUM_POINTS-j;
        m=nhalf-i;
        if ( i <= nquarter )
        {num = (Strasf[i][REAL] * lngammatrasf[j][REAL] +  Strasf[i][IMAG]*lngammatrasf[j][IMAG])/2;
          den =lngammatrasf[j][REAL] *lngammatrasf[j][REAL] + lngammatrasf[j][IMAG]*lngammatrasf[j][IMAG] ;
          num1 = (Strasf[i][IMAG] * lngammatrasf[j][REAL] - Strasf[i][REAL]*lngammatrasf[j][IMAG])/2;
          quot = (corr+num/den)*(corr+num/den) + (num1/den)*(num1/den);
          if (quot< M) {M=quot;};
          }
       else
       {num = (Strasf[m][REAL] * lngammatrasf[k][REAL] +  Strasf[m][IMAG]*lngammatrasf[k][IMAG])/2;
          den =lngammatrasf[k][REAL] *lngammatrasf[k][REAL] + lngammatrasf[k][IMAG]*lngammatrasf[k][IMAG] ;
          num1 = (Strasf[m][IMAG] * lngammatrasf[k][REAL] - Strasf[m][REAL]*lngammatrasf[k][IMAG])/2;
          quot = (corr+num/den)*(corr+num/den) + (num1/den)*(num1/den);
          if (quot< M) {M=quot;};
          }    
       /*printf("%Lf\n", M);*/
        }
  evenresult=sqrtq(M); /* min for even characters */
  clock_t endcompeven = clock(); /* stop clock */
  printf("%s\n\n", ". Done");
  
  /* freeing the fftw_malloc*/ 
  fftwq_free(lngammatrasf);  fftwq_free(Strasf);
  /* clean the fftw setting */ 
  fftwq_cleanup();
 
  
  /******************** output of the results *********************/
  printf("%s\n\n", "*** RESULTS:");
  
   /* print min_q */
  printf("%s", "min_even(");
  printf("%td",q);
  printf("%s", ") = ");
  {
  int width = 46;
  int nquad = quadmath_snprintf (bufRe, sizeof bufRe, "%+-#*.20Qf", width, evenresult);
  if ((size_t) nquad < sizeof bufRe)
   printf ("%s\n", bufRe);
  }
  
      /* print min_q */
  printf("%s", "min_odd(");
  printf("%td",q);
  printf("%s", ") = ");
  {
  int width = 46;
  int nquad = quadmath_snprintf (bufRe, sizeof bufRe, "%+-#*.20Qf", width, oddresult);
  if ((size_t) nquad < sizeof bufRe)
   printf ("%s\n", bufRe);
  }
  
    /* print min_q */
  printf("%s", "min_{chi neq chi_0}|L'/L(1,chi)|(");
  printf("%td",q);
  printf("%s", ") = ");
  {
  int width = 46;
  int nquad = quadmath_snprintf (bufRe, sizeof bufRe, "%+-#*.20Qf", width, fminq(evenresult,oddresult));
  if ((size_t) nquad < sizeof bufRe)
   printf ("%s\n", bufRe);
  }
   
  
  /********************* output of the computation times *********************/
  
  printf("%s\n", "*** TIMES:");
  double time_elapsed_in_sec = (end - start+ end1 - start1 + end2 -start2)/(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 - startplan)/(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 minq - 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 I/O time */
  time_elapsed_in_sec = (endreadonfile - startreadonfile)/(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", "Input data from file - 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 - startinitseq)/(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 - startinitseq + end - start+ end1 - start1 + end2 -start2+ endplan - startplan + endcompeven - startcompeven + endcompodd - startcompodd+ endreadonfile - startreadonfile)/(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:
 
/usr/bin/gcc -o Min-S-fftwq.exe  Min-S-fftwq.c -lfftw3q -lquadmath -lm  (on the optiplex linux box at the dept) 

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