/*CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                             C
C  COMPUTE DIFFERENCES IN NUMBER OF FRACTIONS (0 to 1/I, 1/I to 2/I)	      C
C  08/25/15 (DKC)							      C
C                                                                             C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC*/
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include "table2.h"
int liouvill(unsigned int a, unsigned int *t);
double ramanuj(unsigned int n, unsigned int k, unsigned int *t);
double numdiv(unsigned int a, unsigned int flag);
double mangoldt(unsigned int a, unsigned int *t);
void gauss(unsigned int n, unsigned int k, double *x, double *o);
mertenaw(unsigned int, unsigned int *count, unsigned int I);
//
unsigned int limit=0;	// partial sum flag, normally set to 0
unsigned int norm=0;	// normalize, normally set to 0
unsigned int hflag=0;   // histogram flag
unsigned int MAXN=10000; // 2000000000 maximum
unsigned int BEGINN=2;
unsigned int delta=1;	// subtract slope if set
unsigned int I=100;	// n value (for shorter intervals of Farey points)
unsigned int k=1373;	// for out=38-45, 55 (Ramanujan sum)
unsigned int reim=0;	// for out=46-54, 56,57 output real or imaginary part
unsigned int power=0;	// decimate by J (normally set to 0)
unsigned int scale=1;	// for out=37
unsigned int mflag=0;	// 0,1,2,3, or 4 additional processing
			//   for out=13,37,38,50,61,62 (Fourier coefficients, etc.)
unsigned int again=0;	// for out=37 (additional processing)
unsigned int prime=1;	// for out=41,58 (k should be prime if set)
unsigned int switchi=0; // for out=46-54, 56,57 select Dirichlet characters
unsigned int out=62; // select output
		     // 10, sgn(n(x)-m(x))+1.781072418*sum(sgn(n(x/i)-m(x/i))*i*log(log(i)))
		     // 11, sum of m(x/i)
		     // 12, sum of n(x/i)
		     // 13, sum of n(x/i)-m(x/i)
		     // 14, sum of (n(x/i)-m(x/i))*i
		     // 15, |n(x)-m(x)|+1.781072418*sum(|n(x/i)-m(x/i)|*i*log(log(i)))
		     // 16, sum of (n(x/i)-m(x/i))*(h+exp(h)*log(h))
		     // 17, differences in number of fractions before, after 1/I
		     // 18, number of fractions before 1/I
		     // 19, number of fractions after 1/I
		     // 20, n(x)-m(x)+1.781072418*sum((n(x/i)-m(x/i))*i*log(log(i)))
		     // 21, sum of (n(x/i)-m(x/i))*2*d(i) (see #24)
		     // 22, sum of (n(x/i)-m(x/i))*sigma(i)
		     // 23, sum of (n(x/i)-m(x/i)^2*L(i)
		     // 24, sum of (n(x/i)-m(x/i))*d(i) (see #21)
		     // 25, sum of n(x/i)-m(x/i) when i is a perfect square
		     // 26, sum of (n(x/i)-m(x/i))^2
		     // 27, sum of (n(x/i)-m(x/i))*log(i)
		     // 28, sum of sgn(n(x/i)-m(x/i))
		     // 29, sum of (n(x/i)-m(x/i)*log(i)*d(i)
		     // 30, sum of (n(x/i)-m(x/i)^2*L(i)^2
		     // 31, sum of sgn(n(x/i)-m(x/i))*i
		     // 32, sum of |sgn(n(x/i)-m(x/i))|
		     // 33, sum of (n(x/i)-m(x/i))*L(i)
		     // 34, sum of (n(x/i)-m(x/i))*M(i)
		     // 35, --
		     // 36, sum of (n(x/i)-m(x/i))*sigma2(i)
		     // 37, sum of (n(x/i)-m(x/i))*mangoldt(i)
		     // 38, sum of (n(x/i)-m(x/i))*ramanuj(i)
		     // 39, sum of M(x/i)*ramanuj(i)
		     // 40, sum of ramanuj(x/i)*sigma(i)
		     // 41, sum of ramanuj(x/i)*M(i)
		     // 42, sum of ramanuj(x/i)*log(i)
		     // 43, sum of ramanuj(x/i)
		     // 44, sum of ramanuj(x/i)*2*d(i)
		     // 45, sum of ramanuj(x/i)*sigma2(i)
		     // 46, sum of gauss(x/i)*sigma2(i)
		     // 47, sum of gauss(x/i)
		     // 48, sum of gauss(x/i)*log(i)
		     // 49, sum of gauss(x/i)*sigma(i)
		     // 50, sum of (n(x/i)-m(x/i))*gauss(i)
		     // 51, sum of gauss(x/i)*2*d(i)
		     // 52, sum of gauss(x/i)*log(i)*d(i)
		     // 53, sum of gauss(x/i)*mangoldt(i)
		     // 54, sum of M(x/i)*gauss(i)
		     // 55, sum of ramanuj(x/i)(n(x)-m(x))
		     // 56, sum of gauss(x/i)*(n(x)-m(x))
		     // 57, gauss(x)
		     // 58, sum of ramanuj(x/i)*M(i)*M(i)
		     // 59, ramanuj(x)
		     // 60, sum of M(i)^2 for k divides N/i
		     // 61, equivalent to 60
		     // 62, sum of M(x/i)^2 where i divides x
		     // 63, sum of M(x/i)^2
