/*CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                             C
C  COMPUTE MERTENS FUNCTION						      C
C  08/16/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 merten4d(unsigned int N, unsigned int *count, double *sum,
		unsigned int horder, unsigned int flag) {
unsigned int S[250000];
unsigned int H,K,total;
unsigned int L,I,J,count1,count2,count3,count4,ctemp1,ctemp2,ctemp3,ctemp4;
unsigned int count5,count6,count7,count8,ctemp5,ctemp6,ctemp7,ctemp8;
double temp1,sum1,sum2,temp,delsq;
double pi;
if (N==1) {
   count[0]=0;
   count[1]=0;
   if (horder==0)
      count[2]=1;
   else
      count[2]=0;
   sum[0]=0.5;
   sum[1]=0.0;
   sum[2]=0.0;
   sum[3]=0.0;
   sum[4]=0.0;
   return(1.0);
   }
if (N==2) {
   count[0]=0;
   count[1]=0;
   if (horder==0)
      count[2]=2;
   else
      count[2]=1;
   sum[0]=0.0;
   sum[1]=0.0;
   if (horder==0) {
      sum[2]=0.0;
      sum[3]=0.0;
      }
   else {
      sum[2]=0.5;
      sum[3]=0.5;
      }
   sum[4]=1.414214;
   return(0.0);
   }
if (N==3) {
   count[0]=0;
   count[1]=0;
   if (horder==0)
      count[2]=4;
   else
      count[2]=2;
   sum[0]=0.0;
   sum[1]=-0.5;
   if (horder==0) {
      sum[2]=0.1666667;
      if (flag==0)
	 sum[3]=0.2357023;
      else
	 sum[3]=0.02777778;
      }
   else {
      sum[2]=0.6666667;
      sum[3]=0.7453560;
      }
   sum[4]=2.013841;
   return(-1.0);
   }
if (N==4) {
   count[0]=0;
   count[1]=1;
   if (horder==0)
      count[2]=6;
   else
      count[2]=3;
   sum[0]=0.0;
   sum[1]=-0.5;
   if (horder==0) {
      sum[2]=0.1666667;
      if (flag==0)
	 sum[3]=0.2886751;
      else
	 sum[3]=0.03402069;
      }
   else {
      sum[2]=0.9166667;
      sum[3]=1.0507931;
      }
   sum[4]=2.466441;
   return(-1.0);
   }
//
pi=3.141592654;
if ((N>4)&&(N<33)) {
   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;
   count[2]=total;
   sum[4]=sqrt((delsq+1.0)*(double)total);
   if (flag==0) {
      delsq=delsq*(double)total;
      delsq=sqrt(delsq);
      }
   else
     delsq=sqrt((double)total)*delsq;
   sum[3]=delsq;
   return(2*(sum1+sum2));
   }
if (N>5128) {
   count[0]=0;
   count[1]=0;
   count[2]=0;
   sum[0]=0.0;
   sum[1]=0.0;
   sum[2]=0.0;
   sum[3]=0.0;
   sum[4]=0.0;
   return(0.0);
   }
ctemp1=haros5(N,0,S,0,1,1,N,32);
sum1=0.0;
for (I=1; I<ctemp1; I++) {
   H=S[I]>>16;			 // towards 1/32
   K=S[I]&0xffff;
   temp1=cos(2.0*pi*(double)H/(double)K);
   sum1=sum1+temp1;
   }
L=lagrange1(1,32,N);
count1=haros5(N,0,S,1,32,L>>16,L&0xffff,16);
for (I=1; I<count1; I++) {
   H=S[I]>>16;			 // towards 1/16
   K=S[I]&0xffff;
   temp1=cos(2.0*pi*(double)H/(double)K);
   sum1=sum1+temp1;
   }
count1=count1+ctemp1-2;
L=lagrange1(1,16,N);
ctemp1=haros5(N,0,S,1,16,L>>16,L&0xffff,32);
for (I=1; I<ctemp1; I++) {
   H=S[I]>>16;			 // towards 3/32
   K=S[I]&0xffff;
   temp1=cos(2.0*pi*(double)H/(double)K);
   sum1=sum1+temp1;
   }
count1=count1+ctemp1-1;
L=lagrange1(3,32,N);
ctemp1=haros5(N,0,S,3,32,L>>16,L&0xffff,8);
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;
   }
count1=count1+ctemp1-1;
L=lagrange1(1,8,N);
ctemp1=haros5(N,0,S,1,8,L>>16,L&0xffff,32);
for (I=1; I<ctemp1; I++) {
   H=S[I]>>16;			 // towards 5/32
   K=S[I]&0xffff;
   temp1=cos(2.0*pi*(double)H/(double)K);
   sum1=sum1+temp1;
   }
count1=count1+ctemp1-1;
L=lagrange1(5,32,N);
ctemp1=haros5(N,0,S,5,32,L>>16,L&0xffff,16);
for (I=1; I<ctemp1; I++) {
   H=S[I]>>16;			 // towards 3/16
   K=S[I]&0xffff;
   temp1=cos(2.0*pi*(double)H/(double)K);
   sum1=sum1+temp1;
   }
count1=count1+ctemp1-1;
L=lagrange1(3,16,N);
ctemp1=haros5(N,0,S,3,16,L>>16,L&0xffff,32);
for (I=1; I<ctemp1; I++) {
   H=S[I]>>16;			 // towards 7/32
   K=S[I]&0xffff;
   temp1=cos(2.0*pi*(double)H/(double)K);
   sum1=sum1+temp1;
   }
count1=count1+ctemp1-1;
L=lagrange1(7,32,N);
ctemp1=haros5(N,0,S,7,32,L>>16,L&0xffff,4);
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=count1+ctemp1-2;
count[0]=count1;
sum[0]=sum1;
//
L=lagrange1(1,4,N);
ctemp2=haros5(N,0,S,1,4,L>>16,L&0xffff,32);
sum2=0.0;
for (I=1; I<ctemp2; I++) {
   H=S[I]>>16;		       // towards 9/32
   K=S[I]&0xffff;
   temp1=cos(2.0*pi*(double)H/(double)K);
   sum2=sum2+temp1;
   }
L=lagrange1(9,32,N);
count2=haros5(N,0,S,9,32,L>>16,L&0xffff,16);
for (I=1; I<count2; I++) {
   H=S[I]>>16;		       // towards 5/16
   K=S[I]&0xffff;
   temp1=cos(2.0*pi*(double)H/(double)K);
   sum2=sum2+temp1;
   }
count2=count2+ctemp2-1;
L=lagrange1(5,16,N);
ctemp2=haros5(N,0,S,5,16,L>>16,L&0xffff,32);
for (I=1; I<ctemp2; I++) {
   H=S[I]>>16;		       // towards 11/32
   K=S[I]&0xffff;
   temp1=cos(2.0*pi*(double)H/(double)K);
   sum2=sum2+temp1;
   }
count2=count2+ctemp2-1;
L=lagrange1(11,32,N);
ctemp2=haros5(N,0,S,11,32,L>>16,L&0xffff,8);
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;
   }
count2=count2+ctemp2-1;
L=lagrange1(3,8,N);
ctemp2=haros5(N,0,S,3,8,L>>16,L&0xffff,32);
for (I=1; I<ctemp2; I++) {
   H=S[I]>>16;		       // towards 13/32
   K=S[I]&0xffff;
   temp1=cos(2.0*pi*(double)H/(double)K);
   sum2=sum2+temp1;
   }
count2=count2+ctemp2-1;
L=lagrange1(13,32,N);
ctemp2=haros5(N,0,S,13,32,L>>16,L&0xffff,16);
for (I=1; I<ctemp2; I++) {
   H=S[I]>>16;		       // towards 7/16
   K=S[I]&0xffff;
   temp1=cos(2.0*pi*(double)H/(double)K);
   sum2=sum2+temp1;
   }
count2=count2+ctemp2-1;
L=lagrange1(7,16,N);
ctemp2=haros5(N,0,S,7,16,L>>16,L&0xffff,32);
for (I=1; I<ctemp2; I++) {
   H=S[I]>>16;		       // towards 15/32
   K=S[I]&0xffff;
   temp1=cos(2.0*pi*(double)H/(double)K);
   sum2=sum2+temp1;
   }
count2=count2+ctemp2-1;
L=lagrange1(15,32,N);
ctemp2=haros5(N,0,S,15,32,L>>16,L&0xffff,2);
ctemp2=ctemp2-1;
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=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,32);
for (I=1; I<ctemp1; I++) {
   H=S[I]>>16;			 // towards 1/32
   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,32,N);
count1=haros5(N,0,S,1,32,L>>16,L&0xffff,16);
for (I=1; I<count1; I++) {
   H=S[I]>>16;			 // towards 1/16
   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,16,N);
ctemp2=haros5(N,0,S,1,16,L>>16,L&0xffff,32);
for (I=1; I<ctemp2; I++) {
   H=S[I]>>16;		       // towards 3/32
   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,32,N);
count2=haros5(N,0,S,3,32,L>>16,L&0xffff,8);
for (I=1; I<count2; I++) {
   H=S[I]>>16;		       // towards 1/8
   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;
      }
   }
count2=count2-1;
L=lagrange1(1,8,N);
ctemp3=haros5(N,0,S,1,8,L>>16,L&0xffff,32);
for (I=1; I<ctemp3; I++) {
   H=S[I]>>16;		       // towards 5/32
   K=S[I]&0xffff;
   temp1=(double)(I+ctemp1+count1+ctemp2+count2)/(double)total-(double)H/(double)K;
   delsq=delsq+temp1*temp1;
   if (temp1<0.0)
      temp1=-temp1;
   temp=temp+temp1;
   }
if (horder==0) {
   J=total-ctemp3-count2-ctemp2-count1-ctemp1+1;
   for (I=ctemp3-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;
      }
   }
ctemp3=ctemp3-1;
L=lagrange1(5,32,N);
count3=haros5(N,0,S,5,32,L>>16,L&0xffff,16);
for (I=1; I<count3; I++) {
   H=S[I]>>16;		       // towards 3/16
   K=S[I]&0xffff;
   temp1=(double)(I+ctemp1+count1+ctemp2+count2+ctemp3)/(double)total-(double)H/(double)K;
   delsq=delsq+temp1*temp1;
   if (temp1<0.0)
      temp1=-temp1;
   temp=temp+temp1;
   }
if (horder==0) {
   J=total-count3-ctemp3-count2-ctemp2-count1-ctemp1+1;
   for (I=count3-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;
      }
   }
count3=count3-1;
L=lagrange1(3,16,N);
ctemp4=haros5(N,0,S,3,16,L>>16,L&0xffff,32);
for (I=1; I<ctemp4; I++) {
   H=S[I]>>16;		       // towards 7/32
   K=S[I]&0xffff;
   temp1=(double)(I+ctemp1+count1+ctemp2+count2+ctemp3+count3)/(double)total-(double)H/(double)K;
   delsq=delsq+temp1*temp1;
   if (temp1<0.0)
      temp1=-temp1;
   temp=temp+temp1;
   }
if (horder==0) {
   J=total-ctemp4-count3-ctemp3-count2-ctemp2-count1-ctemp1+1;
   for (I=ctemp4-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;
      }
   }
ctemp4=ctemp4-1;
L=lagrange1(7,32,N);
count4=haros5(N,0,S,7,32,L>>16,L&0xffff,4);
for (I=1; I<count4; I++) {
   H=S[I]>>16;		       // towards 1/4
   K=S[I]&0xffff;
   temp1=(double)(I+ctemp1+count1+ctemp2+count2+ctemp3+count3+ctemp4)/(double)total-(double)H/(double)K;
   delsq=delsq+temp1*temp1;
   if (temp1<0.0)
      temp1=-temp1;
   temp=temp+temp1;
   }
if (horder==0) {
   J=total-count4-ctemp4-count3-ctemp3-count2-ctemp2-count1-ctemp1+1;
   for (I=count4-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;
      }
   }
count4=count4-1;
L=lagrange1(1,4,N);
ctemp5=haros5(N,0,S,1,4,L>>16,L&0xffff,32);
for (I=1; I<ctemp5; I++) {
   H=S[I]>>16;		       // towards 9/32
   K=S[I]&0xffff;
   temp1=(double)(I+ctemp1+count1+ctemp2+count2+ctemp3+count3+ctemp4+count4)/(double)total-(double)H/(double)K;
   delsq=delsq+temp1*temp1;
   if (temp1<0.0)
      temp1=-temp1;
   temp=temp+temp1;
   }
if (horder==0) {
   J=total-ctemp5-count4-ctemp4-count3-ctemp3-count2-ctemp2-count1-ctemp1+1;
   for (I=ctemp5-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;
      }
   }
ctemp4=ctemp4-1;
L=lagrange1(9,32,N);
count5=haros5(N,0,S,9,32,L>>16,L&0xffff,16);
for (I=1; I<count5; I++) {
   H=S[I]>>16;		       // towards 5/16
   K=S[I]&0xffff;
   temp1=(double)(I+ctemp1+count1+ctemp2+count2+ctemp3+count3+ctemp4+count4+ctemp5)/(double)total-(double)H/(double)K;
   delsq=delsq+temp1*temp1;
   if (temp1<0.0)
      temp1=-temp1;
   temp=temp+temp1;
   }
if (horder==0) {
   J=total-count5-ctemp5-count4-ctemp4-count3-ctemp3-count2-ctemp2-count1-ctemp1+1;
   for (I=count5-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;
      }
   }
count5=count5-1;
L=lagrange1(5,16,N);
ctemp6=haros5(N,0,S,5,16,L>>16,L&0xffff,32);
for (I=1; I<ctemp6; I++) {
   H=S[I]>>16;		       // towards 11/32
   K=S[I]&0xffff;
   temp1=(double)(I+ctemp1+count1+ctemp2+count2+ctemp3+count3+ctemp4+count4+ctemp5+count5)/(double)total-(double)H/(double)K;
   delsq=delsq+temp1*temp1;
   if (temp1<0.0)
      temp1=-temp1;
   temp=temp+temp1;
   }
if (horder==0) {
   J=total-ctemp6-count5-ctemp5-count4-ctemp4-count3-ctemp3-count2-ctemp2-count1-ctemp1+1;
   for (I=ctemp6-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;
      }
   }
ctemp6=ctemp6-1;
L=lagrange1(11,32,N);
count6=haros5(N,0,S,11,32,L>>16,L&0xffff,8);
for (I=1; I<count6; I++) {
   H=S[I]>>16;		       // towards 3/8
   K=S[I]&0xffff;
   temp1=(double)(I+ctemp1+count1+ctemp2+count2+ctemp3+count3+ctemp4+count4+ctemp5+count5+ctemp6)/(double)total-(double)H/(double)K;
   delsq=delsq+temp1*temp1;
   if (temp1<0.0)
      temp1=-temp1;
   temp=temp+temp1;
   }
if (horder==0) {
   J=total-count6-ctemp6-count5-ctemp5-count4-ctemp4-count3-ctemp3-count2-ctemp2-count1-ctemp1+1;
   for (I=count6-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;
      }
   }
count6=count6-1;
L=lagrange1(3,8,N);
ctemp7=haros5(N,0,S,3,8,L>>16,L&0xffff,32);
for (I=1; I<ctemp7; I++) {
   H=S[I]>>16;		       // towards 13/32
   K=S[I]&0xffff;
   temp1=(double)(I+ctemp1+count1+ctemp2+count2+ctemp3+count3+ctemp4+count4+ctemp5+count5+ctemp6+count6)/(double)total-(double)H/(double)K;
   delsq=delsq+temp1*temp1;
   if (temp1<0.0)
      temp1=-temp1;
   temp=temp+temp1;
   }
if (horder==0) {
   J=total-ctemp7-count6-ctemp6-count5-ctemp5-count4-ctemp4-count3-ctemp3-count2-ctemp2-count1-ctemp1+1;
   for (I=ctemp7-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;
      }
   }
ctemp7=ctemp7-1;
L=lagrange1(13,32,N);
count7=haros5(N,0,S,13,32,L>>16,L&0xffff,16);
for (I=1; I<count7; I++) {
   H=S[I]>>16;		       // towards 7/16
   K=S[I]&0xffff;
   temp1=(double)(I+ctemp1+count1+ctemp2+count2+ctemp3+count3+ctemp4+count4+ctemp5+count5+ctemp6+count6+ctemp7)/(double)total-(double)H/(double)K;
   delsq=delsq+temp1*temp1;
   if (temp1<0.0)
      temp1=-temp1;
   temp=temp+temp1;
   }
if (horder==0) {
   J=total-count7-ctemp7-count6-ctemp6-count5-ctemp5-count4-ctemp4-count3-ctemp3-count2-ctemp2-count1-ctemp1+1;
   for (I=count7-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;
      }
   }
count7=count7-1;
L=lagrange1(7,16,N);
ctemp8=haros5(N,0,S,7,16,L>>16,L&0xffff,32);
for (I=1; I<ctemp8; I++) {
   H=S[I]>>16;		       // towards 15/32
   K=S[I]&0xffff;
   temp1=(double)(I+ctemp1+count1+ctemp2+count2+ctemp3+count3+ctemp4+count4+ctemp5+count5+ctemp6+count6+ctemp7+count7)/(double)total-(double)H/(double)K;
   delsq=delsq+temp1*temp1;
   if (temp1<0.0)
      temp1=-temp1;
   temp=temp+temp1;
   }
if (horder==0) {
   J=total-ctemp8-count7-ctemp7-count6-ctemp6-count5-ctemp5-count4-ctemp4-count3-ctemp3-count2-ctemp2-count1-ctemp1+1;
   for (I=ctemp8-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;
      }
   }
ctemp8=ctemp8-1;
L=lagrange1(15,32,N);
count8=haros5(N,0,S,15,32,L>>16,L&0xffff,2);
if (horder==0)
   count8=count8-1;
for (I=1; I<count8; I++) {
   H=S[I]>>16;		       // towards 1/2
   K=S[I]&0xffff;
   temp1=(double)(I+ctemp1+count1+ctemp2+count2+ctemp3+count3+ctemp4+count4+ctemp5+count5+ctemp6+count6+ctemp7+count7+ctemp8)/(double)total-(double)H/(double)K;
   delsq=delsq+temp1*temp1;
   if (temp1<0.0)
      temp1=-temp1;
   temp=temp+temp1;
   }
if (horder==0) {
   J=total-count8-ctemp8-count7-ctemp7-count6-ctemp6-count5-ctemp5-count4-ctemp4-count3-ctemp3-count2-ctemp2-count1-ctemp1+1;
   for (I=count8-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;
count[2]=total;
sum[4]=sqrt((delsq+1.0)*(double)total);
if (flag==0) {
   delsq=delsq*(double)total;
   delsq=sqrt(delsq);
   }
else
   delsq=sqrt((double)total)*delsq;
sum[3]=delsq;
return(2*(sum1+sum2));
}