/*CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C C
C COMPUTE CONVOLUTIONS USING GAUSS SUMS ASSOCIATED WITH DIRICHLET CHARACTERS C
C 09/26/15 (DKC) (maxima determined by number of positive divisors) 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=1000000; // maximum N, 4000000 maximum
unsigned int BEGINN=2; // beginning N
unsigned int switchi=5; // select Dirichlet characters
//
void main() {
unsigned int N,i,index,count,maxcnt;
double sum3,sum4,temp;
double pi,resave,imsave,temp1;
double xi[1001*2],outp[2];
double *g;
FILE *Outfp;
Outfp = fopen("out1bzq.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];
}
//
// compute convolutions
//
maxcnt=0;
for (N=BEGINN; N<=MAXN; N++) {
count=2;
for (i=2; i<N; i++) {
if (N==(N/i)*i)
count=count+1;
}
if (count>=maxcnt) {
maxcnt=count;
resave=g[2*N-2];
imsave=g[2*N-1];
sum3=resave*resave;
sum4=imsave*imsave;
for (i=2; i<=N; i++) {
index=N/i;
if (N!=(index*i))
continue;
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;
temp1=resave*resave+imsave*imsave;
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;
}