cmatrix.c

快速傅立叶变换程序代码,学信号的同学,可要注意了

      for ( j = 1 ; j <= Cols ; j ++ ) 
	pikd_r[j] = 0 ;

      for ( num_left = per ; num_left >= 1 ; num_left -- ) {
	do {
	  un = ranf() ; 
	  j = read_int_from_cum ( cum , 1 , Cols , un ) ; 
	} while ( pikd_r[j] ) ; /* while j already picked */
	pikd_r[j] ++ ; 
	if ( num_r[j] > 0 ) {
	  mylist[num_left] = j ; 
	} else { fprintf (stderr , "cum failure %d %d %d\n" ,
			  j , cum[j-1] , cum[j] ) ; 	  num_left ++ ;
	       pause_for_return(); }
      }
      if ( num_left == 0 ) { /* check hd then move on */
	if ( hd > 0 ) {
	  hdfail = 0 ; 
/* 
   for each j on the list
   go and find any other rows that share this j
   check if we have already met this row today. If so, then that's 
   a failure for the traditional setting of hd
*/
	  rr = 0 ; 
	  prevrr = 0 ; 
	  for ( j = 1 ; j <= per && ( hdfail == 0 ) ; j ++ ) {
	    nn = mylist[j] ;
	    for ( u = num_nlist[nn] ; u >= 1 && ( hdfail == 0 ) ; u-- ) {
	      rr ++ ;
	      latest = relatives[rr] = nlist[nn][u] ; 
	      for ( ss = prevrr ; ss >= 1 ; ss -- ) {
		if ( relatives[ss] == latest ) {
		  hdfail = 1 ; 
		  break ; 
		}
	      }
	    }
	    prevrr = rr ; 
	    if ( rr > relmax ) {
	      fprintf ( stderr, "\nwarning!!, relatives overran vector %d %d\n" , rr , relmax ) ;
	      status -- ; 
	    }
	  }
	}
	if ( hd > 0 && hdfail == 1 ) {
	  hfails ++ ;
	  printf ( "h" ) ; fflush(stdout) ;
	  if ( i == 1 ) { /* tough luck, have to accept the failure */
	    fprintf ( stderr, "\nwarning, final hamming check failure at relative %d %d out of %d\n" , rr , ss , relmax ) ;
	    hdfail = 0 ; 
	  }
	  if ( hfailmax && ( hfails > hfailmax) ) { /* tough luck, have to accept the failure */
	    fprintf ( stderr, "\nwarning, hamming check failure owing to boredom %d %d %d\n" , rr , ss , relmax ) ;
	    hdfail = 0 ; 
	  }
	} 
	if ( hdfail == 0 ) {
	  /* final acts before moving on */
	  for ( j = 1 ; j <= per ; j ++ ) {
	    nn = mylist[j] ;
	    add_to_alist ( a , nn , m ) ; 
	    num_r[nn] -- ; 
	  }
	}
      } else {
	printf ( "." ) ; /* fflush(stdout) ;  */
      }
    } while ( num_left != 0 || hdfail != 0 ) ;
  }
  free_ivector ( mylist , -1 , per ) ;
  free_ivector ( num_r , 1 , Cols ) ;
  free_ivector ( relatives , 1 , relmax ) ;
  
  return status ; 
}

int staircase_alist ( alist_matrix *a , int c0 , int c1 , 
				     int per , int Cols , 
				     int Rows ){
/* 
   Constructs a matrix by row, from bottom to top.
   Each row is a step of a staircase. 
*/
  int i , n , m , u ;
  int num_rows_to_do = c1 - c0 ; 
  int status = 0 ; 
  n = 1 ;
  printf ( "Stair:" ) ;
  for ( i = num_rows_to_do , m = Rows - c0 ; i >= 1 ; i-- , m-- ) {
    if ( !(i%1000) ) { printf ( "%d;" , i ) ; fflush(stdout) ; }

    for ( u = 1 ; u <= per ; u++ ) {
      add_to_alist ( a , n , m ) ; 
      if ( u < per ) { n ++ ; if ( n > Cols ) n = 1 ; }
    }
  }
  return status ;
}

int slope_alist ( alist_matrix *a , int c0 , int c1 , 
				     int per , int Cols , 
				     int Rows ){
/* 
   Constructs a matrix by row, from bottom to top.
   Each row is  
   1 1 ----------
   --- 1 1 ------
   ------- 1 1 -- if per should be 2.
*/
  int i , n , m , u ;
  int num_rows_to_do = c1 - c0 ; 
  int status = 0 ; 
  n = 1 ;
  printf ( "Slope:" ) ;
  for ( i = num_rows_to_do , m = Rows - c0 ; i >= 1 ; i-- , m-- ) {
    if ( !(i%1000) ) { printf ( "%d;" , i ) ; fflush(stdout) ; }

    for ( u = 1 ; u <= per ; u++ ) {
      add_to_alist ( a , n , m ) ; 
      n ++ ; if ( n > Cols ) n = 1 ; 
    }
  }
  return status ;
}

