/*CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C C
C INTERPOLATE FRACTIONS AND COMPUTE FUNCTIONS C
C 05/17/15 (DKC) C
C C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC*/
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include "table2.h"
#include "output2.h" // if output5, scale=1000, MAXN=9999
// if output4, scale=1000, MAXN=3999
// if output2, scale=1000000, MAXN=999
// if output1, scale=1000000, MAXN=468 (computed on C64)
double numdiv(unsigned int a, unsigned int flag);
int liouvill(unsigned int a, unsigned int *t);
double mangoldt(unsigned int a, unsigned int *t);
double dirchar(unsigned int a, unsigned int p);
double interp1(double ND, double ND1, unsigned int I,
double init, double *out);
double scale=1000000.0; // for output1.h or output2.h
//double scale=1000.0; // for output4.h or output5.h
unsigned int count=1000; // number of limits
unsigned int MAXN=999; // maximum x
unsigned int delta=0; // if set, add offset (slope for out=0)
unsigned int weight=0; // multiply by h/i
unsigned int prime=101; // prime for Legendre symbol
unsigned int I=1; // m value
unsigned int out=12; // 0 for sum of delta(x/i)
// 1 for sum of delta(x/i)*log(i)
// 2 for sum of delta(x/i)*sigma0(i)
// 3 for sum of delta(x/i)*sigma1(i)
// 4 for sum of delta(x/i)*sigma2(i)
// 5 for sum of delta(x/i)^2
// 6 for sum of delta(x/i)*M(i)
// 7 for sum of delta(x/i)*L(i)^2
// 8 for sum of delta(x/i)^2*L(i)^2
// 9 for sum of delta(x/i)^2*L(i)
// 10 for sum of delta(x/i)*log(i)*d(i)
// 11 for sum of delta(x/i)/i
// 12 for sum of delta(x/i)*mangoldt(i)
// 13 for sum delta(x/i), i perfect square
// 14 for sum delta(x/i)*legendre(i)
// 99 for interpolated limits
double slope[8]={.1704,.1548,.1402,.128,.1179,.1095,.1023,.09617};
//
void main() {
double init,sum11,temp,tsum;
double *output;
unsigned int N,h,i,J,tmp;
int t,sum;
short M[100001];
FILE *Outfp;
Outfp = fopen("out6.dat","w");
if ((count*I)>40000000) {
printf(" array not big enough \n");
goto zskip;
}
if (MAXN>40000000) {
printf(" MAXN too big \n");
goto zskip;
}
if (MAXN>(count*I)) {
printf(" MAXN too big \n");
goto zskip;
}
if ((MAXN>199999)&&(out==12)) {
printf(" MAXN too big \n");
goto zskip;
}
output=(double*) malloc(40000000);
if (output==NULL)
return;
//
// compute Mertens function
//
if (out==6) {
if (MAXN>100000) {
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;
M[N-1]=(short)t;
// printf(" %d \n",M[N-1]);
}
}
for (i=0; i<100000; i++)
output[i]=0.0;
if (I<9)
temp=slope[I-1];
else
temp=0.0;
J=(I+1)/2;
init=0.0;
for (i=0; i<count; i++)
init=interp1((double)frac[i]/scale,(double)frac[i+1]/scale,I,init,
&output[i*I]);
if (out==99) {
for (i=0; i<(count*I); i++) {
printf(" %e \n",output[i]);
fprintf(Outfp," %e \n",output[i]);
}
goto zskip;
}
for (h=2; h<=MAXN; h++) {
sum11=0.0;
tsum=0;
for (i=1; i<=h; i++) {
if (out==0) {
if (weight==0)
sum11=sum11+output[h/i-1];
else
sum11=sum11+output[h/i-1]*(double)(h/i);
}
if (out==1) {
if (delta==0)
sum11=sum11+output[h/i-1]*log((double)i);
else
sum11=sum11+(output[h/i-1]+temp)*log((double)i);
}
if (out==2) {
if (delta==0)
sum11=sum11+output[h/i-1]*2.0*numdiv(i,1);
else
sum11=sum11+(output[h/i-1]+temp)*2.0*numdiv(i,1);
}
if (out==3) {
if (delta==0)
sum11=sum11+output[h/i-1]*numdiv(i,4);
else
sum11=sum11+(output[h/i-1]+temp)*numdiv(i,4);
}
if (out==4) {
if (delta==0)
sum11=sum11+output[h/i-1]*numdiv(i,5);
else
sum11=sum11+(output[h/i-1]+temp)*numdiv(i,5);
}
if (out==5) {
if (delta==0)
sum11=sum11+output[h/i-1]*output[h/i-1];
else
sum11=sum11+(output[h/i-1]+temp)*(output[h/i-1]+temp);
}
if (out==6) {
if (delta==0)
sum11=sum11+output[h/i-1]*(double)M[i-1];
else {
if (weight==0)
sum11=sum11+(output[h/i-1]+temp)*(double)M[i-1];
else
sum11=sum11+(output[h/i-1]+temp)*(double)M[i-1]*(double)(h/i);
}
}
if (out==7) {
t=liouvill(i, table);
tsum=tsum+t;
if (delta==0)
sum11=sum11+output[h/i-1]*(double)tsum*(double)tsum;
else
sum11=sum11+(output[h/i-1]+temp)*(double)tsum*(double)tsum;
}
if (out==8) {
t=liouvill(i, table);
tsum=tsum+t;
if (delta==0)
sum11=sum11+output[h/i-1]*output[h/i-1]*(double)tsum*(double)tsum;
else
sum11=sum11+(output[h/i-1]+temp)*(output[h/i-1]+temp)*(double)tsum*(double)tsum;
}
if (out==9) {
t=liouvill(i, table);
tsum=tsum+t;
if (delta==0)
sum11=sum11+output[h/i-1]*output[h/i-1]*(double)tsum;
else
sum11=sum11+(output[h/i-1]+temp)*(output[h/i-1]+temp)*(double)tsum;
}
if (out==10) {
if (delta==0)
sum11=sum11+output[h/i-1]*log((double)i)*numdiv(i,1);
else
sum11=sum11+(output[h/i-1]+temp)*log((double)i)*numdiv(i,1);
}
if (out==11) {
if (delta==0)
sum11=sum11+output[h/i-1]/(double)i;
else
sum11=sum11+(output[h/i-1]+temp)/(double)i;
}
if (out==12) {
if (delta==0) {
if (weight==0)
sum11=sum11+output[h/i-1]*mangoldt(i,table);
else
sum11=sum11+output[h/i-1]*mangoldt(i,table)*(double)(h/i);
}
else
sum11=sum11+(output[h/i-1]+temp)*mangoldt(i,table);
}
if (out==13) {
tmp=(unsigned int)(sqrt((double)i)+0.001);
if ((tmp*tmp)==i) {
if (delta==0)
sum11=sum11+output[h/i-1];
else {
if (weight==0)
sum11=sum11+output[h/i-1]+temp;
else
sum11=sum11+(output[h/i-1]+temp)*(double)(h/i);
}
}
}
if (out==14) {
// printf(" %d %d %e \n",i,output[h/i-1],dirchar(i,prime));
if (delta==0)
sum11=sum11+output[h/i-1]*dirchar(i,prime);
else
sum11=sum11+(output[h/i-1]+temp)*dirchar(i,prime);
}
}
if (out!=5) {
printf(" %d %e \n",h,-sum11);
fprintf(Outfp," %e\n",-sum11);
}
else {
printf(" %d %e \n",h,sum11);
fprintf(Outfp," %e\n",sum11);
}
}
zskip:
free(output);
fclose(Outfp);
return;
}