﻿ generate polygon using Farey series
/*CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                             C
C  CHECK WHETHER LENGTHS, CUMULATIVE LENGTHS, AND SQRT(LENGTH) U.D. (MOD 1)   C
C  04/20/14 (DKC)							      C
C                                                                             C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC*/
#include <stdio.h>
#include <math.h>
unsigned int lagrange(unsigned int N, unsigned int *R);
//
unsigned int sqrtflg=1;  // compute sum of square roots of lengths
unsigned int flag=5;  // 0 (linear), 1 (quadratic), 2 (cubic), 3 (quartic)
//		      // 4 (sine), 5 (cosine)
void main() {
unsigned int N,S[250000],MAXN;
unsigned int L,I;
double temp1,temp2,func,len,tmpsum,expect,lensum,lenfunc,sqrtsum;
FILE *Outfp;
Outfp = fopen("out1c.dat","w");
//
// ORDER OF FAREY SERIES
//
MAXN=600;
/***************/
/* BEGIN LOOP  */
/***************/
for (N=2; N<=MAXN; N++) {
//
// GENERATE FAREY SERIES (USING HAROS' THEOREM AND LAGRANGE'S ALGORITHM)
//
L=lagrange(N,S);
len=0.0;
tmpsum=0.0;
lensum=0.0;
sqrtsum=0.0;
for (I=0; I<(L-1); I++) {
temp1=((double)S[2*I]-(double)S[2*I+2])*((double)S[2*I]-(double)S[2*I+2]);
temp1=temp1+((double)S[2*I+1]-(double)S[2*I+3])*((double)S[2*I+1]-(double)S[2*I+3]);
len=len+sqrt(temp1);
//
// fractional part of length
//
if (sqrtflg==0)
temp1=sqrt(temp1);
else
temp1=sqrt(sqrt(temp1));
sqrtsum=sqrtsum+temp1;
temp1=temp1-(unsigned int)temp1;
temp2=len-(unsigned int)len;
if (flag==0) {
func=temp1;
lenfunc=temp2;
expect=0.5;
}
if (flag==1) {
func=temp1*temp1;
lenfunc=temp2*temp2;
expect=0.3333333;
}
if (flag==2) {
func=temp1*temp1*temp1;
lenfunc=temp2*temp2*temp2;
expect=0.25;
}
if (flag==3) {
func=temp1*temp1*temp1*temp1;
lenfunc=temp2*temp2*temp2*temp2;
expect=0.20;
}
if (flag==4) {
func=sin(temp1);
lenfunc=sin(temp2);
expect=cos(1.0);
}
if (flag==5) {
func=cos(temp1);
lenfunc=cos(temp2);
expect=sin(1.0);
}
tmpsum=tmpsum+func;
lensum=lensum+lenfunc;
}
len=(len+1.0);
temp2=len-(unsigned int)len;
len=len/(double)L;
if (flag==0)
lenfunc=temp2;
if (flag==1)
lenfunc=temp2*temp2;
if (flag==2)
lenfunc=temp2*temp2*temp2;
if (flag==3)
lenfunc=temp2*temp2*temp2*temp2;
if (flag==4)
lenfunc=sin(temp2);
if (flag==5)
lenfunc=cos(temp2);
lensum=lensum+lenfunc;
lensum=lensum/(double)L;
//
if (flag==5) 	     // cos(0)=1
tmpsum=tmpsum+1.0;
tmpsum=tmpsum/(double)L;
sqrtsum=(sqrtsum+1.0)/(double)L;
printf("%d %d %e %e %e %e %e \n",N,L,len,sqrtsum,lensum,tmpsum,expect);
if (sqrtflg==0)
fprintf(Outfp," %e\n",len);
else
fprintf(Outfp," %e\n",sqrtsum);
}
/***************/
/* END LOOP    */
/***************/
fclose(Outfp);
return;
}