/*CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                             C
C  COMPUTE CONVOLUTIONS USING GAUSS SUMS ASSOCIATED WITH DIRICHLET CHARACTERS C
C  09/26/15 (DKC) (maxima based on sigma0(x))				      C
C                                                                             C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC*/
#include <stdio.h>
#include <math.h>
#include "table2.h"
extern char *malloc();
void gauss(unsigned int n, unsigned int k, double *x, double *o);
//
unsigned int MAXN=2000000;  // maximum N, 4000000 maximum
unsigned int BEGINN=2;	// beginning N
unsigned int switchi=6; // select Dirichlet characters
//
void main() {
unsigned int N,i,index,count,maxcnt,k;
double sum3,sum4,temp;
double pi,resave,imsave,temp1;
double xi[1001*2],outp[2];
double *g;
FILE *Outfp;
Outfp = fopen("out1bzr.dat","w");
pi=3.141592654;
//
if ((switchi<3)||(switchi>19)) {
   printf("k value out of range \n");
   goto zskip;
   }
g=(double*) malloc(160000016);
if (g==NULL) {
   printf("not enough memory \n");
   goto zskip;
   }
for (i=0; i<1000; i++) {
   xi[2*i]=1.0;
   xi[2*i+1]=0.0;
   }
for (i=0; i<MAXN; i++) {
   if (switchi==3) {
      xi[0]=1.0;
      xi[2]=-1.0;
      gauss(i+1,3,xi,outp);
      }
   if (switchi==4) {
      xi[0]=1.0;
      xi[2]=0.0;
      xi[4]=-1.0;
      gauss(i+1,4,xi,outp);
      }
   if (switchi==5) {
      xi[0]=1.0;
      xi[2]=-1.0;
      xi[4]=-1.0;
      xi[6]=1.0;
      gauss(i+1,5,xi,outp);
      }
   if (switchi==6) {
      xi[0]=1.0;
      xi[2]=0.0;
      xi[4]=0.0;
      xi[6]=0.0;
      xi[8]=-1.0;
      gauss(i+1,6,xi,outp);
      }
   if (switchi==7) {
      xi[0]=1.0;
      xi[2]=1.0;
      xi[4]=-1.0;
      xi[6]=1.0;
      xi[8]=-1.0;
      xi[10]=-1.0;
      gauss(i+1,7,xi,outp);
      }
   if (switchi==8) {
      xi[0]=1.0;
      xi[2]=0.0;
      xi[4]=1.0;
      xi[6]=0.0;
      xi[8]=-1.0;
      xi[10]=0.0;
      xi[12]=-1.0;
      gauss(i+1,8,xi,outp);
      }
   if (switchi==9) {
      xi[0]=1.0;
      xi[1]=0.0;
      xi[2]=cos(1.0*pi/3.0);
      xi[3]=sin(1.0*pi/3.0);
      xi[4]=0.0;
      xi[5]=0.0;
      xi[6]=cos(2.0*pi/3.0);
      xi[7]=sin(2.0*pi/3.0);
      xi[8]=-xi[6];
      xi[9]=-xi[7];
      xi[10]=0.0;
      xi[11]=0.0;
      xi[12]=-xi[2];
      xi[13]=-xi[3];
      xi[14]=-1.0;
      xi[15]=0.0;
      gauss(i+1,9,xi,outp);
      }
   if (switchi==10) {
      xi[0]=1.0;
      xi[2]=0.0;
      xi[4]=-1.0;
      xi[6]=0.0;
      xi[8]=0.0;
      xi[10]=0.0;
      xi[12]=-1.0;
      xi[14]=0.0;
      xi[16]=1.0;
      gauss(i+1,10,xi,outp);
      }
   if (switchi==11) {
      xi[0]=1.0;
      xi[2]=-1.0;
      xi[4]=1.0;
      xi[6]=1.0;
      xi[8]=1.0;
      xi[10]=-1.0;
      xi[12]=-1.0;
      xi[14]=-1.0;
      xi[16]=1.0;
      xi[18]=-1.0;
      gauss(i+1,11,xi,outp);
      }
   if (switchi==12) {
      xi[0]=1.0;
      xi[2]=0.0;
      xi[4]=0.0;
      xi[6]=0.0;
      xi[8]=-1.0;
      xi[10]=0.0;
      xi[12]=-1.0;
      xi[14]=0.0;
      xi[16]=0.0;
      xi[18]=0.0;
      xi[20]=1.0;
      gauss(i+1,12,xi,outp);
      }
   if (switchi==13) {
      xi[0]=1.0;
      xi[2]=-1.0;
      xi[4]=1.0;
      xi[6]=1.0;
      xi[8]=-1.0;
      xi[10]=-1.0;
      xi[12]=-1.0;
      xi[14]=-1.0;
      xi[16]=1.0;
      xi[18]=1.0;
      xi[20]=-1.0;
      xi[22]=1.0;
      gauss(i+1,13,xi,outp);
      }
   if (switchi==14) {
      xi[0]=1.0;
      xi[2]=0.0;
      xi[4]=-1.0;
      xi[6]=0.0;
      xi[8]=-1.0;
      xi[10]=0.0;
      xi[12]=0.0;
      xi[14]=0.0;
      xi[16]=1.0;
      xi[18]=0.0;
      xi[20]=1.0;
      xi[22]=0.0;
      xi[24]=-1.0;
      gauss(i+1,14,xi,outp);
      }
   if (switchi==15) {
      xi[0]=1.0;
      xi[2]=1.0;
      xi[4]=0.0;
      xi[6]=1.0;
      xi[8]=0.0;
      xi[10]=0.0;
      xi[12]=-1.0;
      xi[14]=1.0;
      xi[16]=0.0;
      xi[18]=0.0;
      xi[20]=-1.0;
      xi[22]=0.0;
      xi[24]=-1.0;
      xi[26]=-1.0;
      gauss(i+1,15,xi,outp);
      }
   if (switchi==16) {
      xi[0]=1.0;
      xi[1]=0.0;
      xi[2]=0.0;
      xi[3]=0.0;
      xi[4]=0.0;
      xi[5]=1.0;
      xi[6]=0.0;
      xi[7]=0.0;
      xi[8]=0.0;
      xi[9]=-1.0;
      xi[10]=0.0;
      xi[11]=0.0;
      xi[12]=-1.0;
      xi[13]=0.0;
      xi[14]=0.0;
      xi[15]=0.0;
      xi[16]=-1.0;
      xi[17]=0.0;
      xi[18]=0.0;
      xi[19]=0.0;
      xi[20]=0.0;
      xi[21]=-1.0;
      xi[22]=0.0;
      xi[23]=0.0;
      xi[24]=0.0;
      xi[25]=1.0;
      xi[26]=0.0;
      xi[27]=0.0;
      xi[28]=1.0;
      xi[29]=0.0;
      gauss(i+1,16,xi,outp);
      }
   if (switchi==17) {
      xi[0]=1.0;
      xi[2]=1.0;
      xi[4]=-1.0;
      xi[6]=1.0;
      xi[8]=-1.0;
      xi[10]=-1.0;
      xi[12]=-1.0;
      xi[14]=1.0;
      xi[16]=1.0;
      xi[18]=-1.0;
      xi[20]=-1.0;
      xi[22]=-1.0;
      xi[24]=1.0;
      xi[26]=-1.0;
      xi[28]=1.0;
      xi[30]=1.0;
      gauss(i+1,17,xi,outp);
      }
   if (switchi==18) {
      xi[0]=1.0;
      xi[1]=0.0;
      xi[2]=0.0;
      xi[3]=0.0;
      xi[4]=0.0;
      xi[5]=0.0;
      xi[6]=0.0;
      xi[7]=0.0;
      xi[8]=cos(1.0*pi/3.0);
      xi[9]=sin(1.0*pi/3.0);
      xi[10]=0.0;
      xi[11]=0.0;
      xi[12]=cos(2.0*pi/3.0);
      xi[13]=sin(2.0*pi/3.0);
      xi[14]=0.0;
      xi[15]=0.0;
      xi[16]=0.0;
      xi[17]=0.0;
      xi[18]=0.0;
      xi[19]=0.0;
      xi[20]=-xi[12];
      xi[21]=-xi[13];
      xi[22]=0.0;
      xi[23]=0.0;
      xi[24]=-xi[8];
      xi[25]=-xi[9];
      xi[26]=0.0;
      xi[27]=0.0;
      xi[28]=0.0;
      xi[29]=0.0;
      xi[30]=0.0;
      xi[31]=0.0;
      xi[32]=-1.0;
      xi[33]=0.0;
      gauss(i+1,18,xi,outp);
      }
   if (switchi==19) {
      xi[0]=1.0;
      xi[2]=-1.0;
      xi[4]=-1.0;
      xi[6]=1.0;
      xi[8]=1.0;
      xi[10]=1.0;
      xi[12]=1.0;
      xi[14]=-1.0;
      xi[16]=1.0;
      xi[18]=-1.0;
      xi[20]=1.0;
      xi[22]=-1.0;
      xi[24]=-1.0;
      xi[26]=-1.0;
      xi[28]=-1.0;
      xi[30]=1.0;
      xi[32]=1.0;
      xi[34]=-1.0;
      gauss(i+1,19,xi,outp);
      }
   g[2*i]=outp[0];
   g[2*i+1]=outp[1];
   }
if ((switchi==6)||(switchi==10)||(switchi==14)||(switchi==18))
   k=switchi/2;
else
   k=switchi;
//
// compute convolutions
//
maxcnt=0;
for (N=BEGINN; N<=MAXN; N++) {
   resave=g[2*N-2];
   imsave=g[2*N-1];
   temp1=resave*resave+imsave*imsave;
   temp=temp1-(double)k;
   if (temp<0.0)
      temp=-temp;
   if (temp>0.1)
      continue;
   sum3=resave*resave;
   sum4=imsave*imsave;
   count=1;
   for (i=2; i<=N; i++) {
      index=N/i;
      if (N!=(index*i))
	 continue;
      count=count+1;
      outp[0]=g[2*index-2];
      outp[1]=g[2*index-1];
      sum3=sum3+outp[0]*outp[0];
      sum4=sum4+outp[1]*outp[1];
      }
   temp=sum3+sum4;
   if (count>=maxcnt) {
      maxcnt=count;
      printf(" %d %e %e %e %d %e %e %e\n",N,temp,sum3,sum4,count,resave,imsave,temp1);
      fprintf(Outfp," %e, %e, %e, %e, \n",(double)N,temp,(double)count,temp1);
      }
   }
zskip:
fclose(Outfp);
return;
}