/* used by code5 and code6 */
int uneven_sparse_rectangular_alist ( alist_matrix *a , int c0 , int c1 , 
				     double per , int Cols , 
				     int Rows , int hd , int maxattempts ){
/* 
   Constructs a matrix by row, from bottom to top.
   A start can be made a distance c0 from the bottom if desired 
   A stop can be called after only c1 rows from the bottom have been 
   finished also.  

   if hd > 0 then max permitted overlap is hd DIFFERENT FROM WHAT WENT B4
   if a row is hd or closer to another row, then a "h" is printed
   and we restart that row 
*/
  int verbose = 0 ; 
  int status = 0 ; 
  int weight ; /* used to check if we have room to make our t ones */
  int nn , u , latest ; 
  int i , j , num_left  ;
  int *num_c ; /* number left in col */
  int *pikd_c ;   /* whether this col pikd this time */
  int *cum ;   /* cumulative probability */
  int *mylist ; 
  double per_col ;
  int per_col1 , per_col2 , tot_per_col , tot_ones , hdfail = 0 ;
  int extras_r , extras_c ;
  int m ;
  int num_rows_to_do = c1 - c0 ; 
  int *relatives ; 
  int relmax , rr , ss , prevrr ;
  int *num_nlist = a->num_nlist ; 
  int **nlist = a->nlist ; 
  double un ; /* random 0,1 */
  int attempts  ; /* , maxattempts = 20 - number of times to try to get no-overlap */
  static int Attempts = 0 ;  /* not used: counts clangers */
  int Maxattempts = 20000 ; /* maximum total number of clangers allowed  (not used ) */
  int per1 , per2 , this_per ; 

  /* Compute what numbers per row and column */

  if ( num_rows_to_do == 0 ) { fprintf ( stderr , "0 rows to do\n" ) ;
			       return (status);
			     }
  tot_ones = (int)( per * (double) num_rows_to_do ) ;
  per_col = ( per * (double) num_rows_to_do ) / (double) Cols ;  
  per_col1 = tot_ones / Cols ; 
  per_col2 = per_col1 + 1 ; 
  extras_c = tot_ones - per_col1 * Cols ;  /* that's the number of cols that 
					    have to be allocated per_col2 */
  per1 = tot_ones / num_rows_to_do ; 
  per2 = per_col1 + 1 ; 
  extras_r = tot_ones - per1 * num_rows_to_do ;
  /* that's the number of rows that have to be allocated per2 */

  tot_per_col = ( per2 * Rows ) / Cols ;  
  relmax = ( per1 + 1 ) * ( tot_per_col + 2 ) ;  /* was per2 * ( + 1)  */
  relatives = ivector ( 1 , relmax + 10 ) ; 
  
  printf ( "%d %d extras in cols and rows relative to baseline %d. tot_per_col = %d, per_col1 = %d, per_col2 = %d, per2 = %d\n" , extras_c , extras_r , per1 , tot_per_col , per_col1, per_col2 , per2 ) ; 

  num_c = ivector ( 1 , Cols ) ; 
  pikd_c = ivector ( 1 , Cols ) ; 
  cum = ivector ( 0 , Cols ) ; cum[0] = 0 ; 
  mylist =  ivector ( -1 , per2 ) ; 

  printf ( "Start:" ) ; fflush(stdout) ;

  /* allocate the larger numbers of 1s (per_col2) to the cols with small j */

  for ( j = 1 ; j <= Cols ; j++ )  
    num_c[j] = ( j <= extras_c ) ? per_col2 : per_col1 ; 

  if ( verbose ) {
    printf ("initial num_c:\n" ) ; 
    for ( j = 1 ; j <= Cols ; j ++ ) 
      printf ( "%2d " , num_c[j] ) ; 
    printf ("\n" ) ; 
  }

  /* Run up the rows (m) from bottom upwards, starting from Rows - c0 */

  for ( i = num_rows_to_do , m = Rows - c0 ; i >= 1 ; i-- , m-- ) {

    attempts = 0 ; 
    if ( !(i%100) ) { printf ( "%d;" , i ) ; fflush(stdout) ; }
    if ( verbose )  { printf ( " %d-" , i ) ; fflush(stdout) ; }

    weight = 0 ;
    /* compute cumulatives */
    for ( j = 1 ; j <= Cols ; j ++ ) {
      cum[j] = cum[j-1] + num_c[j] ; 
      if ( num_c[j] > 0 ) { weight ++ ; }
    }
    
    this_per = ( i > num_rows_to_do - extras_r ) ? per2 : per1 ;
    if ( weight < this_per ) {
      printf ( "failure caused by greed of earlier generations %d: %d > %d\n" , i , this_per , weight ) ; 
      status -- ; 
      break ;
    }

    /* allocate the larger numbers of 1s (per2) to the rows with large i */

    do { /* generate a row */
      if ( verbose )  { printf ( "%d," , attempts ) ; fflush(stdout) ; }

      for ( j = 1 ; j <= Cols ; j ++ ) {
	pikd_c[j] = 0 ;
      }
      for ( num_left = this_per ; 
	   num_left >= 1 ;
	   num_left -- ) {
	if ( verbose > 1 )  {
	  printf ( "status: %d,%d %d..%d %d..%d\n" , num_left , this_per ,cum[1], cum[Cols] ,num_c[1], num_c[Cols] ) ; 
	  for ( j = 1 ; j <= Cols ; j ++ ) 
	    printf ( "%2d " , j ) ; 
	  printf ( "\n" ); 
	  for ( j = 1 ; j <= Cols ; j ++ ) 
	    printf ( "%2d " , num_c[j] ) ; 
	  printf ( "\n" ); 
	  for ( j = 1 ; j <= Cols ; j ++ ) 
	    printf ( "%2d " , cum[j] ) ; 
	  printf ( "\n" ); 
	  fflush(stdout) ; 
	}

	do { /* we can get stuck here if in fact num is say (2,0,1,0,0,0,0) */
	  un = ranf() ; 
	  j = read_int_from_cum ( cum , 1 , Cols , un ) ; 
	  if ( verbose )  { printf ( "%d." , j ) ; fflush(stdout) ; }
	} while ( pikd_c[j] ) ; /* while j already picked */
	pikd_c[j] ++ ; 
	if ( verbose ) printf ( "\n" ); 
	if ( num_c[j] > 0 ) {
	  mylist[num_left] = j ; 
	} else { fprintf ( stderr , "ERROR cum failure %d %d %d\n" ,
			  j , cum[j-1] , cum[j] ) ; 	  num_left ++ ;
	       pause_for_return(); }
      }
      if ( num_left == 0 ) { /* check hd then move on */
	if ( hd > 0 ) {
	  hdfail = 0 ; 
/* 
   for each j on the list
   go and find any other rows that share this j
   check if we have already met this row today. If so, then that's 
   a failure for the traditional setting of hd
*/
	  rr = 0 ; 
	  prevrr = 0 ; 
	  for ( j = 1 ; j <= this_per && ( hdfail == 0 ) ; j ++ ) {
	    nn = mylist[j] ;
	    for ( u = num_nlist[nn] ; u >= 1 && ( hdfail == 0 ) ; u-- ) {
	      rr ++ ;
	      latest = relatives[rr] = nlist[nn][u] ; 
	      for ( ss = prevrr ; ss >= 1 ; ss -- ) {
		if ( relatives[ss] == latest ) {
		  hdfail = 1 ; 
		  break ; 
		}
	      }
	    }
	    prevrr = rr ; 
	    if ( rr > relmax ) {
	      fprintf ( stderr, "\nWARNING, relatives overran vector %d %d\n" , rr , relmax ) ;
	      status -- ; 
	    }
	  }
	}
	if ( hd > 0 && hdfail == 1 ) {
	  attempts ++ ;
	  if ( i == 1 || (attempts > maxattempts) ) { /* tough luck, have to accept the failure */
	    printf ( "H" ) ;
	    fprintf ( stderr, " warning, hamming relative %d %d out of %d\n" , rr , ss , relmax ) ;
	    hdfail = 0 ; 
	  } else {
	    printf ( "h" ) ;
	  } fflush(stdout) ;
	}
	if ( hdfail == 0 ) {
	  for ( j = 1 ; j <= this_per ; j ++ ) {
	    nn = mylist[j] ;
	    num_c[nn] -- ; 
	  }
	  /* check if we have enough weight for the next row */
	  weight = 0 ;
	  /* compute cumulatives */
	  for ( j = 1 ; j <= Cols ; j ++ ) {
	    if ( num_c[j] > 0 ) { weight ++ ; }
	  } /* check out weight of next row */
	  this_per = ( i-1 > num_rows_to_do - extras_r ) ? per2 : per1 ;
	  if ( (i>1) && ( weight < this_per ) && (!(attempts > maxattempts)) ) {
	    printf ( "s" ) ; /* reject this, do again */
	    /*  I think we can get stuck here.... */
	    for ( j = 1 ; j <= this_per ; j ++ ) {
	      nn = mylist[j] ;
	      num_c[nn] ++ ; /* put them back again! */
	    }
	    hdfail = 1 ; 
	  } else {
	    /* final acts before moving on */
	    for ( j = 1 ; j <= this_per ; j ++ ) {
	      nn = mylist[j] ;
	      add_to_alist ( a , nn , m ) ; 
	    }
	  }
	}
      } else {
	printf ( "?" ) ; /* fflush(stdout) ;  THIS should not happen any more */
      }
    } while ( num_left != 0 || hdfail != 0 ) ;
  }
  if ( verbose )  { printf ( "done uneven\n" ) ; fflush(stdout) ; }
  free_ivector ( mylist , -1 , per2 ) ;
  free_ivector ( num_c , 1 , Cols ) ;
  free_ivector ( pikd_c , 1 , Cols ) ;
  free_ivector ( cum , 0 , Cols ) ;
  free_ivector ( relatives , 1 , relmax ) ;
  printf ( "done\n" ) ; fflush(stdout) ;
  
  return status ; 
}

