cmatrix.c

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

/* cmatrix.c
   
   My routines for handling funny representations of unsigned char 
   matrices and doing modulo 2 arithmetic. 

   Some routines are in r.c, but others are in here 

   See cmatrix.h also

*/

#include "./r.h"
#include "./rand2.h"
#include "./mynr.h"
#include "./cmatrix.h"

void alist_times_cvector_mod2 
( alist_matrix *a , unsigned char *x , unsigned char *y ) {
  int m ; 
  for ( m = 1 ; m <= a->M ; m++ ) {
    y[m] = cdotprod_mod2_index ( a->mlist[m] , x , 1 , a->num_mlist[m] ) ; 
  }
}

/* version of this for sparse x that runs through x
   only doing anything where x is 1
*/
void alist_times_cvector_sparse_mod2 
( alist_matrix *a , unsigned char *x , unsigned char *y ) {
  int n , m , i ; 
  int *nlist ; 

  for ( m = 1 ; m <= a->M ; m++ ) {
    y[m] = 0 ; 
  }
  for ( n = 1 ; n <= a->N ; n++ ) {
    if ( x[n] ) {
      nlist = a->nlist[n] ;
      for ( i = a->num_nlist[n] ; i >= 1 ; i -- ) {
	y[ nlist[i] ] ^= 1 ;
      }
    }
  }
}

/* version of this for sparse x that runs through x
   only doing anything where x is 1
*/
void alist_transpose_cvector_sparse_mod2 
( alist_matrix *a , unsigned char *x , unsigned char *y ) {
  int m , n , i ; 
  int *mlist ; 

  for ( n = 1 ; n <= a->N ; n++ ) {
    y[n] = 0 ; 
  }
  for ( m = 1 ; m <= a->M ; m++ ) {
    if ( x[m] ) {
      mlist = a->mlist[m] ;
      for ( i = a->num_mlist[m] ; i >= 1 ; i -- ) {
	y[ mlist[i] ] ^= 1 ;
      }
    }
  }
}

/* version of this for sparse x that runs through x
   only doing anything where x is 1
*/
int alist_times_ivector_sparse
( alist_matrix *a , int *x , int *y ) {
  int n , m , i , w=0 ; 
  int *nlist ; 

  for ( m = 1 ; m <= a->M ; m++ ) {
    y[m] = 0 ; 
  }
  for ( n = 1 ; n <= a->N ; n++ ) {
    if ( x[n] ) {
      nlist = a->nlist[n] ;
      for ( i = a->num_nlist[n] , w += i ; i >= 1 ; i -- ) {
	y[ nlist[i] ] ++ ;
      }
    }
  }
  return (w); /* weight of y */
}

/* version of this for sparse x that runs through x
   only doing anything where x is 1
*/
int alist_transpose_ivector_sparse
( alist_matrix *a ,  int *x , int *y ) {
  int m , n , i , w=0 ; 
  int *mlist ; 

  for ( n = 1 ; n <= a->N ; n++ ) {
    y[n] = 0 ; 
  }
  for ( m = 1 ; m <= a->M ; m++ ) {
    if ( x[m] ) {
      mlist = a->mlist[m] ;
      for ( i = a->num_mlist[m] , w+=i ; i >= 1 ; i -- ) {
	y[ mlist[i] ] ++ ;
      }
    }
  }
  return w ;
}

void alist_times_dvector
( alist_matrix *a , double *x , double *y ) {
  int n , m , i ; 
  int *nlist ; 

  for ( m = 1 ; m <= a->M ; m++ ) {
    y[m] = 0.0 ; 
  }
  for ( n = 1 ; n <= a->N ; n++ ) {
    nlist = a->nlist[n] ;
    for ( i = a->num_nlist[n] ; i >= 1 ; i -- ) {
	y[ nlist[i] ] += x[n] ;
    }
  }
}

void alist_transpose_dvector
( alist_matrix *a ,  double *x , double *y ) {
  int m , n , i ; 
  int *mlist ; 

  for ( n = 1 ; n <= a->N ; n++ ) {
    y[n] = 0 ; 
  }
  for ( m = 1 ; m <= a->M ; m++ ) {
    mlist = a->mlist[m] ;
    for ( i = a->num_mlist[m] ; i >= 1 ; i -- ) {
      y[ mlist[i] ] += x[m] ;
    }
  }
}

