/* Copyright (C) 2020 Alessandro Languasco */ #include #include #include #include #include #include #include #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 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; long double c; 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; long double pi; Eulergamma=0.57721566490153286060651209008240243104; pi=3.1415926535897932384626433832795028841971693993751; 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 min_{chi neq chi_0} |L'/L(1,chi)| with the T function. [LONG DOUBLE PRECISION - GURU64]\n"); printf("and computed psi values. [IN DOUBLE PRECISION]\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."); clock_t startinitseq = clock(); /* start clock */ /* building psi values sorted using the sequence g^k%q */ printf("%s", "Building the sequence g^k%q and psi(a/q)"); psivalues[0]=gsl_sf_psi( (double) 1/q ); psivalues[nhalf]= psivalues[0] + pi/tanl(pi/q); /*using the reflection formula for psi */ a=1; for (i = 1; i < nhalf; i++){ a = (a*primroot)%q; /* a = g^k%q */ c = (double) a/q; psivalues[i]=gsl_sf_psi(c); psivalues[nhalf+i]= psivalues[i] + pi/tanl(pi * c); /*using the reflection formula for psi */ } clock_t endinitseq = clock(); /* start 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 min over odd characters"); clock_t startcompodd = clock(); /* start clock */ long double M=LDBL_MAX;/* odd characters */ long double num=0; long double num1=0; long double den=0; long double quot=0; long double corr= -logl(q); for (i=1; i <= nhalf; ++i) { //* computes |(a+ib)/(c+id)|= |a+ib|/|c+id| */ j=2*i-1; k=NUM_POINTS-j; if ( j <= nhalf ) {num = Ttrasf[j][REAL] * psitrasf[j][REAL] + Ttrasf[j][IMAG]*psitrasf[j][IMAG]; den =psitrasf[j][REAL]*psitrasf[j][REAL]+psitrasf[j][IMAG]*psitrasf[j][IMAG]; num1 = Ttrasf[j][IMAG] * psitrasf[j][REAL] - Ttrasf[j][REAL]*psitrasf[j][IMAG]; quot = (corr+num/den)*(corr+num/den) + (num1/den)*(num1/den); if (quot < M) {M=quot;}; } else {num = Ttrasf[k][REAL] * psitrasf[k][REAL] + Ttrasf[k][IMAG]*psitrasf[k][IMAG]; den =psitrasf[k][REAL]*psitrasf[k][REAL]+psitrasf[k][IMAG]*psitrasf[k][IMAG]; num1 = Ttrasf[k][IMAG] * psitrasf[k][REAL] - Ttrasf[k][REAL]*psitrasf[k][IMAG]; quot = (corr+num/den)*(corr+num/den) + (num1/den)*(num1/den); if (quot < M) {M=quot;}; } /*printf("%Lf\n", res);*/ } oddresult=sqrtl(M); /* correction factors for odd characters */ clock_t endcompodd = clock(); /* stop clock */ printf("%s\n", ". Done"); printf("%s", "Computing the min over even characters"); clock_t startcompeven = clock(); /* start clock */ M=LDBL_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; if ( j <= nhalf ) {num = Ttrasf[j][REAL] * psitrasf[j][REAL] + Ttrasf[j][IMAG]*psitrasf[j][IMAG]; den =psitrasf[j][REAL]*psitrasf[j][REAL]+psitrasf[j][IMAG]*psitrasf[j][IMAG]; num1 = Ttrasf[j][IMAG] * psitrasf[j][REAL] - Ttrasf[j][REAL]*psitrasf[j][IMAG]; quot = (corr+num/den)*(corr+num/den) + (num1/den)*(num1/den); if (quot < M) {M=quot;}; } else {num = Ttrasf[k][REAL] * psitrasf[k][REAL] + Ttrasf[k][IMAG]*psitrasf[k][IMAG]; den =psitrasf[k][REAL]*psitrasf[k][REAL]+psitrasf[k][IMAG]*psitrasf[k][IMAG]; num1 = Ttrasf[k][IMAG] * psitrasf[k][REAL] - Ttrasf[k][REAL]*psitrasf[k][IMAG]; quot = (corr+num/den)*(corr+num/den) + (num1/den)*(num1/den); if (quot < M) {M=quot;}; } /*printf("%Lf\n", res);*/ } evenresult=sqrtl(M); 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 results *********************/ printf("%s\n\n", "*** RESULTS:"); /* print min_q */ printf("%s", "min_even("); printf("%td",q); printf("%s", ") = "); printf("%Lf", evenresult); printf("\n\n"); /* print min_q */ printf("%s", "min_odd("); printf("%td",q); printf("%s", ") = "); printf("%Lf", oddresult); printf("\n\n"); /* print min_q */ printf("%s", "min_{chi neq chi_0}|L'/L(1,chi)|("); printf("%td",q); printf("%s", ") = "); printf("%Lf", fminl(evenresult,oddresult)); printf("\n\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. "); /* 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 + 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: clang -o Min-T-fftwl.exe Min-T-fftwl.c -lfftw3l -lgsl -lgslcblas -lm clang -o Min-T-fftwl.exe Min-T-fftwl.c -lfftw3l -lgsl -lgslcblas -lm ******************/