int 	  read_int_from_dcum ( double *dcum , int lo , int hi , double u ) {
  /* dcum[lo-1] is expected to exist and be 0
     and dcum[hi] is expected to be 1 ;
     u is expected to be in 0,1
     */
  int i , top = lo  , bot = hi ;
  int done = 0 ; 
  /* top and bot are the extreme eligible values for i */
  /* first guess */
  i = (int) ( (double) ( hi - lo + 1 ) * u ) + lo ;
  if ( i > hi ) { i = hi ; }
  if ( i < lo ) { i = lo ; }

  do {
    if ( ( top <= bot ) ||
	( ( u < dcum[i] ) && ( u > dcum[i-1] ) )
	) { /* we are done */ 
      done = 1 ; 
      break ; 
    } else if  ( u < dcum[i] ) {
      top = i-1 ;
    } else {
      bot = i+1 ;
    }
    i = ( top + bot ) / 2 ; 
  } while ( !done ) ; 


  return i ; 
} 

int 	  read_int_from_cum ( int *cum , int lo , int hi , double u ) {
  /* cum[lo-1] is expected to exist and be 0
     and cum[hi] is expected to be positive ;
     u is expected to be in 0,1
     */
  int i , top = hi  , bot = lo ;
  int done = 0 ; 
  /* top and bot are the extreme eligible values for i */
  int ONE = cum[hi] ; 
  int compare = (int) ( u * (double)(ONE) ) ; 

  /* first guess */
  i = (int) ( (double) ( hi - lo + 1 ) * u ) + lo ;
  if ( i > hi ) { i = hi ; }
  if ( i < lo ) { i = lo ; }

  do {
    if ( ( top <= bot ) ||
	( ( compare < cum[i] ) && ( compare >= cum[i-1] ) )
	) { /* we are done */ 
      done = 1 ; 
      break ; 
    } else if  ( compare < cum[i] ) {
      top = i-1 ;
    } else {
      bot = i+1 ;
    }
    i = ( top + bot ) / 2 ; 
  } while ( !done ) ; 


  return i ; 
} 