unsigned char cdotprod_mod2_index 
( 
 int *a , unsigned char *b , int lo , int hi ) {
  unsigned char ans = 0 ; 
  int i ; 
  for ( i = lo ; i <= hi ; i ++ ) {
    ans ^= b[a[i]] ;
  }
  return ans ; 
}


int cdotprod_index 
( 
 int *a , unsigned char *b , int lo , int hi ) {
  int ans = 0 ; 
  int i ; 
  for ( i = lo ; i <= hi ; i ++ ) {
    ans += b[a[i]] ;
  }
  return ans ;
}

void initialize_alist ( alist_matrix *a , int NN , int MM , 
			      int biggestn , int biggestm ) {
  int mm , nn ; 

  a->N = NN ; a->M = MM ; 
  a->biggest_num_m = a->biggest_num_m_alloc = biggestm ;
  a->biggest_num_n = a->biggest_num_n_alloc = biggestn ;
    
  a->num_mlist = ivector ( 1 , MM ) ; 
  a->mlist = imatrix ( 1 , MM , 1 , a->biggest_num_m ) ; 
  a->num_nlist = ivector ( 1 , NN ) ; 
  a->nlist = imatrix ( 1 , NN , 1 , a->biggest_num_n ) ; 
  a->l_up_to = ivector ( 1 , MM ) ; 
  a->u_up_to = ivector ( 1 , NN ) ; 
  
  for ( mm = 1 ; mm <= MM ; mm++ ) {
    a->num_mlist[mm] = 0 ; 
  }
  for ( nn = 1 ; nn <= NN ; nn++ ) {
    a->num_nlist[nn] = 0 ; 
  }
}

void free_alist ( alist_matrix *a ) {
  int NN = a->N , MM = a->M ; 
  free_ivector ( a->num_mlist , 1 , MM ) ; 
  free_imatrix ( a->mlist , 1 , MM , 1 , a->biggest_num_m_alloc ) ; 
  free_ivector ( a->num_nlist , 1 , NN ) ; 
  free_imatrix ( a->nlist , 1 , NN , 1 , a->biggest_num_n_alloc ) ; 
  free_ivector ( a->l_up_to , 1 , MM ) ; 
  free_ivector ( a->u_up_to , 1 , NN ) ;
}

/* should check to see if we are overflowing */
void safe_add_to_alist ( alist_matrix *a , int nn , int mm ) {
  if ( ( a->num_nlist[nn] < a->biggest_num_n ) && 
       ( a->num_mlist[mm] < a->biggest_num_m ) ) {
    a->num_nlist[nn] ++ ;
    a->num_mlist[mm] ++ ;
    a->mlist[mm][ a->num_mlist[mm] ] = nn ; 
    a->nlist[nn][ a->num_nlist[nn] ] = mm ;
  } else {
    fprintf( stderr, "%d %d overflow\n" , mm , nn ) ; 
  }
}
/* should check to see if we are overflowing */
void add_to_alist ( alist_matrix *a , int nn , int mm ) {
  a->num_nlist[nn] ++ ;
  a->num_mlist[mm] ++ ;
  a->mlist[mm][ a->num_mlist[mm] ] = nn ; 
  a->nlist[nn][ a->num_nlist[nn] ] = mm ;
}

void subtract_from_alist ( alist_matrix *a , int nn , int u ) {
  int m , l , ll ;
  /* identify m and ll */
  m = a->nlist[nn][u] ; 
  ll = 0 ; 
  for ( l = 1 ; l <= a->num_mlist[m] ; l ++ ) {
    if ( a->mlist[m][l] == nn ) {
      ll = l ;
      break ; 
    }
  }
  if ( ll == 0 ) fprintf ( stderr , "problem in alist subtract %d %d\n" , nn , u ) ; 
  /* shuffle back rest of list */
  list_slide_back ( a->mlist[m] , ll , a->num_mlist[m] ) ;
  list_slide_back ( a->nlist[nn] , u , a->num_nlist[nn] ) ;

  a->num_nlist[nn] -- ;
  a->num_mlist[m] -- ;
}

void list_slide_back ( int *list , int bot , int top ) {
  if ( bot < top ) list[bot] = list[top] ;
  list[top] = 0 ; 
}

void subtract_col_from_alist ( alist_matrix *a , int nn ) {
  int u ; 
  for ( u = a->num_nlist[nn] ; u >= 1 ; u -- ) {
    subtract_from_alist ( a , nn , u ) ;
  }
}

