/*
Time-stamp: <2009-11-17 19:08:11 nakata>
逆正弦法則
単純ランダムウォークの正の滞在時間の割合
*/

#include <stdio.h>
#include <math.h>
#include <time.h>
#include <stdlib.h>
#define NUM_TRIAL 100000 // パスの試行回数
#define PATH_LENGTH 3000 // パスの長さ
#define NUM_CLASS 100 // 階級の数 
void srand(unsigned int);

int main(void)
{
  double rnd;
  double p;

  int i, Sn, U, current, trial, num_classes, freq_U[NUM_CLASS+1]={0};
  double fdataU[NUM_TRIAL+1]={0}, haba, katan, jyotan;

  srand((unsigned)time(NULL));

  p = 0.5;  // ランダムウォークの推移確率

  for(trial=0; trial<NUM_TRIAL; trial++){
    Sn = 0; // パス毎にランダムウォークを初期化
    U = 0; // パス毎に正の滞在時間を初期化
    for(i=0;i<PATH_LENGTH;i++){
      rnd = rand()/(1.0+RAND_MAX); 
      current = Sn; // 現在の位置
      if(rnd > p) Sn = Sn - 1;
      else Sn = Sn + 1;
      if(current + Sn > 0)
	U = U + 1;
    }
    fdataU[trial] = (double)U/PATH_LENGTH;
    //    printf("(%d) %f\n",trial, fdataU[trial]);
  }

  // 度数分布を求める
  // [0,1] に分布している fdataU[trial] (TRIAL=0〜NUM_TRIAL-1) に関して
  // num_classes の数の階級を作り 
  // 1/num_classes 刻みで度数分布を配列に入れる
  haba = 1.0/NUM_CLASS; // 階級の幅
  for(trial=0;trial < NUM_TRIAL;trial++){
    if(fdataU[trial] <= haba) {
      freq_U[0]++;
    }
    else{
      katan = 0.0;
      for(i=1; i< NUM_CLASS; i++){
	katan = katan + haba;
	jyotan = katan + haba;
	if(katan < fdataU[trial] && fdataU[trial] <= jyotan){
	  freq_U[i]++;
	  break;
	}
      }
    }
  }
  

  //表示

  printf("############\n#Arcsin Law\n");
  printf("#TRIAL %d\n#PATH LENGTH %d\n############\n", NUM_TRIAL, PATH_LENGTH);
  for(i=0; i< NUM_CLASS; i++){
    //    printf("%f\n", (double)freq_U[i]/NUM_TRIAL);
    printf("%f %f\n", (double)i/NUM_CLASS,(double)freq_U[i]*NUM_CLASS/NUM_TRIAL);
  }
  printf("## % a.exe > arcsin.dat \n");
  printf("#####wgnuplot \n");
  printf("## gnuplot> plot  1/(pi*sqrt(x*(1-x))), \"./arcsin.dat\" \n");
  return 0;
}