/* invert_cmatrix is in r.c */

void invert_left_hand_end_cmatrix ( unsigned char **A , unsigned char **B , int N ) {
/* takes a rectangular matrix (eg) of height N and length longer than N 
   and inverts the left hand bit, modifying the matrix to make 
   sure it is invertible. The inverse ends up in B.
   This can be used on square matrices too.
 memory should be allocated already. 
*/
  int i ;
  cm_inversion pr ;

  allocate_cm_inversion ( A , 1 , N , B , &pr ) ; 
  while ( ( i = invert_cmatrix ( &pr ) ) == 0 )  {
    printf( "m.%d." , pr.i ) ;
    i = modify_cmatrix_row ( &pr ) ;
  }
  free_cm_inversion ( &pr ) ; 
}

/* memory should be allocated already. 
   this creates an A B pair such that A is sparse */
void sparse_invertible_cmatrix (  unsigned char **A , unsigned char **B ,
				int per , int N ) {

  sparse_random_cmatrix2 ( A , per , N  ) ; 
  invert_left_hand_end_cmatrix ( A , B , N ) ;
}

int random_cvector ( unsigned char *v , double f , int lo , int hi ) {
  int j = hi , c = 0 ;

  for (  ; j >= lo ; j -- )  {
    v[j] = ( ranu() < f ) ? 1 : 0 ; 
    c += v[j] ; 
  }
  return c ; 
}

/* n is the required weight */
int fixed_wt_cvector ( unsigned char *v , int n , int lo , int hi ) {
  int N = hi - lo + 1 ; 
  int j = hi ;
  double dj = (double) N ; 
  int status = 0 ; 

  for (  ; j >= lo ; j -- , dj -- )  {
    if ( ranf() <=  (double) n / dj ) {
      n -- ;
      v[j] = 1 ;
    }
    else v[j] = 0 ; 
  }
  if ( n!= 0 ) {
    fprintf ( stderr , "warning: fixed_wt_cvector failed by %d / %d\n" ,
			n , N ) ; 
    status -- ; 
  }
  return status ; 
}