/*CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C C
C COMPUTE UPPER BOUND MINUS LOWER BOUND (AND ESTIMATED VALUES) C
C 08/08/14 (DKC) C
C C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC*/
#include <stdio.h>
#include <math.h>
unsigned int euler(unsigned int a);
void uplolim(unsigned int order, unsigned int *sums);
double mertens6(unsigned int, unsigned int *count, double *sum);
//
unsigned int out=0; // 0, sum of M([x/i]) when 4 divides phi(i)
// 1, estimated sum of M([x/i]) when 4 divides phi(i)
// 2, exceptions (when estimated signs not correct)
// 3, sum of a(x)*m(x/i)/m(x)
// 4, sum of b(x)*n(x/i)/n(x)
// 5, sum of i*a(x)*m(x/i)/m(x)
// 6, sum of i*b(x)*n(x/i)/n(x)
// 7, sum of a(x/i) (set limit to 3)
// 8, sum of b(x/i) (set limit to 3)
// 9, sum of i*a(x/i) (set limit to 3)
// 10, sum of i*b(x/i) (set limit to 3)
// 11, sum of m(x/i)
// 12, sum of n(x/i)
// 13, sum of m(x/i)-n(x/i)
// 14, sum of i*(m(x/i)-n(x/i))
// 15, 1/2 sum of i*M([x/i])
// 16, Mertens function
// 17, differences in number of fractions before, after 1/4
// 18, number of fractions before 1/4
// 19, number of fractions after 1/4
// 20, lower limit of number of fractions
// 21, upper limit of number of fractions
//
double inc=0.0; // esum+increment
unsigned int base=1; // Mertens or 0
unsigned int n=1; // decimation factor
unsigned int limit=0; // skip M values less than or equal to limit
// set to 3 for out=7,8,9,and 10, other values for out=15
//
void main() {
int k;
unsigned int N,MAXN; // MAXN=5128
unsigned int i,j,temp,clsave,crsave,count[2],histop[100],histopn[100];
unsigned int histon[100],histonn[100],histox,histoxn,flag,sums[2];
double sum,sum1,msum,mlsave,mrsave,lratio,rratio,tsum[2],esum,tmpsum;
double lsum,rsum,lsum2,rsum2,sum2,sum3,sum4,sum5,sum6,sum7,sum8,sum9;
double newsum,oldsum,ftemp,upper,lower;
FILE *Outfp;
Outfp = fopen("out1bc.dat","w");
for (i=0; i<100; i++) {
histop[i]=0;
histopn[i]=0;
histon[i]=0;
histonn[i]=0;
}
sums[0]=0;
sums[1]=0;
histox=0;
histoxn=0;
MAXN=5128; // 5128
oldsum=0.0;
for (N=2; N<=MAXN; N++) {
uplolim(N,sums); // compute upper and lower limits of number of fractions
if (out==20)
fprintf(Outfp," %d\n",sums[0]);
if (out==21)
fprintf(Outfp," %d\n",sums[1]);
sum=0.0;
sum1=0.0;
sum2=0.0;
sum3=0.0;
sum4=0.0;
sum5=0.0;
sum6=0.0;
sum7=0.0;
sum8=0.0;
sum9=0.0;
newsum=0.0;
for (i=1; i<=N; i++) {
if ((N/i)<=limit)
break;
msum=mertens6(N/i, count, tsum);
if ((out==16)&&(i==1))
fprintf(Outfp," %e\n",msum);
if ((out==17)&&(i==1))
fprintf(Outfp," %d\n",(int)count[0]-(int)count[1]);
if ((out==18)&&(i==1)) {
fprintf(Outfp," %d\n",count[0]);
if ((count[0]>sums[1])||(count[0]<sums[0])) {
printf("error: N=%d, count=%d, lo=%d, hi=%d \n",N,count[0],sums[0],sums[1]);
goto zskip;
}
}
if ((out==19)&&(i==1)) {
fprintf(Outfp," %d\n",count[1]);
if ((count[1]>sums[1])||(count[1]<sums[0])) {
printf("error: N=%d, count=%d, lo=%d, hi=%d \n",N,count[0],sums[0],sums[1]);
goto zskip;
}
}
if (i==1) {
mlsave=tsum[0];
mrsave=tsum[1];
clsave=count[0];
crsave=count[1];
if (base==0)
esum=msum;
else {
esum=0.0;
lsum=mlsave;
rsum=mrsave;
lsum2=mlsave;
rsum2=mrsave;
}
// printf(" %d %d %e %e \n",clsave,crsave,mlsave,mrsave);
}
//
// compute sums
//
temp=euler(i);
if ((temp/4)*4!=temp) {
sum=sum+msum*(double)i/(double)temp;
sum1=sum1+esum*(double)i/(double)temp;
}
sum2=sum2+tsum[0];
sum3=sum3+tsum[1];
sum4=sum4+tsum[0]*(double)i;
sum5=sum5+tsum[1]*(double)i;
tmpsum=tsum[0]+tsum[1];
if (tmpsum>=0.0)
k=(int)(tmpsum*2.0+0.001);
else
k=(int)(tmpsum*2.0-0.001);
newsum=newsum+(double)(k/2.0)*(double)i;
sum6=sum6+(double)count[0];
sum7=sum7+(double)count[1];
sum8=sum8+(double)count[0]-(double)count[1];
sum9=sum9+(double)i*((double)count[0]-(double)count[1]);
// printf(" %d %d %d %e %e %e \n",i,count[0],N/i,sum4,sum5,2.0*(tsum[0]+tsum[1]));
// printf(" %d %d %d %e %e \n",i,count[0],N/i,sum,msum);
if (i!=1) {
if (clsave!=0)
lratio=(double)count[0]/(double)clsave;
else
lratio=1.0;
if (crsave!=0)
rratio=(double)count[1]/(double)crsave;
else
rratio=1.0;
esum=mlsave*lratio+mrsave*rratio;
esum=esum*2.0;
if ((i-1)==((i-1)/n)*n) {
lsum=lsum+mlsave*lratio;
rsum=rsum+mrsave*rratio;
}
lsum2=lsum2+(double)i*mlsave*lratio;
rsum2=rsum2+(double)i*mrsave*rratio;
// printf(" %e %e %e \n",lsum,rsum,lsum+rsum);
// fprintf(Outfp," i=%d %d %d %e %e %e %e %e %e \n",i,count[0],count[1],tsum[0],tsum[1],
// lratio,rratio,msum,esum);
//
// histogram (for msum>0)
//
if ((msum>0.0)&&(i<=(N/2))) {
j=(unsigned int)(msum+0.0001);
if ((j<100)&&(j!=0)) {
histop[j]=histop[j]+1;
if (esum+inc>0.0)
histopn[j]=histopn[j]+1;
}
printf("N=%d, i=%d, N/i=%d, msum=%e, esum=%e \n",N,i,N/i,msum,esum);
if (out==2)
fprintf(Outfp,"N=%d, i=%d, N/i=%d, msum=%e, esum=%e \n",N,i,N/i,msum,esum);
}
//
// histogram (for msum<0)
//
if (msum<0.0) {
tmpsum=-msum;
j=(unsigned int)(tmpsum+0.0001);
if ((j<100)&&(j!=0)) {
histon[j]=histon[j]+1;
if (esum-inc<0.0)
histonn[j]=histonn[j]+1;
}
}
//
// histogram (for msum==0)
//
if (msum<0.0)
tmpsum=-msum;
else
tmpsum=msum;
j=(unsigned int)(tmpsum+0.0001);
if ((j==0)&&(i<=(N/3))) {
histox=histox+1;
if (esum<0.0)
tmpsum=-esum;
else
tmpsum=esum;
j=(unsigned int)(tmpsum+0.1);
if (j==0)
histoxn=histoxn+1;
}
}
}
if (limit==0) {
if (newsum<=oldsum) {
printf("error 1: N=%d, newsum=%e, oldsum=%e \n",N,newsum,oldsum);
goto zskip;
}
ftemp=sqrt(newsum);
upper=0.3898*(double)N+0.1952+0.50;
lower=upper-1.00;
if ((ftemp>upper)||(ftemp<lower)) {
printf("error: N=%d, newsum=%e, upper=%e, lower=%e \n",N,newsum,upper,lower);
goto zskip;
}
}
if (limit==1) {
if (newsum>oldsum) {
printf("error 1: N=%d, newsum=%e, oldsum=%e \n",N,newsum,oldsum);
goto zskip;
}
if (newsum==oldsum) {
j=N;
while ((j&1)==0)
j=j/2;
if (j!=1) {
printf("error 2: N=%d, newsum=%e, oldsum=%e \n",N,newsum,oldsum);
if ((N!=5107)&&(N!=5113))
goto zskip;
}
}
ftemp=sqrt(-newsum);
upper=0.1885*(double)N+0.09297+1.00;
lower=upper-2.00;
if ((ftemp>upper)||(ftemp<lower)) {
printf("error: N=%d, newsum=%e, upper=%e, lower=%e \n",N,newsum,upper,lower);
goto zskip;
}
}
if (limit==3) {
if (((int)(N-6)/12)*12==(int)(N-6)) {
if (newsum<=oldsum) {
printf("error 1: N=%d, newsum=%e, oldsum=%e \n",N,newsum,oldsum);
goto zskip;
}
}
else {
if (newsum>oldsum) {
printf("error 2: N=%d, newsum=%e, oldsum=%e \n",N,newsum,oldsum);
goto zskip;
}
}
ftemp=sqrt(-newsum);
upper=0.1529*(double)N+0.0746+1.00;
lower=upper-2.00;
if ((ftemp>upper)||(ftemp<lower)) {
printf("error: N=%d, newsum=%e, upper=%e, lower=%e \n",N,newsum,upper,lower);
goto zskip;
}
}
if (limit==4) {
if (((int)N/6)*6==(int)N) {
if (newsum<=oldsum) {
if (N!=(N/5)*5) {
printf("error 1: N=%d, newsum=%e, oldsum=%e \n",N,newsum,oldsum);
goto zskip;
}
}
}
else {
if (newsum>oldsum) {
if (N!=8) {
printf("error 2: N=%d, newsum=%e, oldsum=%e \n",N,newsum,oldsum);
goto zskip;
}
}
}
ftemp=sqrt(-newsum);
upper=0.1332*(double)N+0.06404+1.00;
lower=upper-2.00;
if ((ftemp>upper)||(ftemp<lower)) {
printf("error: N=%d, newsum=%e, upper=%e, lower=%e \n",N,newsum,upper,lower);
goto zskip;
}
}
if (limit==5) {
flag=0;
if ((N&1)==0) {
if (((int)(N+5)==((int)(N+5)/15)*15)||((int)(N-5)==((int)(N-5)/15)*15))
flag=1;
}
if (((int)(N/15)*15==(int)(N))&&((N&1)==1)||(flag==1)) {
if (newsum<=oldsum) {
if (N!=(N/7)*7) {
printf("error 1: N=%d, newsum=%e, oldsum=%e \n",N,newsum,oldsum);
if ((N!=2145)&&(N!=2805)&&(N!=3135)&&(N!=3315)&&(N!=3705)&&(N!=3795))
goto zskip;
}
}
}
else {
if (newsum>oldsum) {
if (N!=8) {
printf("error 2: N=%d, newsum=%e, oldsum=%e \n",N,newsum,oldsum);
goto zskip;
}
}
}
ftemp=sqrt(-newsum);
upper=0.1078*(double)N+0.05121+1.00;
lower=upper-2.00;
if ((ftemp>upper)||(ftemp<lower)) {
printf("error: N=%d, newsum=%e, upper=%e, lower=%e \n",N,newsum,upper,lower);
goto zskip;
}
}
if (limit==6) {
flag=0;
if (((int)N==((int)N/3)*3)||((N&1)==0))
flag=1;
if (N==(N/7)*7)
flag=0;
if ((((int)N/5)*5==(int)N)&&(flag!=0)) {
if (newsum<=oldsum) {
if (N!=5) {
printf("error 1: N=%d, newsum=%e, oldsum=%e \n",N,newsum,oldsum);
if ((N!=2145)&&(N!=2805)&&(N!=3135)&&(N!=3315)&&(N!=3705)&&(N!=3795))
goto zskip;
}
}
}
else {
if (newsum>oldsum) {
if (N!=8) {
printf("error 2: N=%d, newsum=%e, oldsum=%e \n",N,newsum,oldsum);
goto zskip;
}
}
}
ftemp=sqrt(-newsum);
upper=0.09893*(double)N+0.04607+1.00;
lower=upper-2.00;
if ((ftemp>upper)||(ftemp<lower)) {
printf("error: N=%d, newsum=%e, upper=%e, lower=%e \n",N,newsum,upper,lower);
goto zskip;
}
}
if (limit==7) {
ftemp=sqrt(-newsum);
upper=0.086*(double)N+0.03871+1.00;
lower=upper-2.00;
if ((ftemp>upper)||(ftemp<lower)) {
printf("error: N=%d, newsum=%e, upper=%e, lower=%e \n",N,newsum,upper,lower);
goto zskip;
}
}
if (limit==8) {
ftemp=sqrt(-newsum);
upper=0.07587*(double)N+0.0333+1.00;
lower=upper-2.00;
if ((ftemp>upper)||(ftemp<lower)) {
printf("error: N=%d, newsum=%e, upper=%e, lower=%e \n",N,newsum,upper,lower);
goto zskip;
}
}
if (limit==9) {
ftemp=sqrt(-newsum);
upper=0.0677*(double)N+0.02939+1.00;
lower=upper-2.00;
if ((ftemp>upper)||(ftemp<lower)) {
printf("error: N=%d, newsum=%e, upper=%e, lower=%e \n",N,newsum,upper,lower);
goto zskip;
}
}
if (limit==10) {
ftemp=sqrt(-newsum);
upper=0.06441*(double)N+0.02691+1.00;
lower=upper-2.00;
if ((ftemp>upper)||(ftemp<lower)) {
printf("error: N=%d, newsum=%e, upper=%e, lower=%e \n",N,newsum,upper,lower);
goto zskip;
}
}
if (limit==11) {
ftemp=sqrt(-newsum);
upper=0.05907*(double)N+0.02301+1.00;
lower=upper-2.00;
if ((ftemp>upper)||(ftemp<lower)) {
printf("error: N=%d, newsum=%e, upper=%e, lower=%e \n",N,newsum,upper,lower);
goto zskip;
}
}
if (limit==12) {
ftemp=sqrt(-newsum);
upper=0.05455*(double)N+0.01939+1.00;
lower=upper-2.00;
if ((ftemp>upper)||(ftemp<lower)) {
printf("error: N=%d, ftemp=%e, upper=%e, lower=%e \n",N,ftemp,upper,lower);
if ((N!=2263)&&(N!=4199))
goto zskip;
}
}
if (limit==13) {
ftemp=sqrt(-newsum);
upper=0.04862*(double)N+0.01661+1.00;
lower=upper-2.00;
if ((ftemp>upper)||(ftemp<lower)) {
printf("error: N=%d, newsum=%e, upper=%e, lower=%e \n",N,newsum,upper,lower);
goto zskip;
}
}
if (limit==14) {
ftemp=sqrt(-newsum);
upper=0.04511*(double)N+0.01567+1.00;
lower=upper-2.00;
if ((ftemp>upper)||(ftemp<lower)) {
printf("error: N=%d, newsum=%e, upper=%e, lower=%e \n",N,newsum,upper,lower);
goto zskip;
}
}
if (limit==15) {
ftemp=sqrt(-newsum);
upper=0.0436*(double)N+0.01452+1.00;
lower=upper-2.00;
if ((ftemp>upper)||(ftemp<lower)) {
printf("error: N=%d, newsum=%e, upper=%e, lower=%e \n",N,newsum,upper,lower);
goto zskip;
}
}
oldsum=newsum;
sum=-sum;
sum1=-sum1;
if (out==0)
fprintf(Outfp," %e\n",sum);
else {
if (out==1)
fprintf(Outfp," %e\n",sum1);
else {
if (out==3)
fprintf(Outfp," %e\n",lsum);
if (out==4)
fprintf(Outfp," %e\n",rsum);
if (out==5)
fprintf(Outfp," %e\n",lsum2);
if (out==6)
fprintf(Outfp," %e\n",rsum2);
if (out==7)
fprintf(Outfp," %e\n",sum2);
if (out==8)
fprintf(Outfp," %e\n",sum3);
if (out==9)
fprintf(Outfp," %e\n",sum4);
if (out==10)
fprintf(Outfp," %e\n",sum5);
if (out==11)
fprintf(Outfp," %e\n",sum6);
if (out==12)
fprintf(Outfp," %e\n",sum7);
if (out==13)
fprintf(Outfp," %e\n",sum8);
if (out==14)
fprintf(Outfp," %e\n",sum9);
if (out==15)
fprintf(Outfp," %e\n",newsum);
}
}
}
if (out==2)
fprintf(Outfp,"\n");
printf(" %d, %d, %d \n",0,histox,histoxn);
if (out==2)
fprintf(Outfp," %d, %d, %d \n",0,histox,histoxn);
for (i=1; i<20; i++) {
printf(" %d, %d %d %d %d \n",i,histop[i],histopn[i],histon[i],histonn[i]);
if (out==2)
fprintf(Outfp," %d, %d %d %d %d \n",i,histop[i],histopn[i],histon[i],histonn[i]);
}
zskip:
fclose(Outfp);
return;
}