total_iterations += loop;  return ret;}double R(double a, double b)				// function for check node operation{  return 2.0 * atanh(tanh(a/2.0)*tanh(b/2.0));}double sign(double a){  if (a == 0) return 0.0;  if (a > 0) return 1.0;  if (a < 0) return -1.0;}double min(double a, double b){  if (a > b) return b;  else return a;}double minsumR(double a, double b)				// function for check node operation of min-sum algotithm{  return sign(a) * sign(b) * min(fabs(a),fabs(b));}void init_dec(FILE *fp){/**********************************************************************************/// initialization of BP decoder/**********************************************************************************/  int i,j;  int v;  int *counter;  double a,b;  int tmp;// read parity check matrix   fscanf(fp,"%d %d\n",&n,&m);	// n:code length, m:number of parity checks  fscanf(fp,"%d %d\n",&rmax,&cmax); // rmax:maximum number of row weight				    // cmax:maximum number of column weight  row_weight = malloc(sizeof(int)*m); // array for row weight profile  for (i = 0; i <= m-1; i++) fscanf(fp,"%d",&row_weight[i]);  col_weight = malloc(sizeof(int)*n); // array for column weight profile  for (i = 0; i <= n-1; i++) fscanf(fp,"%d",&col_weight[i]);  counter = malloc(sizeof(int)*n);   for (i=0; i <= n-1; i++) counter[i] = 0;    row_list = malloc(sizeof(int*)*m); // parity check matrix in row form  for (i=0;i<=m-1;i++) row_list[i] = malloc(sizeof(int)*rmax);  col_list_r = malloc(sizeof(int*)*n);   for (i=0;i<=n-1;i++) col_list_r[i] = malloc(sizeof(int)*cmax);  col_list_c = malloc(sizeof(int*)*n);   for (i=0;i<=n-1;i++) col_list_c[i] = malloc(sizeof(int)*cmax);// set up for parity check matrix  for (j = 0; j <= m-1; j++) {    for (i = 0; i <= row_weight[j]-1; i++) {      fscanf(fp,"%d",&v);      v--;      row_list[j][i] = v;      col_list_r[v][counter[v]] = j;      col_list_c[v][counter[v]] = i;      counter[v]++;    }  }  punc_ptn = malloc(sizeof(int) * n); // puncture pattern (0:no puncture, 1:puncture)  if (puncture == ON)    for (i = 0; i<=n-1; i++) fscanf(fp,"%d",&punc_ptn[i]);  else     for (i = 0; i<=n-1; i++) punc_ptn[i] = 0;  if (met == ON) {    vnode_type = malloc(sizeof(int)*n);    fscanf(fp,"%d",&tmp);    for (i = 0; i<=m-1; i++) fscanf(fp,"%d",&tmp);    fscanf(fp,"%d",&num_vnode_type);        for (i = 0; i<=n-1; i++) fscanf(fp,"%d",&vnode_type[i]);    vnode_errors = malloc(sizeof(int)*num_vnode_type);    for (i = 0; i <= num_vnode_type-1; i++) vnode_errors[i] = 0;  }  // set up for BP related global variables  alpha = malloc(sizeof(int*)*m); // check to variable message  for (i=0;i<=m-1;i++) alpha[i] = malloc(sizeof(int)*rmax);  beta  = malloc(sizeof(int*)*m); // variable to check message  for (i=0;i<=m-1;i++) beta[i] = malloc(sizeof(int)*rmax);  tmp_bit = malloc(sizeof(int)*n); // tentative decoding result				// set up of table for check node operation  R_tbl = malloc(sizeof(int*) * ASIZE);  for (i = 0; i <= ASIZE-1; i++) {    R_tbl[i] = malloc(sizeof(int) * ASIZE);  }  for (i = LB2; i <= UB2; i++) {    for (j = LB2; j <= UB2; j++) {      a = DELTA * i;      b = DELTA * j;      if (minsum == 1) 	R_tbl[i+UB][j+UB] = quantize(minsumR(a,b));      else	R_tbl[i+UB][j+UB] = quantize(R(a,b));    }  }  qLLR = malloc(sizeof(int)*n);	// quantized LLR  forward = malloc(sizeof(int)*rmax);  backward= malloc(sizeof(int)*rmax);}/**********************************************************************************/// simulation status and result report functions/**********************************************************************************/void report_result2(FILE *out){  fprintf(out,"%4.3f %7.6e %7.6e %d %d %d %d %5.2f %d %d\n",	  snr,	  (double)error_bits/total_bits,	  (double)error_blocks/total_blocks,	   error_bits,	  total_bits,	  error_blocks,	  total_blocks,	  (double)total_iterations/total_blocks,	  n,	  m	  );}void report_result(FILE *out){  fprintf(out,"%4.3f %7.6e %7.6e %d %d %d %d %5.2f var=%4.3f sdv=%4.3f mis=%d seed=%d mitr=%d n=%d m=%d r=%4.3f depth=%d DLETA=%4.3e range=%4.3e config=%s scond=%d serr=%d target=%4.3e log=%d punc=%d met=%d msum=%d norm=%4.3f pchk=%d\n",	  snr,	  (double)error_bits/total_bits,	  (double)error_blocks/total_blocks,	  error_bits,	  total_bits,	  error_blocks,	  total_blocks,	  (double)total_iterations/total_blocks,	  channel_var,	  sqrt(channel_var),	  total_miscorrection,	  seed,	  max_itr,	  n,	  m,	  code_rate,	  quantization_level,	  DELTA,	  UB2*DELTA,	  filename,	  stop_condition,	  stop_errors,	  target_error_prob,	  error_logging,	  puncture,	  met,	  minsum,	  norm_factor,	  parity_check	  );}/**********************************************************************************/// simulation functions/**********************************************************************************/void make_received_word(double *r){  int i;  for (i = 0; i <= n-1; i++) {    r[i] = 1.0-2.0*cword[i] + nrnd(channel_var);  }}void simulation_loop(){  int i,j;  int sum;  int ret;// set up of initial condition  cword = malloc(sizeof(int)*n);  rword = malloc(sizeof(double)*n);  LLR = malloc(sizeof(double)*n);  total_blocks = 0;  total_bits = 0;  error_blocks = 0;  error_bits = 0;  total_iterations = 0;  total_miscorrection = 0;  while(1) {    total_blocks++;    total_bits += n;    if (encoding == 1) {      gen_info_seq(cword);      encode_word(cword);    }     else for (i = 0; i <= n-1;i++) cword[i] = 0;    //for (i =0; i <= n-1; i++) printf("%d",cword[i]);printf("\n");    //printf("parity = %d\n",parity_check_func(cword));    make_received_word(rword);    for (i = 0; i <= n-1; i++) LLR[i] = 2.0 * rword[i]/channel_var;    ret = quantized_BP(LLR);    sum = 0;    for(i = 0; i <= n-1;i++) {      if (encoding == 1) tmp_bit[i] ^= cword[i];      if (punc_ptn[i] == 0) sum += tmp_bit[i];      if (met == ON) 	vnode_errors[vnode_type[i]] += tmp_bit[i];    }    error_bits += sum;    if (sum != 0) {      error_blocks++;      if (error_logging==1) {	for (j=0; j<=n-1;j++) fprintf(error_log,"%d ",tmp_bit[j]);	fprintf(error_log,"(error_weight=%d)\n",sum);      }    }    if ((ret == 0)&&(sum != 0)) total_miscorrection++;    if (display_result == 1) report_result(stderr); // full report    if (display_result == 2) report_result2(stderr); // brief report    if ((error_bits >= stop_errors)&&(stop_condition==BIT)&&(total_blocks>=100)) break;    if ((error_blocks >= stop_errors)&&(stop_condition==BLOCK)) break;    if ((total_bits >= stop_errors/target_error_prob)&&(stop_condition==BIT)) {      fprintf(stdout, "stop by target_error_prob conditon\n");      break;    }    if ((total_blocks >= stop_errors/target_error_prob)&&(stop_condition==BLOCK)) {      fprintf(stdout, "stop by target_error_prob conditon\n");      break;    }    if ((total_bits >= 10/target_error_prob)&&(error_bits==0)&&(stop_condition==BIT)) {      fprintf(stdout, "stop by early stop condition\n");      break;    }    if ((total_blocks >= 10/target_error_prob)&&(error_blocks==0)&&(stop_condition==BLOCK)) {      fprintf(stdout, "stop by early stop condition\n");      break;    }    if ((stop_condition >=2)&&(total_blocks >stop_condition)) break;  }}void terminate_program(char *s){    fprintf(stderr,"%s\n",s);    exit(-1);}void read_config_file(FILE *fp){// set up for simulation parameters  FILE *hfile,*encoder;  char str1[MAX_STR_LEN],str2[MAX_STR_LEN];  int i;  fscanf(fp,"%s",str1);  if (strcmp(str1,"matrix_file")!=0) terminate_program("error in matrix_file\n");  fscanf(fp,"%s",str2);		// parity check matrix file (spmat form)  if ((hfile = fopen(str2,"r")) == NULL) {    fprintf(stderr,"Can't open %s.\n",str2);    exit(-1);  }  fscanf(fp,"%s",str1);  if (strcmp(str1,"code_rate")!=0) terminate_program("error in code_rate\n");  fscanf(fp,"%lf",&code_rate);	// code rate of the target LDPC code  if ((code_rate<0) || (code_rate>1)) terminate_program("invalid code rate\n");  fscanf(fp,"%s",str1);  if (strcmp(str1,"quantize_depth")!=0) terminate_program("error in quantize_depth\n");  fscanf(fp,"%d",&quantization_level); // quantization depth  if ((quantization_level<2) || (quantization_level>20))     terminate_program("invalid quantization_level\n");  ASIZE = (1 << quantization_level)-1; // number of quantization levels  DELTA = 102.4/(ASIZE+1);	// quatization step  LB = -(ASIZE-1)/2;		  UB = (ASIZE-1)/2;  LB2 = -(ASIZE-1)/4;		// lower bound of quantized value  UB2 = (ASIZE-1)/4;		// upper bound of quantized value  fscanf(fp,"%s",str1);  if (strcmp(str1,"max_iteration")!=0) terminate_program("error in max_iteration\n");  fscanf(fp,"%d",&max_itr);	// maximum number of iterations in BP  if (max_itr < 0) terminate_program("invalid max_itr\n");  fscanf(fp,"%s",str1);  if (strcmp(str1,"random_seed")!=0) terminate_program("error in random_seed\n");  fscanf(fp,"%d",&seed);	// seed of random number generator  srand48(seed);		// set up random seed  fscanf(fp,"%s",str1);  if (strcmp(str1,"display")!=0) terminate_program("error in display\n");  fscanf(fp,"%d",&display_result); // flag for display mode of simulation status				   // 0:no display, 1:full display, 2:simple display  fscanf(fp,"%s",str1);  if (strcmp(str1,"stop_condition")!=0) terminate_program("error in stop_condition\n");  fscanf(fp,"%d",&stop_condition); // stop condition of simulation				   // 0:bit errors, 1:block errors, 				   // other positive number: number of loops  if (stop_condition < 0) terminate_program("invalid stop_condition\n");  fscanf(fp,"%s",str1);  if (strcmp(str1,"stop_errors")!=0) terminate_program("error in stop_errors\n");  fscanf(fp,"%d",&stop_errors);	// When the number of errors reaches stop_errors,				// simulation loop is terminated.  if (stop_errors < 0)     terminate_program("invalid stop_errors\n");  fscanf(fp,"%s",str1);  if (strcmp(str1,"target_prob")!=0) terminate_program("error in target_prob\n");  fscanf(fp,"%lf",&target_error_prob); // target error probability  if (target_error_prob < 0.0)     terminate_program("invalid target_error_prob\n");  fscanf(fp,"%s",str1);  if (strcmp(str1,"error_logging")!=0) terminate_program("error in error_loging\n");  fscanf(fp,"%d",&error_logging); // error logging switch 0 or 1  if ((error_logging!=0)&&(error_logging!=1))     terminate_program("invalid error_logging\n");  if(error_logging == 1) {    if ((error_log = fopen("error_log","w")) == NULL) {      fprintf(stderr,"Can't open %s.\n","error_log");      exit(-1);    }    fprintf(stderr,"error logging starts.\n");  }  fscanf(fp,"%s",str1);  if (strcmp(str1,"puncture")!=0) terminate_program("error in puncture\n");  fscanf(fp,"%d",&puncture);	// Are there puncture bits? 0 or 1  if ((puncture!=0)&&(puncture!=1))     terminate_program("invalid puncture\n");  fscanf(fp,"%s",str1);  if (strcmp(str1,"multi_edge_type")!=0) terminate_program("error in multi_edge_type\n");  fscanf(fp,"%d",&met);		// Is it a multi-edge type LDPC code? 0 or 1  if ((met!=0)&&(met!=1))     terminate_program("invalid met\n");  fscanf(fp,"%s",str1);  if (strcmp(str1,"min-sum")!=0) terminate_program("error in min-sum\n");  fscanf(fp,"%d",&minsum);	// switch of min-sum algorithm 0 or 1				// if minsum == 1, then simulation with min-sum algorithm				// otherwise simulation with sum-product algortihm  if ((minsum!=0)&&(minsum!=1))     terminate_program("invalid minsum\n");  fscanf(fp,"%s",str1);  if (strcmp(str1,"norm_factor")!=0) terminate_program("error in norm_factor\n");  fscanf(fp,"%lf",&norm_factor); // normalization factor for min-sum algorithm				 // if normalization == 0, then no normalization				 // See papers of normalized-BP by Prof.Fossorier.  if (norm_factor < 0.0)     terminate_program("invalid norm_factor\n");  fscanf(fp,"%s",str1);  if (strcmp(str1,"parity_check")!=0) terminate_program("error in parity_check\n");  fscanf(fp,"%d",&parity_check); // switch of parity check in BP algortihm, 0 or 1  if ((parity_check!=0)&&(parity_check!=1))     terminate_program("invalid parity_check\n");  #ifdef TEST2  encoding = 1;  printf("encoding is ON\n");  fscanf(fp,"%s",str1);  if ((encoder = fopen(str1,"r")) == NULL) {    fprintf(stderr,"Can't open %s.\n",str1);    exit(-1);  }  read_encfile(encoder);  //print_encfile(encoder);  #endif  init_dec(hfile);		// initialization of decoder}/**********************************************************************************/// main /**********************************************************************************/int main(int argc,char **argv){  FILE *fp;			// pointer for configuration file  int i;#ifdef TEST1/*Test for quantized_BP()-----------------------------------------------------------------6.3.spmat------------------6 33 23 2 31 1 2 2 1 11 2 3 3 44 5 61 0 0 0 0 1------------------log posterior probability ratio (by double-precision BP)    4.0554    0.5884    0.2407    1.0694    1.6060    2.2915    3.9537    0.2780    1.7111    1.7111    1.5387    2.2338    4.3974    1.6925    1.7111    1.7111    2.0840    2.7033    4.3974    1.6925    1.7111    1.7111    2.0840    2.7033    4.3974    1.6925    1.7111    1.7111    2.0840    2.7033*/  double r[] = {1.620803,0.264281,-0.031637,-0.127654,0.746347,1.003543};   int qlevel;  channel_var = 0.794328;  printf("basic test\n");  fp = fopen("6.3.spmat","r");  qlevel = 10;  ASIZE = (1 << qlevel)-1;   DELTA = 102.4/(ASIZE+1);  LB = -(ASIZE-1)/2;		  UB = (ASIZE-1)/2;  LB2 = -(ASIZE-1)/4;		  UB2 = (ASIZE-1)/4;		  puncture = OFF;  met = OFF;  parity_check = OFF;  max_itr = 10;  init_dec(fp);  LLR = malloc(sizeof(int)*6);  for (i = 0;i <= 5; i++) LLR[i] = 2.0*r[i]/channel_var;  quantized_BP(LLR);  exit(0);//-----------------------------------------------------------------#endif// user interface  if (argc != 3) {    printf("usage: awgn_simulator config_file snr\n");				// snr: Eb/N0    exit(-1);  }  if ((fp = fopen(argv[1],"r")) == NULL) {    fprintf(stderr,"Can't open %s.\n",argv[1]);    exit(-1);  }  filename =argv[1];		// filename of configuration file  snr = atof(argv[2]);		// get signal to noise ratio  read_config_file(fp);		// reading configulation file  channel_var = 0.5 * (1.0/pow(10.0,snr/10.0))/code_rate;// simulation and report  simulation_loop();		// main simulation loop  report_result(stdout);	// print out the simulation result  if (error_logging==ON) report_result(error_log);				// for multi-edge type  if (met == ON) {    printf("# number of bit errors for each variable node type\n");    printf("# ");    for (i=0;i<=num_vnode_type-1;i++) printf("%d ",vnode_errors[i]);    printf("\n");    printf("# bit error probability for each variable node type\n");    printf("# ");    for (i=0;i<=num_vnode_type-1;i++) printf("%16.12e ",(double)vnode_errors[i]/total_bits);    printf("\n");  }  return 0;}