maximum.c File Reference

Functions
double	CDF (double x)
double	inverse (double r)
int	main (int argc, char **argv)

Function Documentation

CDF()

double CDF

(

double

x

)

                     {
  return sqrt(3)/(4*M_PI) *( log((x+1)*(x+1)/(x*x-x+1)) +2*sqrt(3) *(atan( (2*x-1)/sqrt(3.0) ) +M_PI/6.0) );
}

inverse()

double inverse

(

double

r

)

                         {
  double min=0, max=0, mid;
  double delta  = 0.1;
  double epsilon=1.0e-6;
  while(CDF(max) < r) max += delta;
  while( fabs(max-min) > epsilon ) {
    mid = (max+min)/2.0;
    if(CDF(mid) >  r) max = mid;
    else min = mid;
  }
  return mid;
}

main()

int main	(	int	argc,
		char **	argv
	)

                               {
  if (argc < 3) {
    printf("Usage: 1.This 2.n 3.ndata\n");
    exit(EXIT_FAILURE);
  }
  Integer_t i, j, count;
  Integer_t const     n = atoi(argv[1]);
  Integer_t const ndata = atoi(argv[2]);
  Integer_t seed=0;
  double max, maxData;
  double rand, rand1, rand2;
 
  Integer_t subMax, itemp;
  Integer_t *histogram;
  double mean, stddev;
  double binwidth, binratio=0.01;
  Integer_t const MAX=1000;
 
  char str[100];
  char *filename = str;
 
  FILE *Infp;
  FILE *Outfp;
  sprintf(str, "Dist_maximum_%lld_raw.txt", n);
  Outfp = fopen(filename, "w");
  if (Outfp == NULL) {
    printf("Can't open a file %s\n", filename);
    exit(EXIT_FAILURE);
  }
  fprintf(Outfp, "#(n=%lld, seed=%lld) Distribution of maximum of |x| \n\n", n, seed);
 
  init_genrand(seed);
  maxData = NAN;
  mean    = 0;
  stddev  = 0;
  for(i = 0;i < ndata; ++i) {
    max = NAN;
    for(j = 0;j < n; ++j) {
      rand1 = genrand_real3();
      rand  = inverse(rand1);
      max = fmax(max, rand);
      // max = rand;
    }
    subMax  = (Integer_t)(max/binwidth) +1;
    if(subMax > MAX) {
      i -=1;
      continue;
    }
    maxData = fmax(maxData, max);
    mean   += max;
    stddev += max*max;
    fprintf(Outfp, "%lld %le\n", i, max);
  }
  mean   /= (double)ndata;
  stddev /= (double)ndata;
  stddev  = sqrt( stddev -mean*mean );
  binwidth *= stddev;
  // %binwidth = 1.0;
  binwidth = 0.01;
  subMax  = (Integer_t)(maxData/binwidth) +1;
  if(subMax > MAX) subMax = MAX;
  subMax=MAX;
  // printf("n=%lld, subMax=%lld, MAX=%lld, maxData=%le, binwidth=%lf\n", n, subMax, MAX, maxData, binwidth);
  histogram = alloc_ivector(subMax); for(i = 0;i < subMax; ++i) histogram[i] = 0;
 
  sprintf(str, "Dist_maximum_%lld_raw.txt", n);
  Infp = fopen(filename, "r");
  if (Infp == NULL) {
    printf("Can't open a file %s\n", filename);
    exit(EXIT_FAILURE);
  }
  skip_header(Infp);
 
  count=0;
  while ( fscanf(Infp, "%lld %le", &i, &max) != EOF ) {
    itemp = (Integer_t)(max/binwidth);
    if(itemp >= subMax) continue;
    if(itemp < 0) continue;
    histogram[itemp] += 1;
    count += 1;
  }
 
  sprintf(str, "Dist_maximum_%lld.txt", n);
  Outfp = fopen(filename, "w");
  if (Outfp == NULL) {
    printf("Can't open a file %s\n", filename);
    exit(EXIT_FAILURE);
  }
  fprintf(Outfp, "# mean=%le, stddev=%le\n\n", mean, stddev);
  fprintf(Outfp, "# 1.Maximum value 2.Probability density\n\n");
  for(i = 0;i < subMax; ++i) {
    fprintf(Outfp, "%le %le\n", binwidth*(i+0.5), histogram[i]/(binwidth*count));
  }
  printf("atan=%lf, %lf\n", atan( -1/sqrt(3.0) ), -M_PI/6);
 
  return 0;
}