/*CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C C
C COMPUTE MERTENS FUNCTION C
C 05/30/14 (DKC) C
C C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC*/
#include <math.h>
unsigned int haros5(unsigned int N, unsigned int M, unsigned int *R,
unsigned int H, unsigned int K, unsigned int HP,
unsigned int KP, unsigned int D);
unsigned int lagrange1(unsigned int N, unsigned int D, unsigned int O);
double merten4b(unsigned int N, unsigned int *count, double *sum,
unsigned int horder) {
unsigned int S[250000];
unsigned int H,K,total;
unsigned int L,I,J,count1,count2,ctemp1,ctemp2;
double temp1,sum1,sum2,temp,delsq;
double pi;
if (N==1) {
count[0]=0;
count[1]=0;
sum[0]=0.5;
sum[1]=0.0;
sum[2]=0.0;
sum[3]=0.0;
return(1.0);
}
if (N==2) {
count[0]=0;
count[1]=0;
sum[0]=0.0;
sum[1]=0.0;
sum[2]=0.0;
sum[3]=0.0;
return(0.0);
}
if (N==3) {
count[0]=0;
count[1]=0;
sum[0]=0.0;
sum[1]=-0.5;
if (horder==0) {
sum[2]=0.1666667;
sum[3]=0.2357023;
}
else {
sum[2]=0.6666667;
sum[3]=0.7453558;
}
return(-1.0);
}
if (N==4) {
count[0]=0;
count[1]=1;
sum[0]=0.0;
sum[1]=-0.5;
if (horder==0) {
sum[2]=0.1666667;
sum[3]=0.2886751;
}
else {
sum[2]=0.9166666;
sum[3]=1.0507931;
}
return(-1.0);
}
//
pi=3.141592654;
if ((N>4)&&(N<9)) {
ctemp1=haros5(N,0,S,0,1,1,N,4);
ctemp1=ctemp1-1;
sum1=0.0;
for (I=1; I<ctemp1; I++) {
H=S[I]>>16; // towards 1/4
K=S[I]&0xffff;
temp1=cos(2.0*pi*(double)H/(double)K);
sum1=sum1+temp1;
}
count1=ctemp1-1;
count[0]=count1;
sum[0]=sum1;
//
L=lagrange1(1,4,N);
ctemp2=haros5(N,0,S,1,4,L>>16,L&0xffff,2);
ctemp2=ctemp2-1;
sum2=0.0;
for (I=1; I<ctemp2; I++) {
H=S[I]>>16; // towards 1/2
K=S[I]&0xffff;
temp1=cos(2.0*pi*(double)H/(double)K);
sum2=sum2+temp1;
}
count2=ctemp2-1;
count[1]=count2;
sum[1]=sum2;
//
total=count1+count2;
if (horder==0)
total=total*2+4;
else
total=total+2;
temp=0.0; // |delta| measure
delsq=0.0; // delta^2 measure
ctemp1=haros5(N,0,S,0,1,1,N,2);
for (I=1; I<ctemp1; I++) {
H=S[I]>>16; // towards 1/2
K=S[I]&0xffff;
temp1=(double)I/(double)total-(double)H/(double)K;
delsq=delsq+temp1*temp1;
if (temp1<0.0)
temp1=-temp1;
temp=temp+temp1;
}
if (horder==0) {
J=total-ctemp1+1;
for (I=ctemp1-1; I>0; I--) {
H=S[I]>>16;
K=S[I]&0xffff;
temp1=(double)J/(double)total-(double)(K-H)/(double)K;
delsq=delsq+temp1*temp1;
if (temp1<0.0)
temp1=-temp1;
temp=temp+temp1;
J=J+1;
}
}
sum[2]=temp;
delsq=delsq*(double)total;
delsq=sqrt(delsq);
sum[3]=delsq;
return(2*(sum1+sum2));
}
if (N>2560) {
count[0]=0;
count[1]=0;
sum[0]=0.0;
sum[1]=0.0;
sum[2]=0.0;
sum[3]=0.0;
return(0.0);
}
ctemp1=haros5(N,0,S,0,1,1,N,8);
sum1=0.0;
for (I=1; I<ctemp1; I++) {
H=S[I]>>16; // towards 1/8
K=S[I]&0xffff;
temp1=cos(2.0*pi*(double)H/(double)K);
sum1=sum1+temp1;
}
L=lagrange1(1,8,N);
count1=haros5(N,0,S,1,8,L>>16,L&0xffff,4);
for (I=1; I<count1; I++) {
H=S[I]>>16; // towards 1/4
K=S[I]&0xffff;
temp1=cos(2.0*pi*(double)H/(double)K);
sum1=sum1+temp1;
}
count1=count1+ctemp1-3;
count[0]=count1;
sum[0]=sum1;
//
L=lagrange1(1,4,N);
ctemp2=haros5(N,0,S,1,4,L>>16,L&0xffff,8);
sum2=0.0;
for (I=1; I<ctemp2; I++) {
H=S[I]>>16; // towards 3/8
K=S[I]&0xffff;
temp1=cos(2.0*pi*(double)H/(double)K);
sum2=sum2+temp1;
}
L=lagrange1(3,8,N);
count2=haros5(N,0,S,3,8,L>>16,L&0xffff,2);
count2=count2-1;
for (I=1; I<count2; I++) {
H=S[I]>>16; // towards 1/2
K=S[I]&0xffff;
temp1=cos(2.0*pi*(double)H/(double)K);
sum2=sum2+temp1;
}
count2=count2+ctemp2-2;
count[1]=count2;
sum[1]=sum2;
total=count1+count2;
if (horder==0)
total=total*2+4;
else
total=total+2;
//
// compute Franel/Landau measures
//
temp=0.0; // |delta| measure
delsq=0.0; // delta^2 measure
ctemp1=haros5(N,0,S,0,1,1,N,8);
for (I=1; I<ctemp1; I++) {
H=S[I]>>16; // towards 1/8
K=S[I]&0xffff;
temp1=(double)I/(double)total-(double)H/(double)K;
delsq=delsq+temp1*temp1;
if (temp1<0.0)
temp1=-temp1;
temp=temp+temp1;
}
if (horder==0) {
J=total-ctemp1+1;
for (I=ctemp1-1; I>0; I--) {
H=S[I]>>16;
K=S[I]&0xffff;
temp1=(double)J/(double)total-(double)(K-H)/(double)K;
delsq=delsq+temp1*temp1;
if (temp1<0.0)
temp1=-temp1;
temp=temp+temp1;
J=J+1;
}
}
ctemp1=ctemp1-1;
L=lagrange1(1,8,N);
count1=haros5(N,0,S,1,8,L>>16,L&0xffff,4);
for (I=1; I<count1; I++) {
H=S[I]>>16; // towards 1/4
K=S[I]&0xffff;
temp1=(double)(I+ctemp1)/(double)total-(double)H/(double)K;
delsq=delsq+temp1*temp1;
if (temp1<0.0)
temp1=-temp1;
temp=temp+temp1;
}
if (horder==0) {
J=total-count1-ctemp1+1;
for (I=count1-1; I>0; I--) {
H=S[I]>>16;
K=S[I]&0xffff;
temp1=(double)J/(double)total-(double)(K-H)/(double)K;
delsq=delsq+temp1*temp1;
if (temp1<0.0)
temp1=-temp1;
temp=temp+temp1;
J=J+1;
}
}
count1=count1-1;
L=lagrange1(1,4,N);
ctemp2=haros5(N,0,S,1,4,L>>16,L&0xffff,8);
for (I=1; I<ctemp2; I++) {
H=S[I]>>16; // towards 3/8
K=S[I]&0xffff;
temp1=(double)(I+ctemp1+count1)/(double)total-(double)H/(double)K;
delsq=delsq+temp1*temp1;
if (temp1<0.0)
temp1=-temp1;
temp=temp+temp1;
}
if (horder==0) {
J=total-ctemp2-count1-ctemp1+1;
for (I=ctemp2-1; I>0; I--) {
H=S[I]>>16;
K=S[I]&0xffff;
temp1=(double)J/(double)total-(double)(K-H)/(double)K;
delsq=delsq+temp1*temp1;
if (temp1<0.0)
temp1=-temp1;
temp=temp+temp1;
J=J+1;
}
}
ctemp2=ctemp2-1;
L=lagrange1(3,8,N);
count2=haros5(N,0,S,3,8,L>>16,L&0xffff,2);
if (horder==0)
count2=count2-1;
for (I=1; I<count2; I++) {
H=S[I]>>16; // towards 1/2
K=S[I]&0xffff;
temp1=(double)(I+ctemp1+count1+ctemp2)/(double)total-(double)H/(double)K;
delsq=delsq+temp1*temp1;
if (temp1<0.0)
temp1=-temp1;
temp=temp+temp1;
}
if (horder==0) {
J=total-count2-ctemp2-count1-ctemp1+1;
for (I=count2-1; I>0; I--) {
H=S[I]>>16;
K=S[I]&0xffff;
temp1=(double)J/(double)total-(double)(K-H)/(double)K;
delsq=delsq+temp1*temp1;
if (temp1<0.0)
temp1=-temp1;
temp=temp+temp1;
J=J+1;
}
}
sum[2]=temp;
delsq=delsq*(double)total;
delsq=sqrt(delsq);
sum[3]=delsq;
return(2*(sum1+sum2));
}