//
void main() {
int t,tsum,sum;
unsigned int N,i,temp,count[2],iters,J,JP,first,kp;
double sum6,sum7,sum8,sum9,sum10,sum11,sum12,sum13,sum14,sum15,sum16,sum17,sum18;
double sum1,sum2,sum3,sum4,sum5,h,temp1,resum,imsum;
unsigned int histo[301],maxt,saveN[4],savet[4];
double pi,mean,std,msum,msum2,dsum,d2sum,*output,*outman,*manout;
double gre[10000],gim[10000],xi[14],xip[22],xiq[26],outp[2];
double xi1[14],xi2[14],xi3[14],xi4[14],xi5[14],xi6[14];
//short M[60001];
short *M;
FILE *Outfp;
Outfp = fopen("out1bz.dat","w");
pi=3.141592654;
if (I<4) {
   printf("I too small \n");
   goto zskip;
   }
if ((out==61)&&(MAXN<k)) {
   printf("error 61: \n");
   goto zskip;
   }
if ((out==37)||(out==38)||(out==13)||(out==50)) {
   output=(double*) malloc(1000001);
   if (output==NULL)
      return;
   outman=(double*) malloc(1000001);
   if (outman==NULL)
      return;
   manout=(double*) malloc(1000001);
   if (outman==NULL)
      return;
   outman[0]=0.0;
   }
J=(I+1)/2;
JP=1;
if (power!=0)
   JP=J;
//
// compute Mertens function
//
if ((out==34)||(out==39)||(out==41)||(out==54)||(out==58)||(out==60)||(out==61)||(out==62)||(out==63)) {
   M=(short*) malloc(2000001);
   if (M==NULL)
      return;
   if (MAXN>1000000) {
      printf("not enough memory allocated");
      goto zskip;
      }
   M[0]=1;
   for (N=2; N<=(MAXN+1); N++) {
      sum=0;
      for (i=2; i<=(N/3); i++)
	 sum=sum+M[N/i-1];
      sum=sum+(N+1)/2;
      t=1-sum;
      if ((t>32767)||(t<-32768)) {
	 printf("overflow \n");
	 goto zskip;
	 }
      M[N-1]=(short)t;
//	  printf(" %d \n",M[N-1]);
      }
   }
for (i=0; i<301; i++)
   histo[i]=0;
//
// Dirichlet characters mod 7
//
xi[0]=1.0;
xi[1]=0.0;
xi[2]=cos(2.0*pi/3.0);
xi[3]=sin(2.0*pi/3.0);
xi[4]=cos(pi/3.0);
xi[5]=sin(pi/3.0);
xi[6]=-xi[4];
xi[7]=-xi[5];
xi[8]=-xi[2];
xi[9]=-xi[3];
xi[10]=-1.0;
xi[11]=0.0;
xi[12]=0.0;
xi[13]=0.0;
//
xi1[0]=1.0;
xi1[1]=0.0;
xi1[2]=1.0;
xi1[3]=0.0;
xi1[4]=1.0;
xi1[5]=0.0;
xi1[6]=1.0;
xi1[7]=0.0;
xi1[8]=1.0;
xi1[9]=0.0;
xi1[10]=1.0;
xi1[11]=0.0;
xi1[12]=0.0;
xi1[13]=0.0;
//
xi2[0]=1.0;
xi2[1]=0.0;
xi2[2]=1.0;
xi2[3]=0.0;
xi2[4]=-1.0;
xi2[5]=0.0;
xi2[6]=1.0;
xi2[7]=0.0;
xi2[8]=-1.0;
xi2[9]=0.0;
xi2[10]=-1.0;
xi2[11]=0.0;
xi2[12]=0.0;
xi2[13]=0.0;
//
xi3[0]=1.0;
xi3[1]=0.0;
xi3[2]=cos(2.0*pi/3.0);
xi3[3]=sin(2.0*pi/3.0);
xi3[4]=cos(pi/3.0);
xi3[5]=sin(pi/3.0);
xi3[6]=-xi[4];
xi3[7]=-xi[5];
xi3[8]=-xi[2];
xi3[9]=-xi[3];
xi3[10]=-1.0;
xi3[11]=0.0;
xi3[12]=0.0;
xi3[13]=0.0;
//
xi4[0]=1.0;
xi4[1]=0.0;
xi4[2]=cos(2.0*pi/3.0);
xi4[3]=sin(2.0*pi/3.0);
xi4[4]=-cos(pi/3.0);
xi4[5]=-sin(pi/3.0);
xi4[6]=xi[4];
xi4[7]=xi[5];
xi4[8]=xi[2];
xi4[9]=xi[3];
xi4[10]=1.0;
xi4[11]=0.0;
xi4[12]=0.0;
xi4[13]=0.0;
//
xi5[0]=1.0;
xi5[1]=0.0;
xi5[2]=-cos(pi/3.0);
xi5[3]=-sin(pi/3.0);
xi5[4]=cos(2.0*pi/3.0);
xi5[5]=sin(2.0*pi/3.0);
xi5[6]=xi[4];
xi5[7]=xi[5];
xi5[8]=xi[2];
xi5[9]=xi[3];
xi5[10]=1.0;
xi5[11]=0.0;
xi5[12]=0.0;
xi5[13]=0.0;
//
xi6[0]=1.0;
xi6[1]=0.0;
xi6[2]=-cos(pi/3.0);
xi6[3]=-sin(pi/3.0);
xi6[4]=-cos(2.0*pi/3.0);
xi6[5]=-sin(2.0*pi/3.0);
xi6[6]=-xi[4];
xi6[7]=-xi[5];
xi6[8]=-xi[2];
xi6[9]=-xi[3];
xi6[10]=-1.0;
xi6[11]=0.0;
xi6[12]=0.0;
xi6[13]=0.0;
//
//for (i=0; i<7; i++)
// printf(" %e %e \n",xi[2*i],xi[2*i+1]);
//
// Dirichlet character mod 11
//
xip[0]=1.0;
xip[1]=0.0;
xip[2]=cos(pi/5.0);
xip[3]=sin(pi/5.0);
xip[4]=-cos(3.0*pi/5.0);
xip[5]=-sin(3.0*pi/5.0);
xip[6]=cos(2.0*pi/5.0);
xip[7]=sin(2.0*pi/5.0);
xip[8]=cos(4.0*pi/5.0);
xip[9]=sin(4.0*pi/5.0);
xip[10]=-xip[8];
xip[11]=-xip[9];
xip[12]=-xip[6];
xip[13]=-xip[7];
xip[14]=-xip[4];
xip[15]=-xip[5];
xip[16]=-xip[2];
xip[17]=-xip[3];
xip[18]=-1.0;
xip[19]=0.0;
xip[20]=0.0;
xip[21]=0.0;
//
// Dirichlet character mod 13
//
xiq[0]=1.0;
xiq[1]=0.0;
xiq[2]=-cos(pi/6.0);
xiq[3]=-sin(pi/6.0);
xiq[4]=cos(2.0*pi/3.0);
xiq[5]=sin(2.0*pi/3.0);
xiq[6]=cos(pi/3.0);
xiq[7]=sin(pi/3.0);
xiq[8]=0.0;
xiq[9]=1.0;
xiq[10]=-cos(5.0*pi/6.0);
xiq[11]=-sin(5.0*pi/6.0);
xiq[12]=-xiq[10];
xiq[13]=-xiq[11];
xiq[14]=-xiq[8];
xiq[15]=-xiq[9];
xiq[16]=-xiq[6];
xiq[17]=-xiq[7];
xiq[18]=-xiq[4];
xiq[19]=-xiq[5];
xiq[20]=-xiq[2];
xiq[21]=-xiq[3];
xiq[22]=-1.0;
xiq[23]=0.0;
xiq[24]=0.0;
xiq[25]=0.0;
//
iters=0;
dsum=0.0;
d2sum=0.0;
first=0;
maxt=0;
saveN[0]=0;
saveN[1]=0;
saveN[2]=0;
saveN[3]=0;
for (N=BEGINN; N<=MAXN; N++) {
   if (N!=(N/JP)*JP)
      continue;
   if ((out!=60)&&(out!=61)&&(out!=62))
      printf(" %d \n",N);
   if ((out==17)||(out==18)||(out==19))
      mertenaw(N, count, I);
   if (out==17) {
      t=(int)count[0]-(int)count[1];
      if (hflag==2) {
	 dsum=dsum+(double)t;
	 d2sum=d2sum+(double)t*(double)t;
	 iters=iters+1;
	 std=d2sum-dsum*dsum/(double)iters;
	 if (iters!=1)
	    std=sqrt(std/(double)(iters-1));
	 else
	    std=0.0;
	 fprintf(Outfp," %e \n",std);
	 }
      if (hflag==0)
	 fprintf(Outfp," %d\n",t);
      if (hflag==1) {
	 if ((t<-150)||(t>150)) {
	    printf("difference out of range \n");
	    goto zskip;
	    }
	 histo[t+150]=histo[t+150]+1;
	 }
      continue;
      }
   if (out==18) {
      fprintf(Outfp," %d\n",count[0]);
      continue;
      }
   if (out==19) {
      fprintf(Outfp," %d\n",count[1]);
      continue;
      }
   if ((out==11)||(out==12))  {
      sum6=0.0;
      sum7=0.0;
      for (i=1; i<=N; i++) {
	 mertenaw(N/i, count, I);
	 sum6=sum6+(double)count[0];
	 sum7=sum7+(double)count[1];
	 }
      if (out==11)
	 fprintf(Outfp," %e\n",sum6);
      if (out==12)
	 fprintf(Outfp," %e\n",sum7);
      continue;
      }
   if ((out==13)||(out==14)) {
      sum8=0.0;
      sum9=0.0;
      for (i=1; i<=N; i++) {
	 mertenaw(N/i, count, I);
	 if (delta==0)
	    sum8=sum8+(double)count[1]-(double)count[0];	      // h(x)
	 else
	    sum8=sum8+(double)count[0]-(double)count[1]+(double)(J-1)/(double)I; // h(x)
	 sum9=sum9+(double)i*((double)count[1]-(double)count[0]);
	 }
      if (out==13) {
	 if (norm==0) {
	    if (mflag==0)
	       fprintf(Outfp," %e\n",sum8);
	    else
	       outman[N-1]=sum8;
	    }
	 else
	    fprintf(Outfp," %e\n",sum8/(double)N);
	 }
      if (out==14)
	 fprintf(Outfp," %e\n",sum9);
      continue;
      }
   if (out==10) {
      sum2=0.0;
      for (i=2; i<=N; i++) {
	 mertenaw(N/i, count, I);
	 temp1=(double)count[1]-(double)count[0];
	 if (temp1<0.0)
	    sum2=sum2-(double)i*log(log((double)i));
	 if (temp1>0.0)
	    sum2=sum2+(double)i*log(log((double)i));
	 }
      sum2=1.781072418*sum2;
      mertenaw(N, count, I);
      temp1=(double)count[1]-(double)count[0];
      if (temp1<0.0)
	 sum2=sum2-1;
      if (temp1>0.0)
	 sum2=sum2+1;
      fprintf(Outfp," %e\n",sum2);
      continue;
      }
   if (out==15) {
      sum2=0.0;
      for (i=2; i<=N; i++) {
	 mertenaw(N/i, count, I);
	 temp1=(double)count[1]-(double)count[0];
	 if (temp1<0.0)
	    temp1=-temp1;
	 sum2=sum2+temp1*(double)i*log(log((double)i));
	 }
      mertenaw(N, count, I);
      temp1=(double)count[1]-(double)count[0];
      if (temp1<0.0)
	 temp1=-temp1;
      sum2=1.781072418*sum2+temp1;
      fprintf(Outfp," %e\n",sum2);
      continue;
      }
   if (out==16) {
      sum1=0.0;
      h=0.0;
      for (i=1; i<=N; i++) {
	 mertenaw(N/i, count, I);
	 h=h+1.0/(double)i;
	 sum1=sum1+((double)count[1]-(double)count[0])*(h+exp(h)*log(h));
	 }
      fprintf(Outfp," %e\n",sum1);
      continue;
      }
   if (out==20) {
      sum2=0.0;
      for (i=2; i<=N; i++) {
	 mertenaw(N/i, count, I);
	 sum2=sum2+((double)count[1]-(double)count[0])*(double)i*log(log((double)i));
	 }
      mertenaw(N, count, I);
      sum2=1.781072418*sum2+(double)count[1]-(double)count[0];
      fprintf(Outfp," %e\n",sum2);
      continue;
      }
   if (out==21) {
      sum3=0.0;
      for (i=1; i<=N; i++) {
	 mertenaw(N/i, count, I);
	 if (delta==0)
	    sum3=sum3+((double)count[1]-(double)count[0])*2.0*numdiv(i,1);
	 else
	    sum3=sum3+((double)count[0]-(double)count[1]+(double)(J-1)/(double)I)*2.0*numdiv(i,1);
	 }
      fprintf(Outfp," %e\n",sum3);
      continue;
      }
   if (out==22) {
      sum4=0.0;
      for (i=1; i<=N; i++) {
	 mertenaw(N/i, count, I);
	 if (delta==0)
	    sum4=sum4+((double)count[1]-(double)count[0])*numdiv(i,4);
	 else
	    sum4=sum4+((double)count[0]-(double)count[1]+(double)(J-1)/(double)I)*numdiv(i,4);
	 }
      fprintf(Outfp," %e\n",sum4);
      continue;
      }
   if (out==23) {
      sum5=0.0;
      tsum=0;
      for (i=1; i<=N; i++) {
	 mertenaw(N/i, count, I);
	 t=liouvill(i, table);
	 tsum=tsum+t;
	 sum5=sum5+((double)count[1]-(double)count[0])*((double)count[1]-(double)count[0])*(double)tsum;
	 }
      if (norm!=0)
	 sum5=sum5/(double)N;
      fprintf(Outfp," %e\n",sum5);
      continue;
      }
   if (out==24) {
      sum14=0.0;
      for (i=1; i<=N; i++) {
	 mertenaw(N/i, count, I);
	 if (delta==0)
	    sum14=sum14+((double)count[1]-(double)count[0])*numdiv(i,1);
	 else
	    sum14=sum14+((double)count[0]-(double)count[1]+(double)(J-1)/(double)I)*numdiv(i,1);
	 }
      fprintf(Outfp," %e\n",sum14);
      continue;
      }
   if (out==25) {
      sum13=0.0;
      for (i=1; i<=N; i++) {
	 temp=(unsigned int)(sqrt((double)i)+0.001);
	 if ((temp*temp)==i) {
	    mertenaw(N/i, count, I);
	    if (delta==0)
	       sum13=sum13+(double)count[1]-(double)count[0];
	    else
	       sum13=sum13+(double)count[0]-(double)count[1]+(double)(J-1)/(double)I;
	    }
	 }
      fprintf(Outfp," %e\n",sum13);
      continue;
      }
   if (out==26) {
      sum10=0.0;
      for (i=1; i<=N; i++) {
	 if ((N/i)<=limit)
	    break;
	 mertenaw(N/i, count, I);
	 if (delta==0)
	    sum10=sum10+((double)count[1]-(double)count[0])*((double)count[1]-(double)count[0]);
	 else
	    sum10=sum10+((double)count[0]-(double)count[1]+(double)(J-1)/(double)I)*((double)count[0]-(double)count[1]+(double)(J-1)/(double)I);
	 }
      fprintf(Outfp," %e\n",sum10);
      continue;
      }
   if (out==27) {
      sum11=0.0;
      for (i=1; i<=N; i++) {
	 mertenaw(N/i, count, I);
	 if (delta==0)
	    sum11=sum11+((double)count[1]-(double)count[0])*log((double)i);
	 else
	    sum11=sum11+((double)count[0]-(double)count[1]+(double)(J-1)/(double)I)*log((double)i);
	 }
      fprintf(Outfp," %e\n",sum11);
      continue;
      }
   if (out==28) {
      sum12=0.0;
      for (i=1; i<=N; i++) {
	 mertenaw(N/i, count, I);
	 if (count[0]>count[1])       // sum of sgn(n(x/i)-m(x/i))
	    sum12=sum12-1.0;
	 if (count[0]<count[1])
	    sum12=sum12+1.0;
	 }
      if (norm==0)
	 fprintf(Outfp," %e\n",sum12);
      else
	 fprintf(Outfp," %e\n",sum12/(double)N);
      continue;
      }
   if (out==29) {
      sum18=0.0;
      for (i=1; i<=N; i++) {
	 mertenaw(N/i, count, I);
	 if (delta==0)
	    sum18=sum18+((double)count[1]-(double)count[0])*log((double)i)*numdiv(i,1);
	 else
	    sum18=sum18+((double)count[0]-(double)count[1]+(double)(J-1)/(double)I)*log((double)i)*numdiv(i,1);
	 }
      fprintf(Outfp," %e\n",sum18);
      continue;
      }
   if (out==30) {
      sum17=0.0;
      tsum=0;
      for (i=1; i<=N; i++) {
	 mertenaw(N/i, count, I);
	 t=liouvill(i, table);
	 tsum=tsum+t;
	 sum17=sum17+((double)count[1]-(double)count[0])*((double)count[1]-(double)count[0])*(double)tsum*(double)tsum;
	 }
      if (norm!=0)
	 sum17=sum17/(double)N;
      fprintf(Outfp," %e\n",sum17);
      continue;
      }
   if (out==31) {
      sum15=0.0;
      for (i=1; i<=N; i++) {
	 mertenaw(N/i, count, I);
	 if (count[0]>count[1])        // sum of sgn(n(x/i)-m(x/i))*i
	    sum15=sum15-(double)i;
	 if (count[0]<count[1])
	    sum15=sum15+(double)i;
	 }
      fprintf(Outfp," %e\n",sum15);
      continue;
      }
   if (out==32) {
      sum16=0.0;
      for (i=1; i<=N; i++) {
	 mertenaw(N/i, count, I);
	 if (count[0]!=count[1])      // sum of |sgn(n(x/i)-m(x/i))|
	    sum16=sum16+1.0;
	 }
      if (norm==0)
	 fprintf(Outfp," %e\n",sum16);
      else
	 fprintf(Outfp," %e\n",sum16/(double)N);
      }
   if (out==33) {
      sum17=0.0;
      tsum=0;
      for (i=1; i<=N; i++) {
	 mertenaw(N/i, count, I);
	 t=liouvill(i, table);
	 tsum=tsum+t;
	 sum17=sum17+((double)count[1]-(double)count[0])*(double)tsum;
	 }
      if (norm!=0)
	 sum17=sum17/(double)N;
      fprintf(Outfp," %e\n",sum17);
      continue;
      }
   if (out==34) {
      sum18=0.0;
      for (i=1; i<=N; i++) {
	 mertenaw(N/i, count, I);
	 if (delta==0)
	    sum18=sum18+((double)count[1]-(double)count[0])*(double)M[i-1];
	 else {
	    if (delta==1)
	       sum18=sum18+((double)count[0]-(double)count[1]+(double)(J-1)/(double)I)*(double)M[i-1];
	    else {
	       if (delta==2)
		  sum18=sum18+(double)count[0]*(double)M[i-1];
	       else {
		  if (delta==3)
		     sum18=sum18+(double)count[1]*(double)M[i-1];
		  else {
		     if (delta==4)
			sum18=sum18+((double)count[1]-(double)count[0])*(double)M[N/i-1];
		     else
			sum18=sum18+((double)count[0]-(double)count[1]+(double)(J-1)/(double)I)*(double)M[N/i-1];
		     }
		  }
	       }
	    }
	 }
      fprintf(Outfp," %e\n",sum18);
      continue;
      }
   if (out==36) {
      sum4=0.0;
      for (i=1; i<=N; i++) {
	 mertenaw(N/i, count, I);
	 if (delta==0)
	    sum4=sum4+((double)count[1]-(double)count[0])*numdiv(i,5);
	 else
	    sum4=sum4+((double)count[0]-(double)count[1]+(double)(J-1)/(double)I)*numdiv(i,5);
	 }
      fprintf(Outfp," %e\n",sum4);
      continue;
      }
   if (out==37) {
      sum4=0.0;
      for (i=1; i<=N; i++) {
	 mertenaw(N/i, count, I);
	 if (delta==0) {
	    if (scale==0)
	       sum4=sum4+((double)count[1]-(double)count[0])*mangoldt(i,table);
	    else
	       sum4=sum4+(((double)count[0]-(double)count[1])/(double)I)*mangoldt(i,table);
	    }
	 else
	    sum4=sum4+((double)count[0]-(double)count[1]+(double)(J-1)/(double)I)*mangoldt(i,table);
	 }
      if (mflag==0)
	 fprintf(Outfp," %e\n",sum4);
      else {
	 outman[N-1]=sum4;
	 mertenaw(N, count, I);
	 output[N-1]=((double)count[0]-(double)count[1])/(double)I;
	 }
      continue;
      }
   if (out==38) {
      sum4=0.0;
      for (i=1; i<=N; i++) {
	 mertenaw(N/i, count, I);
	 if (delta==0)
	    sum4=sum4+((double)count[1]-(double)count[0])*ramanuj(i,k,table);
	 else
	    sum4=sum4+((double)count[0]-(double)count[1]+(double)(J-1)/(double)I)*ramanuj(i,k,table);
	 }
      if (mflag==0)
	 fprintf(Outfp," %e\n",sum4);
      outman[N-1]=sum4;
      if (N>(I*k)) {
	 if (delta==0) {
	    if (sum4!=outman[N-I*k-1]) {
	       printf("error: %e %e \n",sum4,outman[N-I*k-1]);
	       goto zskip;
	       }
	    }
	 else {
	    temp1=sum4-outman[N-I*k-1];
	    if (temp1<0.0)
	       temp1=-temp1;
	    if (temp1>0.0000001) {
	       printf("warning: %e %e \n",sum4,outman[N-I*k-1]);
	       if (first==0)
		  first=1;
	       else {
		  printf("error: %e %e \n",sum4,outman[N-I*k-1]);
		  goto zskip;
		  }
	       }
	    }
	 }
      continue;
      }
   if (out==39) {
      sum4=0.0;
      for (i=1; i<=N; i++)
	 sum4=sum4+(double)M[N/i-1]*ramanuj(i,k,table);
      fprintf(Outfp," %e\n",sum4);
      continue;
      }
   if (out==40) {
      sum4=0.0;
      for (i=1; i<=N; i++)
	 if (delta==0)
	    sum4=sum4+ramanuj(N/i,k,table)*numdiv(i,4);
	 else
	    sum4=sum4+(ramanuj(N/i,k,table)+1.0)*numdiv(i,4);
      fprintf(Outfp," %e\n",sum4);
      continue;
      }
   if (out==41) {
      sum4=0.0;
      for (i=1; i<=N; i++) {
	 if (delta==0)
	    sum4=sum4+ramanuj(N/i,k,table)*(double)M[i-1];
	 else
	    sum4=sum4+(ramanuj(N/i,k,table)+1.0)*(double)M[i-1];
	 }
      if ((delta!=0)&&(prime!=0)) {
	 t=(int)sum4;
	 if (t<0)
	    t=-t;
	 if ((unsigned int)t!=((unsigned int)t/k)*k) {
	    printf("error \n",t);
	    goto zskip;
	    }
	 }
      fprintf(Outfp," %e\n",sum4);
      continue;
      }
   if (out==42) {
      sum4=0.0;
      for (i=1; i<=N; i++)
	 if (delta==0)
	    sum4=sum4+ramanuj(N/i,k,table)*log((double)i);
	 else {
	    if (delta==1)
	       sum4=sum4+(ramanuj(N/i,k,table)+1.0)*log((double)i);
	    else {
	       if (delta==2) {
		  t=((i+1)/(k+1)-1)*(k+1);
		  sum4=sum4+(ramanuj(N/i,k,table)+(double)t/(double)i)*log((double)i);
		  }
	       else
		  sum4=sum4+(ramanuj(N/i,k,table)+(i+1)/(k+1))*log((double)i);
	       }
	    }
      fprintf(Outfp," %e\n",sum4);
      continue;
      }
   if (out==43) {
      sum4=0.0;
      for (i=1; i<=N; i++)
	 sum4=sum4+ramanuj(N/i,k,table);
      fprintf(Outfp," %e\n",sum4);
      continue;
      }
   if (out==44) {
      sum3=0.0;
      for (i=1; i<=N; i++) {
	 if (delta==0)
	    sum3=sum3+ramanuj(N/i,k,table)*2.0*numdiv(i,1);
	 else
	    sum3=sum3+(ramanuj(N/i,k,table)+1.0)*2.0*numdiv(i,1);
	 }
      fprintf(Outfp," %e\n",sum3);
      continue;
      }
   if (out==45) {
      sum3=0.0;
      for (i=1; i<=N; i++) {
	 if (delta==0)
	    sum3=sum3+ramanuj(N/i,k,table)*numdiv(i,5);
	 else
	    sum3=sum3+(ramanuj(N/i,k,table)+1.0)*numdiv(i,5);
	 }
      fprintf(Outfp," %e\n",sum3);
      continue;
      }
   if (out==46) {
      sum3=0.0;
      for (i=1; i<=N; i++) {
	 if (switchi==0)
	    gauss(N/i,7,xi,outp);
	 if (switchi==1)
	    gauss(N/i,11,xip,outp);
	 if (switchi==2)
	    gauss(N/i,13,xiq,outp);
	 if (delta==0) {
	    if (reim==0)
	       sum3=sum3+outp[0]*numdiv(i,5);
	    else
	       sum3=sum3+outp[1]*numdiv(i,5);
	    }
	 else {
	    if (reim==0) {
	       if (switchi==0)
		  sum3=sum3+(outp[0]+0.9076)*numdiv(i,5);
	       if (switchi==1)
		  sum3=sum3+(outp[0]+0.4207)*numdiv(i,5);
	       if (switchi==2)
		  sum3=sum3+(outp[0]+1.247)*numdiv(i,5);
	       }
	    else {
	       if (switchi==0)
		  sum3=sum3+(outp[1]-0.8163)*numdiv(i,5);
	       if (switchi==1)
		  sum3=sum3+(outp[1]-2.127)*numdiv(i,5);
	       if (switchi==2)
		  sum3=sum3+(outp[1]-0.08855)*numdiv(i,5);
	       }
	    }
	 }
      fprintf(Outfp," %e\n",sum3);
      continue;
      }
   if (out==47) {
      sum3=0.0;
      for (i=1; i<=N; i++) {
	 if (switchi==0)
	    gauss(N/i,7,xi3,outp);
	 if (switchi==1)
	    gauss(N/i,11,xip,outp);
	 if (switchi==2)
	    gauss(N/i,13,xiq,outp);
	 if (reim==0)
	    sum3=sum3+outp[0];
	 else
	    sum3=sum3+outp[1];
	 }
      fprintf(Outfp," %e\n",sum3);
      continue;
      }
   if (out==48) {
      sum3=0.0;
      for (i=1; i<=N; i++) {
	 if (switchi==0)
	    gauss(N/i,7,xi,outp);
	 if (switchi==1)
	    gauss(N/i,11,xip,outp);
	 if (switchi==2)
	    gauss(N/i,13,xiq,outp);
	 if (delta==0) {
	    if (reim==0)
	       sum3=sum3+outp[0]*log((double)i);
	    else
	       sum3=sum3+outp[1]*log((double)i);
	    }
	 if (delta==1) {
	    if (reim==0) {
	       if (switchi==0)
		  sum3=sum3+(outp[0]+0.9076)*log((double)i);
	       if (switchi==1)
		  sum3=sum3+(outp[0]+0.4207)*log((double)i);
	       if (switchi==2)
		  sum3=sum3+(outp[0]+1.247)*log((double)i);
	       }
	    else {
	       if (switchi==0)
		  sum3=sum3+(outp[1]-0.8163)*log((double)i);
	       if (switchi==1)
		  sum3=sum3+(outp[1]-2.127)*log((double)i);
	       if (switchi==2)
		  sum3=sum3+(outp[1]-0.08855)*log((double)i);
	       }
	    }
	 if (delta==2) {
	    if (reim==0)
	       sum3=sum3+(outp[0]+0.9076)*log((double)i);
	    else
	       sum3=sum3+(outp[1]+0.9076)*log((double)i);
	    }
	 if (delta==3) {
	    if (reim==0)
	       sum3=sum3+(outp[0]-0.8163)*log((double)i);
	    else
	       sum3=sum3+(outp[1]-0.8163)*log((double)i);
	    }
	 }
      fprintf(Outfp," %e\n",sum3);
      continue;
      }
   if (out==49) {
      sum3=0.0;
      for (i=1; i<=N; i++) {
	 if (switchi==0)
	    gauss(N/i,7,xi,outp);
	 if (switchi==1)
	    gauss(N/i,11,xip,outp);
	 if (switchi==2)
	    gauss(N/i,13,xiq,outp);
	 if (delta==0) {
	    if (reim==0)
	       sum3=sum3+outp[0]*numdiv(i,4);
	    else
	       sum3=sum3+outp[1]*numdiv(i,4);
	    }
	 else {
	    if (reim==0) {
	       if (switchi==0)
		  sum3=sum3+(outp[0]+0.9076)*numdiv(i,4);
	       if (switchi==1)
		  sum3=sum3+(outp[0]+0.4207)*numdiv(i,4);
	       if (switchi==2)
		  sum3=sum3+(outp[0]+1.247)*numdiv(i,4);
	       }
	    else {
	       if (switchi==0)
		  sum3=sum3+(outp[1]-0.8163)*numdiv(i,4);
	       if (switchi==1)
		  sum3=sum3+(outp[1]-2.127)*numdiv(i,4);
	       if (switchi==2)
		  sum3=sum3+(outp[1]-0.08855)*numdiv(i,4);
	       }
	    }
	 }
      fprintf(Outfp," %e\n",sum3);
      continue;
      }
   if (out==50) {
      sum4=0.0;
      for (i=1; i<=N; i++) {
	 mertenaw(N/i, count, I);
	 if (switchi==0)
	    gauss(i,7,xi,outp);
	 if (switchi==1)
	    gauss(i,11,xip,outp);
	 if (switchi==2)
	    gauss(i,13,xiq,outp);
	 if (delta==0) {
	    if (reim==0)
	       sum4=sum4+((double)count[1]-(double)count[0])*outp[0];
	    else
	       sum4=sum4+((double)count[1]-(double)count[0])*outp[1];
	    }
	 else {
	    if (reim==0)
	       sum4=sum4+((double)count[0]-(double)count[1]+(double)(J-1)/(double)I)*outp[0];
	    else
	       sum4=sum4+((double)count[0]-(double)count[1]+(double)(J-1)/(double)I)*outp[1];
	    }
	 }
      fprintf(Outfp," %e\n",sum4);
      continue;
      }
   if (out==51) {
      sum3=0.0;
      for (i=1; i<=N; i++) {
	 if (switchi==0)
	    gauss(N/i,7,xi,outp);
	 if (switchi==1)
	    gauss(N/i,11,xip,outp);
	 if (switchi==2)
	    gauss(N/i,13,xiq,outp);
	 if (delta==0) {
	    if (reim==0)
	       sum3=sum3+outp[0]*2.0*numdiv(i,1);
	    else
	       sum3=sum3+outp[1]*2.0*numdiv(i,1);
	    }
	 else {
	    if (reim==0) {
	       if (switchi==0)
		  sum3=sum3+(outp[0]+0.9076)*2.0*numdiv(i,1);
	       if (switchi==1)
		  sum3=sum3+(outp[0]+0.4207)*2.0*numdiv(i,1);
	       if (switchi==2)
		  sum3=sum3+(outp[0]+1.247)*2.0*numdiv(i,1);
	       }
	    else {
	       if (switchi==0)
		  sum3=sum3+(outp[1]-0.8163)*2.0*numdiv(i,1);
	       if (switchi==1)
		  sum3=sum3+(outp[1]-2.127)*2.0*numdiv(i,1);
	       if (switchi==2)
		  sum3=sum3+(outp[1]-0.08855)*2.0*numdiv(i,1);
	       }
	    }
	 }
      fprintf(Outfp," %e\n",sum3);
      continue;
      }
   if (out==52) {
      sum3=0.0;
      for (i=1; i<=N; i++) {
	 if (switchi==0)
	    gauss(N/i,7,xi,outp);
	 if (switchi==1)
	    gauss(N/i,11,xip,outp);
	 if (switchi==2)
	    gauss(N/i,13,xiq,outp);
	 if (delta==0) {
	    if (reim==0)
	       sum3=sum3+outp[0]*log((double)i)*numdiv(i,1);
	    else
	       sum3=sum3+outp[1]*log((double)i)*numdiv(i,1);
	    }
	 else {
	    if (reim==0) {
	       if (switchi==0)
		  sum3=sum3+(outp[0]+0.9076)*log((double)i)*numdiv(i,1);
	       if (switchi==1)
		  sum3=sum3+(outp[0]+0.4207)*log((double)i)*numdiv(i,1);
	       if (switchi==2)
		  sum3=sum3+(outp[0]+1.247)*log((double)i)*numdiv(i,1);
	       }
	    else {
	       if (switchi==0)
		  sum3=sum3+(outp[1]-0.8163)*log((double)i)*numdiv(i,1);
	       if (switchi==1)
		  sum3=sum3+(outp[1]-2.127)*log((double)i)*numdiv(i,1);
	       if (switchi==2)
		  sum3=sum3+(outp[1]-0.08855)*log((double)i)*numdiv(i,1);
	       }
	    }
	 }
      fprintf(Outfp," %e\n",sum3);
      continue;
      }
   if (out==53) {
      sum3=0.0;
      for (i=1; i<=N; i++) {
	 if (switchi==0)
	    gauss(N/i,7,xi,outp);
	 else
	    gauss(N/i,11,xip,outp);
	 if (delta==0) {
	    if (reim==0)
	       sum3=sum3+outp[0]*mangoldt(i,table);
	    else
	       sum3=sum3+outp[1]*mangoldt(i,table);
	    }
	 else {
	    if (reim==0) {
	       if (switchi==0)
		  sum3=sum3+(outp[0]+0.9076)*mangoldt(i,table);
	       else
		  sum3=sum3+(outp[0]+0.4207)*mangoldt(i,table);
	       }
	    else {
	       if (switchi==0)
		  sum3=sum3+(outp[1]-0.8163)*mangoldt(i,table);
	       else
		  sum3=sum3+(outp[1]-2.127)*mangoldt(i,table);
	       }
	    }
	 }
      fprintf(Outfp," %e\n",sum3);
      continue;
      }
   if (out==54) {
      sum3=0.0;
      for (i=1; i<=N; i++) {
	 if (switchi==0)
	    gauss(N/i,7,xi,outp);
	 else
	    gauss(N/i,11,xip,outp);
	 if (delta==0) {
	    if (reim==0)
	       sum3=sum3+outp[0]*(double)M[i-1];
	    else
	       sum3=sum3+outp[1]*(double)M[i-1];
	    }
	 else {
	    if (reim==0) {
	       if (switchi==0)
		  sum3=sum3+(outp[0]+0.9076)*M[i-1];
	       else
		  sum3=sum3+(outp[0]+0.4207)*M[i-1];
	       }
	    else {
	       if (switchi==0)
		  sum3=sum3+(outp[1]-0.8163)*M[i-1];
	       else
		  sum3=sum3+(outp[1]-2.127)*M[i-1];
	       }
	    }
	 }
      fprintf(Outfp," %e\n",sum3);
      continue;
      }
   if (out==55) {
      sum4=0.0;
      for (i=1; i<=N; i++) {
	 mertenaw(i, count, I);
	 if (delta==0)
	    sum4=sum4+ramanuj(N/i,k,table)*((double)count[0]-(double)count[1]);
	 else
	    sum4=sum4+(ramanuj(N/i,k,table)+1.0)*((double)count[0]-(double)count[1]);
	 }
      fprintf(Outfp," %e\n",sum4);
      continue;
      }
   if (out==56) {
      sum3=0.0;
      for (i=1; i<=N; i++) {
	 mertenaw(i, count, I);
	 if (switchi==0)
	    gauss(N/i,7,xi,outp);
	 else
	    gauss(N/i,11,xip,outp);
	 if (delta==0) {
	    if (reim==0)
	       sum3=sum3+outp[0]*((double)count[0]-(double)count[1]);
	    else
	       sum3=sum3+outp[1]*((double)count[0]-(double)count[1]);
	    }
	 else {
	    if (reim==0) {
	       if (switchi==0)
		  sum3=sum3+(outp[0]+0.9076)*((double)count[0]-(double)count[1]);
	       else
		  sum3=sum3+(outp[0]+0.4207)*((double)count[0]-(double)count[1]);
	       }
	    else {
	       if (switchi==0)
		  sum3=sum3+(outp[1]-0.8163)*((double)count[0]-(double)count[1]);
	       else
		  sum3=sum3+(outp[1]-2.127)*((double)count[0]-(double)count[1]);
	       }
	    }
	 }
      fprintf(Outfp," %e\n",sum3);
      continue;
      }
   if (out==57) {
      if (switchi==0)
	 gauss(N,7,xi,outp);
      if (switchi==1)
	 gauss(N,11,xip,outp);
      if (switchi==2)
	 gauss(N,13,xiq,outp);
      if (reim==0)
	 sum3=outp[0];
      else
	 sum3=outp[1];
      fprintf(Outfp," %e\n",sum3);
      continue;
      }
   if (out==58) {
      sum4=0.0;
      for (i=1; i<=N; i++) {
	 if (delta==0)
	    sum4=sum4+ramanuj(N/i,k,table)*(double)M[i-1]*(double)M[i-1];
	 else
	    sum4=sum4+(ramanuj(N/i,k,table)+1.0)*(double)M[i-1]*(double)M[i-1];
	 }
      if ((delta!=0)&&(prime!=0)) {
	 t=(int)sum4;
	 if (t<0) {
	    printf("errror 0: \n",t);
	    goto zskip;
	    }
	 if (N<=(k*k)) {
	    if ((unsigned int)t>=maxt) {
	       savet[3]=savet[2];
	       savet[2]=savet[1];
	       savet[1]=savet[0];
	       savet[0]=t;
	       maxt=(unsigned int)t;
	       saveN[3]=saveN[2];
	       saveN[2]=saveN[1];
	       saveN[1]=saveN[0];
	       saveN[0]=N;
	       }
	    }
	 if (N==(k*k)) {
	    printf(" N=%d %d %d %d \n",saveN[0],saveN[1],saveN[2],saveN[3]);
	    printf(" max=%d %d %d %d \n",savet[0]/k,savet[1]/k,savet[2]/k,savet[3]/k);
	    break;
	    }
	 }
      fprintf(Outfp," %e\n",sum4);
      continue;
      }
   if (out==59) {
      sum4=ramanuj(N,k,table);
      fprintf(Outfp," %e\n",sum4);
      continue;
      }
   if (out==60) {
      if (N!=(N/k)*k)
	 continue;
      sum4=0.0;
      for (i=1; i<=N; i++) {
	 if ((N/i)==((N/i)/k)*k)
	    sum4=sum4+(double)M[i-1]*(double)M[i-1];
	 }
      t=(int)sum4;
      if (N<=(k*k)) {
	 if ((unsigned int)t>=maxt) {
	    savet[3]=savet[2];
	    savet[2]=savet[1];
	    savet[1]=savet[0];
	    savet[0]=t;
	    maxt=(unsigned int)t;
	    saveN[3]=saveN[2];
	    saveN[2]=saveN[1];
	    saveN[1]=saveN[0];
	    saveN[0]=N;
	    }
	 }
      if (N==(k*k)) {
	 printf(" N=%d %d %d %d \n",saveN[0],saveN[1],saveN[2],saveN[3]);
	 printf(" max=%d %d %d %d \n",savet[0],savet[1],savet[2],savet[3]);
	 break;
	 }
      fprintf(Outfp," %e\n",sum4);
      continue;
      }
   if (out==61) {
      sum4=0.0;
      for (i=1; i<=N; i++) {
	 if (N==(N/i)*i)
	    sum4=sum4+(double)M[N/i-1]*(double)M[N/i-1];
	 }
      if (mflag!=0) {
	 t=(int)sum4;
	 if (N<k) {
	    if ((unsigned int)t>=maxt) {
	       savet[3]=savet[2];
	       savet[2]=savet[1];
	       savet[1]=savet[0];
	       savet[0]=t;
	       maxt=(unsigned int)t;
	       saveN[3]=saveN[2];
	       saveN[2]=saveN[1];
	       saveN[1]=saveN[0];
	       saveN[0]=N;
	       }
	    }
	 if (N==k) {
	    printf(" N=%d %d %d %d \n",saveN[0],saveN[1],saveN[2],saveN[3]);
	    printf(" max=%d %d %d %d \n",savet[0],savet[1],savet[2],savet[3]);
	    break;
	    }
	 }
      fprintf(Outfp," %e\n",sum4);
      continue;
      }
   if (out==62) {
      sum4=0.0;
      for (i=1; i<=N; i++) {
	 if (N==(N/i)*i)
	    sum4=sum4+(double)M[N/i-1]*(double)M[N/i-1];
	 }
      if (mflag==4)
	 fprintf(Outfp," %e \n",sum4);
      else {
	 t=(int)sum4;
	 if ((unsigned int)t>maxt) {
	    maxt=(unsigned int)t;
	    if ((N>3)&&((double)maxt<=(double)N/log((double)N))) {
	       printf("error 62: %d %d %e \n",N,maxt,(double)N/log((double)N));
	       goto zskip;
	       }
	    printf(" %d %d \n",N,maxt);
	    if (mflag==0)
	       fprintf(Outfp," %d \n",maxt);
	    if (mflag==1)
	       fprintf(Outfp," %d \n",N);
	    if (mflag==2)
	       fprintf(Outfp," %e \n",(double)N/log((double)N)/(double)maxt);
	    if (mflag==3)
	       fprintf(Outfp," %e \n",(double)maxt/(double)N);
	    }
	 }
      continue;
      }
   if (out==63) {
      sum4=0.0;
      for (i=1; i<=N; i++)
	 sum4=sum4+(double)M[N/i-1]*(double)M[N/i-1];
      fprintf(Outfp," %e \n",sum4);
      continue;
      }
   }
//
// additional processing
//
if ((out==13)&&(mflag!=0)) {
   kp=I;
   if (MAXN<=kp) {
      printf("MAXN too small \n");
      goto zskip;
      }
   printf("\n");
   outman[0]=outman[kp];
   for (N=1; N<=kp; N++) {
      sum4=0.0;
      sum5=0.0;
      temp1=-2.0*pi*(double)(N-1)/(double)kp;
      for (i=0; i<kp; i++) {
	 sum4=sum4+outman[i]*cos((double)i*temp1);
	 sum5=sum5+outman[i]*sin((double)i*temp1);
	 }
      printf(" %d %e %e \n",N,sum4/(double)kp,sum5/(double)kp);
      if (mflag==1) {
	 if (reim==0)
	    fprintf(Outfp," %e\n",sum4/(double)kp);
	 else
	    fprintf(Outfp," %e\n",sum5/(double)kp);
	 }
      else {
	 gre[N-1]=sum4/(double)kp;
	 gim[N-1]=sum5/(double)kp;
	 }
      }
   if (mflag==1)
      goto zskip;
   printf("\n");
   for (N=1; N<=kp; N++) {
      sum4=0.0;
      sum5=0.0;
      temp1=2.0*pi*(double)(N-1)/(double)kp;
      for (i=0; i<kp; i++) {
	 sum4=sum4+gre[i]*cos((double)i*temp1)-gim[i]*sin((double)i*temp1);
	 sum5=sum5+gim[i]*cos((double)i*temp1)+gre[i]*sin((double)i*temp1);
	 }
      printf(" %d %e %e \n",N,sum4,sum5);
      if (reim==0)
	 fprintf(Outfp," %e\n",sum4);
      else
	 fprintf(Outfp," %e\n",sum5);
      }
   goto zskip;
   }
if ((out==38)&&(mflag!=0)) {
   kp=k*I;
   outman[0]=outman[kp];
   if (mflag==3)
      goto askip;
   if (mflag==4)
      goto bskip;
   if (MAXN<=kp) {
      printf("MAXN too small \n");
      goto zskip;
      }
   printf("\n");
   resum=0.0;
   imsum=0.0;
   for (N=1; N<=kp; N++) {
      sum4=0.0;
      sum5=0.0;
      temp1=-2.0*pi*(double)(N-1)/(double)kp;
      for (i=0; i<kp; i++) {
	 sum4=sum4+outman[i]*cos((double)i*temp1);
	 sum5=sum5+outman[i]*sin((double)i*temp1);
	 }
      printf(" %d %e %e \n",N,sum4/(double)kp,sum5/(double)kp);
      if (mflag==1) {
	 if (reim==0)
	    fprintf(Outfp," %e\n",sum4/(double)kp);
	 else
	    fprintf(Outfp," %e\n",sum5/(double)kp);
	 }
      else {
	 gre[N-1]=sum4/(double)kp;
	 gim[N-1]=sum5/(double)kp;
	 }
      resum=resum+sum4;
      imsum=imsum+sum5;
      }
   printf(" sum=%e %e \n",resum/(double)kp,imsum/(double)kp);
   if (mflag==1)
      goto zskip;
   printf("\n");
   for (N=1; N<=kp; N++) {
      sum4=0.0;
      sum5=0.0;
      temp1=2.0*pi*(double)(N-1)/(double)kp;
      for (i=0; i<kp; i++) {
	 sum4=sum4+gre[i]*cos((double)i*temp1)-gim[i]*sin((double)i*temp1);
	 sum5=sum5+gim[i]*cos((double)i*temp1)+gre[i]*sin((double)i*temp1);
	 }
      printf(" %d %e %e \n",N,sum4,sum5);
      if (reim==0)
	 fprintf(Outfp," %e\n",sum4);
      else
	 fprintf(Outfp," %e\n",sum5);
      }
   goto zskip;
//
// Ramanujan's sum
//
askip:
//   printf(" %e \n",outman[0]);
//   goto zskip;
   for (N=2; N<=MAXN; N++) {
      sum4=0.0;
      for (i=1; i<=N; i++) {
	 if (delta==0)
	    sum4=sum4+ramanuj(N/i,k,table)*outman[i-1];
	 else
	    sum4=sum4+(ramanuj(N/i,k,table)+1.0)*outman[i-1];
	 }
      fprintf(Outfp," %e\n",sum4);
      }
   goto zskip;
//
// Gauss sum
//
bskip:
   for (N=2; N<=MAXN; N++) {
      sum3=0.0;
      for (i=1; i<=N; i++) {
	 if (switchi==0)
	    gauss(N/i,7,xi,outp);
	 if (switchi==1)
	    gauss(N/i,11,xip,outp);
	 if (switchi==2)
	    gauss(N/i,13,xiq,outp);
	 if (delta==0) {
	    if (reim==0)
	       sum3=sum3+outp[0]*outman[i-1];
	    else
	       sum3=sum3+outp[1]*outman[i-1];
	    }
	 else {
	    if (reim==0) {
	       if (switchi==0)
		  sum3=sum3+(outp[0]+0.9076)*outman[i-1];
	       if (switchi==1)
		  sum3=sum3+(outp[0]+0.4207)*outman[i-1];
	       if (switchi==2)
		  sum3=sum3+(outp[0]+1.247)*outman[i-1];
	       }
	    else {
	       if (switchi==0)
		  sum3=sum3+(outp[1]-0.8163)*outman[i-1];
	       if (switchi==1)
		  sum3=sum3+(outp[1]-2.127)*outman[i-1];
	       if (switchi==2)
		  sum3=sum3+(outp[1]-0.08855)*outman[i-1];
	       }
	    }
	 }
      fprintf(Outfp," %e\n",sum3);
      }
   goto zskip;
   }
if ((out==37)&&(mflag!=0)) {
   outman[0]=0.0;
   output[0]=0.0;
   for (N=2; N<=MAXN; N++) {
      sum4=0.0;
      for (i=1; i<=N; i++) {
	 if (delta==0) {
	    if (mflag==1)
	       sum4=sum4+outman[N/i-1]*mangoldt(i,table);
	    if (mflag==2)
	       sum4=sum4+outman[i-1]*mangoldt(N/i,table);
	    if (mflag==3)
	       sum4=sum4+outman[N/i-1]*output[i-1];
	    if (mflag==4)
	       sum4=sum4+output[N/i-1]*outman[i-1];
	    }
	 else {
	    if (mflag==1)
	       sum11=sum11+(outman[N/i-1]+temp)*mangoldt(i,table);
	    else
	       sum11=sum11+(outman[i-1]+temp)*mangoldt(N/i,table);
	    }
//	 printf(" %d %e %e \n",i,outman[i-1],mangoldt(N/i,table));
	 }
      printf(" %d %e \n",N,sum4);
      if (again==0)
	 fprintf(Outfp," %e\n",sum4);
      else
	 manout[N-1]=sum4;
      }
   if (again==0)
      goto zskip;
   manout[0]=0.0;
   for (N=2; N<=MAXN; N++) {
      sum4=0.0;
      for (i=1; i<=N; i++) {
	 if (delta==0) {
	    if (mflag==1)
	       sum4=sum4+manout[N/i-1]*mangoldt(i,table);
	    if (mflag==2)
	       sum4=sum4+manout[i-1]*mangoldt(N/i,table);
	    if (mflag==3)
	       sum4=sum4+manout[N/i-1]*output[i-1];
	    if (mflag==4)
	       sum4=sum4+output[N/i-1]*manout[i-1];
	    }
	 else {
	    if (mflag==1)
	       sum11=sum11+(manout[N/i-1]+temp)*mangoldt(i,table);
	    else
	       sum11=sum11+(manout[i-1]+temp)*mangoldt(N/i,table);
	    }
//	 printf(" %d %e %e \n",i,manout[i-1],mangoldt(N/i,table));
	 }
      printf(" %d %e \n",N,sum4);
      fprintf(Outfp," %e\n",sum4);
      }
   }
if ((out==17)&&(hflag==1)) {
   mean=0.0;
   msum2=0.0;
   for (i=0; i<301; i++) {
      fprintf(Outfp," %d\n",histo[i]);
      mean=mean+((double)i-150.0)*histo[i];
//    printf(" %d %d %d %e\n",(int)(i-150),histo[i],(int)(i-150)*histo[i],mean);
      msum2=msum2+((double)i-150.0)*histo[i]*((double)i-150.0);
      }
   msum=mean;
   mean=mean/(double)(MAXN-1);
   std=(double)(MAXN-1)*msum2-msum*msum;
   std=sqrt(std/((double)(MAXN-1)*(double)(MAXN-2)));
   printf("mean=%e, std=%e \n",mean,std);
   }
zskip:
if ((out==37)||(out==38)||(out==13)||(out==50)) {
   free(output);
   free(outman);
   free(manout);
   }
if ((out==34)||(out==39)||(out==41)||(out==54)||(out==58)||(out==60)||(out==61)||(out==62)||(out==63))
   free(M);
fclose(Outfp);
return;
}