/*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;
}