void kill_blank_col_from_alist ( alist_matrix *a , int nn ) {
  int u , N = a->N , m , l ; 
  int count ; 

  a->same_length = 0 ; 

  if ( a->num_nlist[nn] ) fprintf ( stderr , "warning, nonempty column %d\n" , nn ) ; 
  if ( nn == a->N ) { /* no need to do anything */
  } else {            /* move last column into this one's spot */
     u = a->num_nlist[N] ;
     a->num_nlist[nn] = u ; 
     for ( ; u >= 1 ; u -- ) { 
       m = a->nlist[N][u] ;
       a->nlist[nn][u] = m ; 
       /* get item m, and change his entry 
	  mlist[m][l] from N to nn */
       count = 0 ; 
       for ( l = a->num_mlist[m] ; l >= 1 ; l -- ) {
	 if ( a->mlist[m][l] == N ) { a->mlist[m][l] = nn ; count++ ; }
       }
       if ( count != 1 ) fprintf ( stderr , "warning, expected 1 count, but at %d it's %d\n" , nn , count ) ; 
     }
  }
  a->N -- ; 
}

void write_alist ( FILE *fp , alist_matrix *a ) {
  /* this assumes that mlist and nlist have the form of a rectangular
     matrix in the file; if lists have unequal lengths, then the 
     entries should be present but are ignored
     */
  int N = a->N , M = a->M ; 

  fprintf ( fp , "%d %d\n" , N , M ) ; 
  fprintf ( fp , "%d %d\n" , a->biggest_num_n , a->biggest_num_m ) ; 
  write_ivector ( fp , a->num_nlist , 1 , N ) ; 
  write_ivector ( fp , a->num_mlist , 1 , M ) ; 
  write_imatrix ( fp , a->nlist , 1 , N , 1 , a->biggest_num_n ) ; 
  write_imatrix ( fp , a->mlist , 1 , M , 1 , a->biggest_num_m ) ; 
}  

void write_alist_transpose ( FILE *fp , alist_matrix *a ) {
  int N = a->N , M = a->M ; 

  fprintf ( fp , "%d %d\n" , M , N ) ; 
  fprintf ( fp , "%d %d\n" , a->biggest_num_m , a->biggest_num_n ) ; 
  write_ivector ( fp , a->num_mlist , 1 , M ) ; 
  write_ivector ( fp , a->num_nlist , 1 , N ) ; 
  write_imatrix ( fp , a->mlist , 1 , M , 1 , a->biggest_num_m ) ; 
  write_imatrix ( fp , a->nlist , 1 , N , 1 , a->biggest_num_n ) ; 
}  

int read_allocate_alist ( alist_matrix *a , char *file ) { 
  /* this assumes that mlist and nlist have the form of a rectangular
     matrix in the file; if lists have unequal lengths, then the 
     entries should be present but are ignored
     */
  int status = 0 ; 
  int N , M ; 
  int biggestn , biggestm ;
  FILE    *fp;
  
  fp = fopen( file, "r" );
  if( !fp )   fprintf( stderr, "No such file: %s\n", file ), exit(0);

  do {
    if ( fscanf(fp,"%d %d " , &N , &M ) == EOF ) {
      status = -1 ;       break ;
    }
    fprintf(stderr,"-");  fflush(stderr);
    if ( fscanf(fp,"%d %d " , &biggestn , &biggestm ) == EOF ) {
      status = -1 ;       break ;
    }

    initialize_alist ( a , N , M ,  biggestn ,  biggestm ) ;

    status += fread_ivector ( a->num_nlist , 1 , N , fp ) ; 
    status += fread_ivector ( a->num_mlist , 1 , M , fp ) ; 
    status += fread_imatrix ( a->nlist , 1 , N , 1 , biggestn , fp ) ;
    status += fread_imatrix ( a->mlist , 1 , M , 1 , biggestm , fp ) ;
    /*    fprintf(stderr,";");  fflush(stderr); */
  } while ( 0 ) ; 

  if ( status < 0 ) { fprintf(stderr,"\nread_allocate_alist status %d\n",status); }
  fprintf(stderr,"Z");  fflush(stderr);
  return status ; 
}

