﻿ compute Mertens function, etc.
```/*CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                             C
C  COMPUTE MERTENS FUNCTION (sum of |sgn(M([x/i]))| or |sgn(M([x/i]))|i)      C
C  12/21/14 (DKC)							      C
C                                                                             C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC*/
#include <math.h>
#include <stdio.h>
unsigned int out=1;  // 0 for Mertens, 1 for sum function, 2 for histogram
unsigned int max=1000; // maximum x value
unsigned int n=1;	// decimation factor (normally set to 1)
unsigned int bins=200;	// number of histogram bins (for out=2)
unsigned int flag1=0;	// select 1 or i
unsigned int sflag=0;	// square root flag
unsigned int limit=0;	// partial sum flag, normally set to 0
void main() {
short m[500000];
int sum,t;
unsigned int i,j,index,hisdif[1001],width,savsum,savind,numone;
double delta,newsum,oldsum,maxdel,maxsav[10],temp;
FILE *Outfp;
Outfp = fopen("out7sx.dat","w");
if (max>500000) {
printf("x too large \n");
goto zskip;
}
if ((out==2)&&(sflag==0)) {
printf("parameter error \n");
goto zskip;
}
//
// compute Mertens function
//
m[0]=1;
if (out==0) {
// printf(" %d \n",1);
fprintf(Outfp," %d \n",1);
}
for (index=2; index<=max; index++) {
sum=0;
for (i=2; i<=(index/3); i++)
sum=sum+m[index/i-1];
sum=sum+(index+1)/2;
t=1-sum;
if ((t>32767)||(t<-32768)) {
printf("overflow \n");
goto zskip;
}
m[index-1]=(short)t;
if (out==0) {
//    printf(" %d \n",1-sum);
fprintf(Outfp," %d \n",1-sum);
}
}
if (out==0)
return;
//
// compute sum or histogram
//
if (out==2) {
for (i=0; i<1001; i++)
hisdif[i]=0;
for (i=0; i<10; i++)
maxsav[i]=0.0;
numone=0;
oldsum=0.0;
maxdel=0.0;
savsum=0;
savind=0;
width=1;
}
for (index=2; index<=max; index++) {
newsum=0.0;
for (i=n; i<=index; i+=n) {
if ((index/i)<=limit)
break;
temp=(double)m[index/i-1];
if (temp<0.0)
temp=-temp;
if (temp>0.0) {
if (flag1==0)
newsum=newsum+1.0;
else
newsum=newsum+(double)i;
}
}
if (sflag!=0)
newsum=sqrt(newsum);
if (out==1) {
printf(" %d %e \n",index,newsum);
fprintf(Outfp," %e\n",newsum);
}
else {
delta=newsum-oldsum;
if (delta<0.0) {
delta=-delta;
if ((flag1!=0)&&(limit==0)) {
printf("negative step: %d \n",index);
goto zskip;
}
}
if (delta>1.0)
numone=numone+1;
if (delta>maxdel) {
for (j=9; j>0; j--)
maxsav[j]=maxsav[j-1];
maxsav[0]=delta;
maxdel=delta;
savsum=sum;
savind=index;
}
oldsum=newsum;
if (delta>0.00001) {
width=width+1;
if (width!=(width/n)*n) {
printf("error: width=%d \n",width);
goto zskip;
}
width=0;
j=(unsigned int)(delta*(double)bins+0.5);
if ((limit==0)&&(flag1!=0)) {
if (j<=(bins+100))
hisdif[j]=hisdif[j]+1;
}
else {
if (j<=bins)
hisdif[j]=hisdif[j]+1;
}
}
else
width=width+1;
}
}
if (out==2) {
printf("maximum difference=%e, sum=%d, index=%d \n",maxdel,savsum,savind);
printf("number greater than one=%d \n",numone);
for (i=0; i<10; i++)
printf(" %e \n",maxsav[i]);
j=0;
if (limit==0) {
for (i=0; i<=(bins+100); i++) {
fprintf(Outfp," %d \n",hisdif[i]);
j=j+hisdif[i];
}
}
else {
for (i=0; i<=bins; i++) {
fprintf(Outfp," %d \n",hisdif[i]);
j=j+hisdif[i];
}
}
printf("sum=%d \n",j);
}
zskip:
return;
}
```