/*CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                             C
C  COMPUTE SUM OF EULER'S PHI FUNCTION                                        C
C  12/23/14 (DKC)							      C
C                                                                             C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC*/
#include <stdio.h>
#include <math.h>
unsigned int euclid(unsigned int d, unsigned int e);
//
unsigned int order=50000;
unsigned int out=2;   // losum, hisum, hisum-losum, sqrt((hisum-losum)), sqrt(A(x))
unsigned int flag=0;  // for out=4, sqrt(A(x)), A(x)/sqrt(x(x+1)(2x+1)/6),
		      // or g(x)/sqrt(x(x+1)(2x+1)/6)
unsigned int flag4=0; // scale for out=4 (zero, plus, or minus)
unsigned int sflag=1; // square root flag (for out=4, normally set to 0)
unsigned int flag2=1; // for out=2, if set compute 1.5x/log(x)-(hisum-losum)
//
void main() {
unsigned int g,h,j,k,l,r;
unsigned int sum,losum,hisum,count;
double misum,esum,edelt;
FILE *Outfp;
Outfp = fopen("out1r1.dat","w");
losum=0;
hisum=1;   // changed from 0 on 08/19/14
misum=1.0;
count=0;
for (g=2; g<=order; g++) {
   if (out==4) {
      k=g/6;
      r=g-k*6;
      if (k==0) {
	 if (r==1)
	    esum=1.0;
	 if (r==2)
	    esum=2.0;
	 if (r==3)
	    esum=4.0;
	 if (r==4)
	    esum=6.0;
	 if (r==5)
	    esum=11.0;
	 }
      else {
	 esum=12.0*k+23.0*((double)k*(k-1))/2;
	 edelt=0.0;
	 if (r==1)
	    edelt=7.0+6.0*(k-1);
	 if (r==2)
	    edelt=11.0+9.0*(k-1);
	 if (r==3)
	    edelt=17.0+13.0*(k-1);
	 if (r==4)
	    edelt=22.0+16.0*(k-1);
	 if (r==5)
	    edelt=33.0+22.0*(k-1);
	 esum=esum+edelt;
	 }
      }
   sum=1;
   for (j=2; j<g; j++) {
      k=euclid(g,j);
      if (k==1)
	 sum=sum+1;
      }
   misum=misum+(double)sum;
   l=sum/4;
   h=l;
   if (sum!=l*4)
      h=h+1;
   else
      count=count+1;
   losum=losum+l;
   hisum=hisum+h;
   printf(" %d %d %d %d %e %d \n",g,sum,losum,hisum,misum,(int)(hisum-losum));
   if (out==0)
      fprintf(Outfp," %d\n",losum);
   else {
      if (out==1)
	 fprintf(Outfp," %d\n",hisum);
      else {
	 if (out==2) {
	    if (flag2==0)
	       fprintf(Outfp," %d\n",hisum-losum);
	    else
	       fprintf(Outfp," %e\n",1.5*(double)g/log((double)g)-(hisum-losum));
	    }
	 else {
	    if (out==3)
	       fprintf(Outfp," %e\n",sqrt((double)(hisum-losum)));
	    else {
	       if (flag==0) {
		  if (flag4==0) {
		     if (sflag!=0)
			fprintf(Outfp," %e\n",sqrt(misum));
		     else
			fprintf(Outfp," %e\n",misum);
		     }
		  if (flag4==1)
		     fprintf(Outfp," %e\n",0.5513*(double)g+0.2767+0.5);
		  if (flag4==2)
		     fprintf(Outfp," %e\n",0.5513*(double)g+0.2767-0.5);
		  }
	       else {
		  if (flag==1)
		     fprintf(Outfp," %e\n",misum/sqrt((double)g*(double)(g+1)*(double)(2*g+1)/6.0));
		  else
		     fprintf(Outfp," %e\n",esum/sqrt((double)g*(double)(g+1)*(double)(2*g+1)/6.0));
		  }
	       }
	    }
	 }
      }
   }
printf("count=%d, ratio=%e \n",count,(double)count/(double)order);
fclose(Outfp);
return;
}