int code0_matrix_alist ( alist_matrix *a , int N , int M , int tr ) {
  int m , n , num_left ; 
  int status = 0 ; 
  
  if ( M < N ) {
    fprintf ( stderr , "code0 matrix problem with N and M\n" ) ; 
    pause_for_return () ;
    exit (0) ;
  }
  for ( n = 1 ; n <= N ; n++ ) {
    m = n ; 
    add_to_alist ( a , n , m ) ;
  }
  for ( m = 1 ; m <= M ; m++ ) { /* exact number per row */
    n = N ;
    num_left = tr ;
    for ( ; n >= 1 ; n-- ) {
      if ( num_left > 0 && ( ranf() <=  (double) num_left / (double) n ) ) {
	add_to_alist ( a , n , m ) ;
	num_left -- ;
      } 
    }
  }
  status += finish_off_alist ( a ) ;
  return status ; 
}

/*
   square matrices only 
*/
int cmatrix_2_alist ( unsigned char **A , alist_matrix *a , int NN ) {
  int mm , nn , uu ; 
  int status = 0 ; 
  
  for ( nn = 1 ; nn <= NN ; nn++ ) {
    for ( uu = 1 , mm = 1 ; mm <= NN ; mm++ ) {
      if ( A[mm][nn] ) {
	a->nlist[nn][uu] = mm ;
	uu ++ ; 
	a->num_mlist[mm] ++ ;
	a->mlist[mm][ a->num_mlist[mm] ] = nn ; 
      }
    }
    uu -- ;
    a->num_nlist[nn] = uu ;
  }
  status += finish_off_alist ( a ) ;
  return status ; 
}

int gen_cmatrix_2_alist ( unsigned char **A , alist_matrix *a , int M , int N , int transpose ) {
  int mm , nn , uu ; 
  int status = 0 ; 
  
  for ( nn = 1 ; nn <= N ; nn++ ) {
    for ( uu = 1 , mm = 1 ; mm <= M ; mm++ ) {
      if ( ( transpose && A[nn][mm] ) || ( !transpose && A[mm][nn] ) ) {
	a->nlist[nn][uu] = mm ;
	uu ++ ; 
	a->num_mlist[mm] ++ ;
	a->mlist[mm][ a->num_mlist[mm] ] = nn ; 
      }
    }
    uu -- ;
    a->num_nlist[nn] = uu ;
  }
  status += finish_off_alist ( a ) ;
  return status ; 
}

/* 
   puts them side by side if vertical = 0, 
   or on top of each other if 1.
                                                */
int cmatrices_2_alist ( unsigned char **A ,  unsigned char **B ,
			       alist_matrix *a , int NN , int vertical ) {
  int mm , nn , mat , tn , tm , moffset ; 
  unsigned char **C ;
  int status = 0 ; 

  tn = 1 ;  moffset = 0 ; 
  for ( C = A , mat = 1 ; 
       mat <= 2 ; 
       ) {
    for ( nn = 1 ; nn <= NN ; nn++ , tn++ ) {
      for ( mm = 1 ; mm <= NN ; mm++ ) {
	if ( C[mm][nn] ) {
	  tm = mm + moffset ;
	  add_to_alist ( a , tn , tm ) ; 
	}
      }
    }
    /* end of loop */
    mat++ ; C = B ; 
    if ( vertical ) {
      tn = 1 ;
      moffset = NN ;
    }
  }
  status += finish_off_alist ( a ) ;
  return status ; 
}

/*
   puts identity matrix alongside the provided matrix 
   [ C I ]
*/
int cmatrix_andI_2_alist ( unsigned char **A ,  
			       alist_matrix *a , int NN , int vertical ) {
  int mm , nn ; 
  int status = 0 ; 

  for ( nn = 1 ; nn <= NN ; nn++ ) {
    for ( mm = 1 ; mm <= NN ; mm++ ) {
      if ( A[mm][nn] ) {
	add_to_alist ( a , nn , mm ) ; 
      }
    }
  }
  if ( !vertical ) {
    for ( mm = 1, nn=NN+1 ; mm <= NN ; mm++ , nn++ ) {
      add_to_alist ( a , nn , mm ) ;
    }
  } else {
    for ( mm = NN+1, nn=1 ; nn <= NN ; mm++ , nn++ ) {
      add_to_alist ( a , nn , mm ) ;
    }
  }

  status += finish_off_alist ( a ) ;
  return status ; 
}

int finish_off_alist ( alist_matrix *a ) {
  int nn , mm , uu ; 
  int tbm = 0  , tbn = 0  ; /* keeps track of the true biggest_m , n */
  int status = 0 ; 

  for ( mm  = 1 ; mm <= a->M ; mm++ ) {