/*
    Copyright (C) 2019 Alessandro Languasco
*/

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

#define REAL 0
#define IMAG 1


void getTvalues( long double* Tvalues, ptrdiff_t q) {
    FILE *myFile;
    myFile = fopen("precomp_T.res", "r");
    ptrdiff_t i;
    ptrdiff_t dummy;
    
    
    fscanf(myFile, "%tu", &dummy);
    
    for (i = 0; i < q-1; i++){
        fscanf(myFile, "%Lf", &Tvalues[i]);
    } 
    
    fclose(myFile);
    
}

void getpsivalues(long double* psivalues, ptrdiff_t q) {
    FILE *myFile;
    myFile = fopen("precomp_psi.res", "r");
    ptrdiff_t i;
    ptrdiff_t dummy;
 
    fscanf(myFile, "%tu", &dummy);
        
    for (i = 0; i < q-1; i++){
        fscanf(myFile, "%Lf", &psivalues[i]); 
    } 
    
    fclose(myFile);
}

 

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

 

int main() {
    ptrdiff_t a; 
	ptrdiff_t i;
	ptrdiff_t j;
	ptrdiff_t k;
    ptrdiff_t q; 
    ptrdiff_t qcheck; 
    ptrdiff_t primroot; 
    ptrdiff_t NUM_POINTS;
    ptrdiff_t nhalf;
    long double Eulergamma;
    Eulergamma=0.57721566490153286060651209008240243104;
    long double res;
    long double evenresult; 
    long double oddresult;

    /* input of a prime number from keyboard */
   /* printf("Give me a prime:\n");
    scanf("%ld", &q);*/
    
    printf("********** A. LANGUASCO  *********** \n");
    printf("Computation of the Euler-Kronecker constant with the T function. [LONG DOUBLE PRECISION - GURU64]\n");
    printf("and precomputed psi values.\n");
    printf("Acquiring q from file: q = ");
    FILE *myFile;
    myFile = fopen("precomp_T.res", "r"); 
    
    if (myFile == NULL){
        printf("\n ******  ERROR: File precomp_T.res doesn't exist\n");
        exit (0);
    }
    fscanf(myFile, "%tu", &q);
    printf("%td\n",q);
    fclose(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, "%td", &primroot);
    fclose(myFile);
    
    NUM_POINTS=q-1;
    nhalf=NUM_POINTS/2;  
    
     	/*BEGIN GURU 64 PLAN PREPARATION */

    long double *Tvalues;
    Tvalues = (long double*) fftwl_malloc(sizeof(long double) * NUM_POINTS+2 ); /* for an INPLACE TRANSFORMS; values should be of the same dim of the output ; 2*(nhalf+1) */
    fftwl_complex *Ttrasf;
   /* Ttrasf = (fftwl_complex*) fftwl_malloc(sizeof(fftwl_complex) * (nhalf+1));*/
    Ttrasf = (fftwl_complex*)Tvalues;  /* for an INPLACE TRANSFORM */


 	printf("%s", "Start FFT plans preparation"); 
    clock_t startplan = clock();  /* start clock */
 	int rank=1;
	fftwl_iodim64 *dims=malloc(rank*sizeof(fftwl_iodim64));
	dims[0].n= NUM_POINTS;  
	dims[0].is=1;
    dims[0].os=1;       
    int howmany_rank=1;
    fftwl_iodim64 *howmany_dims=malloc(howmany_rank*sizeof(fftwl_iodim64)); 
    howmany_dims[0].n=1;
    howmany_dims[0].is=1;
    howmany_dims[0].os=1;   
	fftwl_plan plan_T = fftwl_plan_guru64_dft_r2c(rank, /*rank=1 means monodimensional vector*/
									  dims, /* NUM_POINTS */
									  howmany_rank,
									  howmany_dims,
                                      Tvalues,  /*input*/
                                      Ttrasf,   /*output*/ 
                                      FFTW_ESTIMATE); /*flags*/
                                      
    if (plan_T ==NULL){fprintf(stderr,"T plan creation failed\n");exit(1);}
    
    
        
    long double *psivalues;
    psivalues = (long double*) fftwl_malloc(sizeof(long double) * NUM_POINTS+2); /* for an INPLACE TRANSFORMS; values should be of the same dim of the output ; 2*(nhalf+1) */
    fftwl_complex *psitrasf;
   /* Ttrasf = (fftwl_complex*) fftwl_malloc(sizeof(fftwl_complex) * (nhalf+1));*/
    psitrasf = (fftwl_complex*)psivalues;  /* for an INPLACE TRANSFORM */
    
    fftwl_plan plan_psi = fftwl_plan_guru64_dft_r2c(rank, /*rank=1 means monodimensional vector*/
									  dims, /* NUM_POINTS */
									  howmany_rank,
									  howmany_dims,
                                      psivalues,  /*input*/
                                      psitrasf,   /*output*/ 
                                      FFTW_ESTIMATE); /*flags*/
                                      
    if (plan_psi ==NULL){fprintf(stderr," psi plan creation failed\n");exit(1);}
    
     /* the two plan differ just for the output array; for T it is not a generic one*/
    clock_t endplan = clock(); /* stop clock */
    printf("%s\n", ". Done"); 
    /*END GURU 64 PLAN PREPARATION */
    
                                      
    
    printf("%s", "Acquiring precomputed values for T(a/q)");
    clock_t startreadonfile = clock(); /* start clock */
    /* acquire Tvalues sorted using the sequence g^k%q */
    getTvalues(Tvalues,q);
    clock_t endreadonfile = clock(); /* end clock */
    printf("%s\n", ". Done");
    
    printf("%s", "Start first FFT transform");  
    clock_t start = clock();  /* start clock */
                                      
    fftwl_execute(plan_T); /*FFT of Tvalues */
    /* prnt(result,nhalf);   result contains the DFT of T*/ 

     /* destroy the FFTW plan */ 
    fftwl_destroy_plan(plan_T);
    /*prnt(Ttrasf,NUM_POINTS); */
    
 	clock_t end = clock(); /* stop clock */ 
 	 printf("%s\n", " -- End first FFT transform."); 
    
     printf("%s", "Acquiring precomputed values for psi(a/q)");
    clock_t startreadonfilepsi = clock(); /* start clock */
    /* acquire psi values sorted using the sequence g^k%q */
    getpsivalues(psivalues, q);  
    clock_t endreadonfilepsi = clock(); /* end clock */  
    printf("%s\n", ". Done");
    
 	/*prntr(values,NUM_POINTS); */
 
    printf("%s", "Start second FFT transform"); 
 	clock_t start1 = clock();  /* start clock */
     
    fftwl_execute(plan_psi);
    /* prnt(result,nhalf);  result contains the DFT of psi */
     clock_t end1 = clock(); /* stop clock */
    printf("%s\n", " -- End second FFT transform."); 
    
          /* destroy the FFTW plan */ 
    fftwl_destroy_plan(plan_psi); 
      
    /* psitrasf contains the DFT of psi; due to the simmetries the output has half elements than the input */
    /* psitrasf[i][REAL]=result[NUM_POINTS-i][REAL];  psitrasf[i][IMAG]=-result[NUM_POINTS-i][IMAG] */
    
    printf("%s", "Computing the sum over odd characters"); 
  
  clock_t startcompodd = clock(); /* start clock */
  res=0;  /* odd characters */
  for (i=1;  i <= nhalf; ++i) { /* computes the real part of (a+ib)/(c+id); the imaginary part is 0 by def */
    j=2*i-1;
    k=NUM_POINTS-j;
    if ( j <= nhalf )
    {res = res +
     (Ttrasf[j][REAL]*psitrasf[j][REAL]+Ttrasf[j][IMAG]*psitrasf[j][IMAG])
     /(psitrasf[j][REAL]*psitrasf[j][REAL]+psitrasf[j][IMAG]*psitrasf[j][IMAG]);
    }
   else
   {res = res +
    (Ttrasf[k][REAL]*psitrasf[k][REAL]+Ttrasf[k][IMAG]*psitrasf[k][IMAG])
    /(psitrasf[k][REAL]*psitrasf[k][REAL]+psitrasf[k][IMAG]*psitrasf[k][IMAG]);
    }
    /*printf("%Lf\n", res);*/
  }
  oddresult=res; /* correction factors for odd characters */
  
  clock_t endcompodd = clock(); /* stop clock */
  printf("%s\n", ". Done");
  
  
  printf("%s", "Computing the sum over even characters"); 
  clock_t startcompeven = clock(); /* start clock */
  res=0;/* even characters */
  for (i=1;  i <= nhalf-1; ++i) { 
        j=2*i;
        k=NUM_POINTS-j; 
        if ( j <= nhalf )
        { res = res +
         (Ttrasf[j][REAL]*psitrasf[j][REAL]+Ttrasf[j][IMAG]*psitrasf[j][IMAG])
         /(psitrasf[j][REAL]*psitrasf[j][REAL]+psitrasf[j][IMAG]*psitrasf[j][IMAG]);
          }
       else
       { res = res +
         (Ttrasf[k][REAL]*psitrasf[k][REAL]+Ttrasf[k][IMAG]*psitrasf[k][IMAG])
         /(psitrasf[k][REAL]*psitrasf[k][REAL]+psitrasf[k][IMAG]*psitrasf[k][IMAG]);
          }    
       /*printf("%Lf\n", res);*/
        }
  evenresult=res;
  clock_t endcompeven = clock(); /* stop clock */
  printf("%s\n\n", ". Done");
  
  
  
    /* freeing the fftw_malloc*/ 
    fftwl_free(psitrasf); fftwl_free(Ttrasf); 
	/* clean the fftw setting */ 
    fftwl_cleanup();
 
    
    /* output of the result */
    printf("%s\n\n", "*** RESULT:");
    printf("%s", "EK(");
    printf("%td",q);
    printf("%s", ") = ");
     /*  summing with the correction factor: - (q-2) logl(q) + Eulergamma */
	printf("%Lf\n", Eulergamma-logl(q)*(q-2)+ evenresult + oddresult);
    printf("\n");
    
     /* print EK_q^+ */
  printf("%s", "EK(");
  printf("%td",q);
  printf("%s", ")^+ = ");
  printf("%Lf", Eulergamma-logl(q)*(q-3)/2+ evenresult);
  printf("\n\n");
  
  /* print EK_q- EK_q^+ */
  printf("%s", "EK(");
  printf("%td",q);
  printf("%s", ") - ");
  printf("%s", "EK(");
  printf("%td",q);
  printf("%s", ")^+ = ");
  printf("%Lf\n", -logl(q)*(q-1)/2+ oddresult);
  printf("\n");
    
    printf("%s\n", "*** TIMES:");
    /* output of the computation time */
     
     double time_elapsed_in_sec = (end - start+ end1 - start1)/(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", "FFT - computation time: ");
    printf("%f", min);
    printf("%s", " min. ");
    printf("%f", sec);
    printf("%s", " sec. ");
    printf("%f", msec);
    printf("%s\n", " msec. "); 
    
     /* output of the 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 - computation time:: ");
    printf("%f", min);
    printf("%s", " min. ");
    printf("%f", sec);
    printf("%s", " sec. ");
    printf("%f", msec);
    printf("%s\n", " msec. ");
    
    /* output of the computation time */
     time_elapsed_in_sec = (endcompodd - startcompodd +endcompeven - startcompeven)/(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 EKq - 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", "Acquiring precomputed data from file - 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 = (end - start+ end1 - start1 + endplan - startplan + endcompodd - startcompodd +endcompeven - startcompeven + 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:
 
 clang -o EK-T-fftwl-psiprec.exe EK-T-fftwl-psiprec.c -lfftw3l   -lm
 
 clang -o EK-T-fftwl-psiprec.exe EK-T-fftwl-psiprec.c -lfftw3l  -lm
  

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