//
// main.c
// Test2
//
// Created by lianhua on 2014/9/20.
// Copyright (c) 2014年 lianhua. All rights reserved.
//
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
void binary(int N, const int n, int *b_vec)
{
int *ptr ;
int i ;
ptr
= (int*)malloc(sizeof(int)*n
); for(i = 1; i <= n ; i++){
ptr[n-i] = N%2 ;
N = N/2;
}
for(i = 0; i < n; ++i){
b_vec[i]=ptr[i];
}
}
void vectorM(int vec[], const int n, int m, int DM[n][m]){
int i,j;
int Coding_array[n] ;
for(j = 0; j < m; j++){
binary(vec[j],n,Coding_array);
for(i = 0; i<n; i++){
int *ptr = Coding_array;
DM[i][j] = 2*ptr[i]-1;
}
}
}
void Ltri1(const int m, int Matrix[m][m]){
int i, j ;
for(i = 0; i < m; i++){
for(j = 0; j < m; j++){
if(i <= j){
Matrix[i][j] = 0;
}else{
Matrix[i][j] = 1;
}
}
}
}
void TransM(const int m, const int n, int Mi[][m], int Mo[][m]){
int h, i, j ;
for(i = 0; i<m;i++){
for(j = 0; j<m; j++){
Mo[i][j] = 0 ;
for(h = 0; h<n; h++){
Mo[i][j] += Mi[h][i]*Mi[h][j];
}
}
}
}
static int n_step(int n)
{
if (n <= 1) return 1;
return n * n_step(n - 1);
}
double Es2(const int n, const int m, int vec[]){
int i, j;
int k = 0;
float a,b ;
int DM[n][m] ;
int TM[m][m] ;
int SM[m][m] ;
int VL[m][m] ;
vectorM(vec,n,m,DM);
Ltri1(m,TM);
TransM(m,n,DM,SM) ;
for(i = 0; i<m;i++){
for(j = 0; j<m; j++){
VL[i][j] = 0 ;
VL[i][j] = SM[i][j]*TM[i][j];
}
}
for(i = 0; i<m;i++){
for(j = 0; j<m; j++){
k = k+(VL[i][j]*VL[i][j]);
}
}
b = (n_step(m)/(n_step(2)*n_step(m-2))) ;
a = k/b ;
return a;
}
int find_index(int vec[], int a){
int j = 0;
while(vec[j]!=a){
j++;
}
return j;
}
int find_indexf(float vec[],float a){
int j = 0;
while(vec[j]!=a){
j++;
}
return j;
}
void array_delete(int vec[], int m, int b){
int j;
for(j = b; j< m-1;j++){
vec[j] = vec[j+1];
}
}
void array_insert(int vec[], int m, int p, int value){
int j ;
for(j = m-1; j>p; j--){
vec[j+1] = vec[j];
}
vec[p] = value;
}
void Delcol(const int n, const int m ,int vec[]){
int i, j;
int DM[n][m] ;
int TM[m][m] ;
int SM[m][m] ;
int delta_s2[m] ;
vectorM(vec,n,m,DM);
Ltri1(m,TM);
TransM(m,n,DM,SM) ;
for(i = 0; i<m;i++){
for(j = 0; j<m; j++){
if(SM[i][j]==n) SM[i][j] = 0 ;
else SM[i][j] = SM[i][j];
}
}
for(i = 0; i < m;i++){
int sum = 0;
for(j = 0; j<m; j++){
sum += SM[i][j]*SM[i][j] ;
}
delta_s2[i] = sum ;
}
int max = delta_s2[0] ;
for(i = 0; i < m;i++){
if(delta_s2[i] > max) max = delta_s2[i];
}
int idx = find_index(delta_s2,max);
array_delete(vec,m,idx);
}
void addcol(int vec_old[], int vec_ch[], int m, const int n, int rLB, int k){
int h, i, j ;
m = m-rLB+k;
int m_star = rLB-k ;
for(i = 0; i<m;i++){
int idx ;
idx = find_index(vec_ch,vec_old[i]);
array_delete(vec_ch,m,idx);
}
int DM1[n][m] ;
vectorM(vec_old,n,m,DM1);
int DM2[n][m_star] ;
vectorM(vec_ch,n,m_star,DM2);
int Mo[m][m_star] ;
for(i = 0; i<m;i++){
for(j = 0; j<m_star; j++){
Mo[i][j] = 0 ;
for(h = 0; h<n; h++){
Mo[i][j] += DM1[h][i]*DM2[h][j];
}
}
}
int delta_s2[m_star];
for(i = 0; i < m_star;i++){
int sum = 0;
for(j = 0; j<m; j++){
sum += Mo[j][i]*Mo[j][i] ;
}
delta_s2[i] = sum ;
}
int min = delta_s2[0] ;
for(i = 0; i < m_star;i++){
if(delta_s2[i] < min) min = delta_s2[i];
}
int idx = find_index(delta_s2,min);
array_insert(vec_old,m,m,vec_ch[idx]);
}
int main(int argc, char *argv[])
{
int i;
int m = 6 ;
int n = 4 ;
int rLB = 3;
int vec_old[3] = {6,12,5};
int vec_ch[6] = {10,9,3,6,12,5};
for(i = 0; i<rLB; i++) addcol(vec_old, vec_ch, m, n, rLB, i);
}
Ly8KLy8gIG1haW4uYwovLyAgVGVzdDIKLy8KLy8gIENyZWF0ZWQgYnkgbGlhbmh1YSBvbiAyMDE0LzkvMjAuCi8vICBDb3B5cmlnaHQgKGMpIDIwMTTlubQgbGlhbmh1YS4gQWxsIHJpZ2h0cyByZXNlcnZlZC4KLy8KCiNpbmNsdWRlIDxzdGRpby5oPgojaW5jbHVkZSA8c3RkbGliLmg+CiNpbmNsdWRlIDx0aW1lLmg+CiNpbmNsdWRlIDxtYXRoLmg+Cgp2b2lkIGJpbmFyeShpbnQgTiwgY29uc3QgaW50IG4sIGludCAqYl92ZWMpCnsKICAgIGludCAqcHRyIDsKICAgIGludCBpIDsKICAgIHB0ciA9IChpbnQqKW1hbGxvYyhzaXplb2YoaW50KSpuKTsKICAgIGZvcihpID0gMTsgaSA8PSBuIDsgaSsrKXsKICAgICAgICBwdHJbbi1pXSA9IE4lMiA7CiAgICAgICAgTiA9IE4vMjsKICAgIH0KICAgIGZvcihpID0gMDsgaSA8IG47ICsraSl7CiAgICAgICAgYl92ZWNbaV09cHRyW2ldOwogICAgfQogICAgZnJlZShwdHIpIDsKfQoKdm9pZCB2ZWN0b3JNKGludCB2ZWNbXSwgY29uc3QgaW50IG4sIGludCBtLCBpbnQgRE1bbl1bbV0pewogICAgaW50IGksajsKICAgIGludCBDb2RpbmdfYXJyYXlbbl0gOwogICAgZm9yKGogPSAwOyBqIDwgbTsgaisrKXsKICAgICAgICBiaW5hcnkodmVjW2pdLG4sQ29kaW5nX2FycmF5KTsKICAgICAgICBmb3IoaSA9IDA7IGk8bjsgaSsrKXsKICAgICAgICAgICAgaW50ICpwdHIgPSBDb2RpbmdfYXJyYXk7CiAgICAgICAgICAgIERNW2ldW2pdID0gMipwdHJbaV0tMTsKICAgICAgICB9CiAgICB9Cn0KCnZvaWQgTHRyaTEoY29uc3QgaW50IG0sIGludCBNYXRyaXhbbV1bbV0pewogICAgaW50IGksIGogOwogICAgZm9yKGkgPSAwOyBpIDwgbTsgaSsrKXsKICAgICAgICBmb3IoaiA9IDA7IGogPCBtOyBqKyspewogICAgICAgICAgICBpZihpIDw9IGopewogICAgICAgICAgICAgICAgTWF0cml4W2ldW2pdID0gMDsKICAgICAgICAgICAgfWVsc2V7CiAgICAgICAgICAgICAgICBNYXRyaXhbaV1bal0gPSAxOwogICAgICAgICAgICB9CiAgICAgICAgfQogICAgfQp9Cgp2b2lkIFRyYW5zTShjb25zdCBpbnQgbSwgY29uc3QgaW50IG4sIGludCBNaVtdW21dLCBpbnQgTW9bXVttXSl7CiAgICBpbnQgaCwgaSwgaiA7CiAgICBmb3IoaSA9IDA7IGk8bTtpKyspewogICAgICAgIGZvcihqID0gMDsgajxtOyBqKyspewogICAgICAgICAgICBNb1tpXVtqXSA9IDAgOwogICAgICAgICAgICBmb3IoaCA9IDA7IGg8bjsgaCsrKXsKICAgICAgICAgICAgICAgIE1vW2ldW2pdICs9IE1pW2hdW2ldKk1pW2hdW2pdOwogICAgICAgICAgICB9CiAgICAgICAgfQogICAgfQp9CgoKc3RhdGljIGludCBuX3N0ZXAoaW50IG4pCnsKICAgIGlmIChuIDw9IDEpIHJldHVybiAxOwogICAgcmV0dXJuIG4gKiBuX3N0ZXAobiAtIDEpOwp9Cgpkb3VibGUgRXMyKGNvbnN0IGludCBuLCBjb25zdCBpbnQgbSwgaW50IHZlY1tdKXsKICAgIGludCBpLCBqOwogICAgaW50IGsgPSAwOwogICAgZmxvYXQgYSxiIDsKICAgIGludCBETVtuXVttXSA7CiAgICBpbnQgVE1bbV1bbV0gOwogICAgaW50IFNNW21dW21dIDsKICAgIGludCBWTFttXVttXSA7CiAgICB2ZWN0b3JNKHZlYyxuLG0sRE0pOwogICAgTHRyaTEobSxUTSk7CiAgICBUcmFuc00obSxuLERNLFNNKSA7CiAgICBmb3IoaSA9IDA7IGk8bTtpKyspewogICAgICAgIGZvcihqID0gMDsgajxtOyBqKyspewogICAgICAgICAgICBWTFtpXVtqXSA9IDAgOwogICAgICAgICAgICBWTFtpXVtqXSA9IFNNW2ldW2pdKlRNW2ldW2pdOwogICAgICAgIH0KICAgIH0KICAgIGZvcihpID0gMDsgaTxtO2krKyl7CiAgICAgICAgZm9yKGogPSAwOyBqPG07IGorKyl7CiAgICAgICAgICAgIGsgPSBrKyhWTFtpXVtqXSpWTFtpXVtqXSk7CiAgICAgICAgfQogICAgfQogICAgYiA9IChuX3N0ZXAobSkvKG5fc3RlcCgyKSpuX3N0ZXAobS0yKSkpIDsKICAgIGEgPSBrL2IgOwogICAgcmV0dXJuIGE7Cn0KCmludCBmaW5kX2luZGV4KGludCB2ZWNbXSwgaW50IGEpewogICAgaW50IGogPSAwOwogICAgd2hpbGUodmVjW2pdIT1hKXsKICAgICAgICBqKys7CiAgICB9CiAgICByZXR1cm4gajsKfQoKaW50IGZpbmRfaW5kZXhmKGZsb2F0IHZlY1tdLGZsb2F0IGEpewogICAgaW50IGogPSAwOwogICAgd2hpbGUodmVjW2pdIT1hKXsKICAgICAgICBqKys7CiAgICB9CiAgICByZXR1cm4gajsKfQoKCnZvaWQgYXJyYXlfZGVsZXRlKGludCB2ZWNbXSwgaW50IG0sIGludCBiKXsKICAgIGludCBqOwogICAgZm9yKGogPSBiOyBqPCBtLTE7aisrKXsKICAgICAgICB2ZWNbal0gPSB2ZWNbaisxXTsKICAgIH0KfQoKdm9pZCBhcnJheV9pbnNlcnQoaW50IHZlY1tdLCBpbnQgbSwgaW50IHAsIGludCB2YWx1ZSl7CiAgICBpbnQgaiA7CiAgICBmb3IoaiA9IG0tMTsgaj5wOyBqLS0pewogICAgICAgIHZlY1tqKzFdID0gdmVjW2pdOwogICAgfQogICAgdmVjW3BdID0gdmFsdWU7Cn0KCnZvaWQgRGVsY29sKGNvbnN0IGludCBuLCBjb25zdCBpbnQgbSAsaW50IHZlY1tdKXsKICAgIGludCBpLCBqOwogICAgaW50IERNW25dW21dIDsKICAgIGludCBUTVttXVttXSA7CiAgICBpbnQgU01bbV1bbV0gOwogICAgaW50IGRlbHRhX3MyW21dIDsKICAgIHZlY3Rvck0odmVjLG4sbSxETSk7CiAgICBMdHJpMShtLFRNKTsKICAgIFRyYW5zTShtLG4sRE0sU00pIDsKICAgIGZvcihpID0gMDsgaTxtO2krKyl7CiAgICAgICAgZm9yKGogPSAwOyBqPG07IGorKyl7CiAgICAgICAgICAgIGlmKFNNW2ldW2pdPT1uKSBTTVtpXVtqXSA9IDAgOwogICAgICAgICAgICBlbHNlIFNNW2ldW2pdID0gU01baV1bal07CiAgICAgICAgfQogICAgfQogICAgZm9yKGkgPSAwOyBpIDwgbTtpKyspewogICAgICAgIGludCBzdW0gPSAwOwogICAgICAgIGZvcihqID0gMDsgajxtOyBqKyspewogICAgICAgICAgICBzdW0gKz0gU01baV1bal0qU01baV1bal0gOwogICAgICAgIH0KICAgICAgICBkZWx0YV9zMltpXSA9IHN1bSA7CiAgICB9CiAgICBpbnQgbWF4ID0gZGVsdGFfczJbMF0gOwogICAgZm9yKGkgPSAwOyBpIDwgbTtpKyspewogICAgICAgIGlmKGRlbHRhX3MyW2ldID4gbWF4KSBtYXggPSBkZWx0YV9zMltpXTsKICAgIH0KICAgIGludCBpZHggPSBmaW5kX2luZGV4KGRlbHRhX3MyLG1heCk7CiAgICBhcnJheV9kZWxldGUodmVjLG0saWR4KTsKfQoKdm9pZCBhZGRjb2woaW50IHZlY19vbGRbXSwgaW50IHZlY19jaFtdLCBpbnQgbSwgY29uc3QgaW50IG4sIGludCByTEIsIGludCBrKXsKICAgIGludCBoLCBpLCBqIDsKICAgIG0gPSBtLXJMQitrOwogICAgaW50IG1fc3RhciA9IHJMQi1rIDsKICAgIGZvcihpID0gMDsgaTxtO2krKyl7CiAgICAgICAgaW50IGlkeCA7CiAgICAgICAgaWR4ID0gZmluZF9pbmRleCh2ZWNfY2gsdmVjX29sZFtpXSk7CiAgICAgICAgYXJyYXlfZGVsZXRlKHZlY19jaCxtLGlkeCk7CiAgICB9CiAgICBpbnQgRE0xW25dW21dIDsKICAgIHZlY3Rvck0odmVjX29sZCxuLG0sRE0xKTsKICAgIGludCBETTJbbl1bbV9zdGFyXSA7CiAgICB2ZWN0b3JNKHZlY19jaCxuLG1fc3RhcixETTIpOwogICAgaW50IE1vW21dW21fc3Rhcl0gOwogICAgZm9yKGkgPSAwOyBpPG07aSsrKXsKICAgICAgICBmb3IoaiA9IDA7IGo8bV9zdGFyOyBqKyspewogICAgICAgICAgICBNb1tpXVtqXSA9IDAgOwogICAgICAgICAgICBmb3IoaCA9IDA7IGg8bjsgaCsrKXsKICAgICAgICAgICAgICAgIE1vW2ldW2pdICs9IERNMVtoXVtpXSpETTJbaF1bal07CiAgICAgICAgICAgIH0KICAgICAgICB9CiAgICB9CiAgICBpbnQgZGVsdGFfczJbbV9zdGFyXTsKICAgIGZvcihpID0gMDsgaSA8IG1fc3RhcjtpKyspewogICAgICAgIGludCBzdW0gPSAwOwogICAgICAgIGZvcihqID0gMDsgajxtOyBqKyspewogICAgICAgICAgICBzdW0gKz0gTW9bal1baV0qTW9bal1baV0gOwogICAgICAgIH0KICAgICAgICBkZWx0YV9zMltpXSA9IHN1bSA7CiAgICB9CiAgICBpbnQgbWluID0gZGVsdGFfczJbMF0gOwogICAgZm9yKGkgPSAwOyBpIDwgbV9zdGFyO2krKyl7CiAgICAgICAgaWYoZGVsdGFfczJbaV0gPCBtaW4pIG1pbiA9IGRlbHRhX3MyW2ldOwogICAgfQogICAgaW50IGlkeCA9IGZpbmRfaW5kZXgoZGVsdGFfczIsbWluKTsKICAgIGFycmF5X2luc2VydCh2ZWNfb2xkLG0sbSx2ZWNfY2hbaWR4XSk7Cn0KCgppbnQgbWFpbihpbnQgYXJnYywgY2hhciAqYXJndltdKQp7CiAgICAKICAgIGludCBpOwogICAgaW50IG0gPSA2IDsKICAgIGludCBuID0gNCA7CiAgICBpbnQgckxCID0gMzsKICAgIGludCB2ZWNfb2xkWzNdID0gezYsMTIsNX07CiAgICBpbnQgdmVjX2NoWzZdID0gezEwLDksMyw2LDEyLDV9OwogICAgZm9yKGkgPSAwOyBpPHJMQjsgaSsrKSBhZGRjb2wodmVjX29sZCwgdmVjX2NoLCBtLCBuLCByTEIsIGkpOwogICAgcHJpbnRmKCIlM2QiLHZlY19vbGQpOwogICAgcHJpbnRmKCJcbiIpOwoKICAgIAp9Cg==