shrink.c

   1 #include <stdio.h>
   2 #include <stdlib.h>
   3 #include <string.h>
   4 #include <math.h>
   5 #include <complex.h>
   6
   7 #include "matrix.h"
   8 #include "shrink.h"
   9
  10 // linked list
  11 struct Node {
  12   int data;
  13   struct Node* next;
  14 };
  15
  16 void printLinkedList(struct Node *node);
  17 void append(struct Node** head_ref, int new_data);
  18 void freeList(struct Node** head_ref);
  19
  20 /****************************************************************************
  21  * Attempts to shrinks the dimension of the affine space by one dimension.
  22  * The affine space is specified by
  23  * the first 'k' columns of the 'n'x'n' matrix 'G', translated by 'h',
  24  * that have an inner product with 'xi' of 'alpha'. 'GBar' is 'G^(-1)^T',
  25  * and ('Q', 'D', 'J') specify the quadratic form on this affine space.
  26  * ---
  27  * If result is an EMPTY affine space then return value is 'E'.
  28  * If result is the SAME affine space then return value is 'A'.
  29  * If result is a SUCCESSfully shrunk affine space then return value is 'U'.
  30  * ***
  31  * Note that D, and J are assumed to be nx1 and nxn even though they are
  32  * actually kx1 and kxk (extra elements are ignored).
  33  ****************************************************************************/
  34
  35 char shrink(int n, int *k, int *h, int **G, int **GBar, int *Q, int **D, int ***J, int *xi, int alpha) {
  36   // Note that D and J must be dynamically allocated arrays as they can change size in this function. This is also the reason that we need their base pointer and are greater in pointer depth by one.
  37
  38   int a, b, c, d;
  39   int beta;
  40
  41   struct Node* setS = NULL;
  42   struct Node* setWalker = NULL;
  43
  44   for(a=0; a<*k; a++) {
  45     if(dotProductMod(xi, G[a], n, 2) == 1)
  46       append(&setS, a);
  47   }
  48
  49   beta = (alpha + dotProductMod(xi, h, n, 2))%2;
  50
  51   if(setS == NULL) {
  52     if(beta == 1) {
  53       freeList(&setS);
  54       return 'E'; // EMPTY
  55     }
  56     else {// beta = 0
  57       freeList(&setS);
  58       return 'A'; // SAME
  59     }
  60   }
  61
  62   // we define tmpVec, tmpMatrix, R and RT now for later use since we are sure that *k!=0 (since otherwise setS==NULL) and so they will be finite-sized
  63   int* tmpVec; // temporary vector used in swapping G and GBar rows and in Eq. 49 for updating D
  64   int** tmpMatrix; // temporary matrix used to store intermediate value of J = R times J times R^T
  65   tmpMatrix = calloc(*k, sizeof(int*));
  66   for(a=0; a<*k; a++)
  67     tmpMatrix[a] = calloc(*k, sizeof(int));
  68
  69   int** R;  // basis change matrix R and R^T
  70   int** RT;
  71   R = calloc(*k, sizeof(int*));
  72   RT = calloc(*k, sizeof(int*));
  73   // initially set it to the k-dimensional identity matrix
  74   for(a=0; a<*k; a++) {
  75     R[a] = calloc(*k, sizeof(int));
  76     R[a][a] = 1;
  77     RT[a] = calloc(*k, sizeof(int));
  78     RT[a][a] = 1;
  79   }
  80
  81   int i = setS->data; // pick first element in setS to be added to the new space's shift vector
  82   // remove from setS
  83   if(setS->next != NULL) {
  84     setWalker = setS->next;
  85     setS->data = setS->next->data;
  86     setS->next = setS->next->next;
  87     free(setWalker);
  88   } else
  89     setS = NULL;
  90
  91   int tmp;
  92   setWalker = setS;
  93   while(setWalker != NULL) {
  94     // update G rows, g
  95     for(a=0; a<n; a++)
  96       G[setWalker->data][a] = (G[setWalker->data][a] + G[i][a])%2;
  97
  98     // update D and J
  99     R[setWalker->data][i] = 1;
 100     RT[i][setWalker->data] = 1;
 101
 102     tmpVec = calloc(*k, sizeof(int));
 103     for(c=0;c<*k;c++) {
 104       tmp = 0;
 105     /*   tmpVec[c] = 0; */
 106     /*   for(a=0;a<*k;a++) */
 107     /*  tmpVec[c] = tmpVec[c] + R[c][a]*((*D)[a]); */
 108     /*   for(a=0;a<*k;a++) */
 109     /*  for(b=a+1;b<*k;b++) */
 110     /*    tmpVec[c] = tmpVec[c] + (*J)[a][b]*R[c][a]*R[c][b]; */
 111     /*   tmpVec[c] = tmpVec[c]%8; */
 112       if(c!=i)
 113         tmp += R[c][c]*((*D)[c]);
 114       tmp += R[c][i]*((*D)[i]);
 115       if(i>c)
 116         tmp += ((*J)[c][i])*R[c][c]*R[c][i];
 117       else if(i<c)
 118         tmp += ((*J)[i][c])*R[c][i]*R[c][c];
 119       tmpVec[c] = tmp%8;
 120     }
 121     memcpy(*D, tmpVec, sizeof(int)*(*k));
 122     free(tmpVec);
 123
 124     // Eq. 50 update of J
 125     for(c=0;c<*k;c++) {
 126       for(d=0;d<*k; d++) {
 127         tmp = 0;
 128         if(c!=i) {
 129           if(d!=i) {
 130             tmp += R[c][c]*((*J)[c][d])*R[d][d];
 131             tmp += R[c][c]*((*J)[c][i])*R[d][i];
 132             tmp += R[c][i]*((*J)[i][d])*R[d][d];
 133           }else {
 134             tmp += R[c][c]*((*J)[c][i])*R[d][i];
 135           }
 136         } else {
 137           if(d!=i) {
 138             tmp += R[c][i]*((*J)[i][d])*R[d][d];
 139           }
 140         }
 141         tmp += R[c][i]*((*J)[i][i])*R[d][i];
 142         tmp %= 8;
 143         tmpMatrix[c][d] = tmp;
 144       }
 145     }
 146     for(c=0; c<*k; c++)
 147       memcpy(((*J)[c]), tmpMatrix[c], sizeof(int)*(*k));
 148     /* multMatrixMod(*J,RT,tmpMatrix,*k,*k,*k,*k,8); */
 149     /* multMatrixMod(R,tmpMatrix,*J,*k,*k,*k,*k,8); */
 150
 151     if(setWalker->data != i) {
 152       R[setWalker->data][i] = 0; // reset R
 153       RT[i][setWalker->data] = 0;
 154     }
 155
 156     // update GBar rows, gbar
 157     for(a=0; a<n; a++)
 158       GBar[i][a] = (GBar[i][a] + GBar[setWalker->data][a])%2;
 159
 160     setWalker = setWalker->next;
 161   }
 162
 163   // Swap 'i' and 'k'
 164   R[*k-1][*k-1] = 0; R[i][i] = 0;
 165   R[i][*k-1] = 1; R[*k-1][i] = 1;
 166   RT[*k-1][*k-1] = 0; RT[i][i] = 0;
 167   RT[i][*k-1] = 1; RT[*k-1][i] = 1;
 168
 169   tmpVec = G[i];
 170   G[i] = G[*k-1];
 171   G[*k-1] = tmpVec;
 172   tmpVec = GBar[i];
 173   GBar[i] = GBar[*k-1];
 174   GBar[*k-1] = tmpVec;
 175
 176   // update D and J
 177   //(Note: can't quite directly use previous faster algorithm since R[c][c]!=1 forall c. However, this can still be brought down by a factor of N by the same algorithm appropriately tweaked.)
 178   if(i!=(*k-1)) {
 179     tmpVec = calloc(*k, sizeof(int));
 180     for(c=0;c<*k;c++) {
 181       tmp = 0;
 182     /* tmpVec[c] = 0; */
 183     /* for(a=0;a<*k;a++) */
 184     /*   tmpVec[c] = tmpVec[c] + R[c][a]*((*D)[a]); */
 185     /* for(a=0;a<*k;a++) { */
 186     /*   for(b=a+1;b<*k;b++) */
 187     /*  tmpVec[c] = tmpVec[c] + (*J)[a][b]*R[c][a]*R[c][b]; */
 188     /* } */
 189     /* tmpVec[c] = tmpVec[c]%8; */
 190       tmp += R[c][c]*((*D)[c]); // iff c!=i and c!=(*k-1) this can be non-zero
 191       tmp += R[c][*k-1]*((*D)[*k-1]); // iff c=i this can be non-zero
 192       tmp += R[c][i]*((*D)[i]); // iff c=*k-1 this can be non-zero
 193
 194     /* if(i>c) */
 195     /*   tmp += ((*J)[c][i])*R[c][c]*R[c][i]; */
 196     /* else if(i<c) */
 197     /*   tmp += ((*J)[i][c])*R[c][i]*R[c][c]; */
 198     /* if((*k-1)>c) */
 199     /*   tmp += ((*J)[c][*k-1])*R[c][c]*R[c][*k-1]; */
 200     /* else if((*k-1)<c) // pretty sure this can't happen */
 201     /*   tmp += ((*J)[*k-1][c])*R[c][*k-1]*R[c][c]; */
 202     /* if((*k-1)>i) */
 203     /*   tmp += ((*J)[i][*k-1])*R[c][i]*R[c][*k-1]; */
 204     /* else if((*k-1)<i) */
 205     /*   tmp += ((*J)[*k-1][i])*R[c][*k-1]*R[c][i]; */
 206       tmpVec[c] = tmp%8;
 207     }
 208     memcpy(*D, tmpVec, sizeof(int)*(*k));
 209     free(tmpVec);
 210
 211     // Eq. 50 update of J
 212     for(c=0;c<*k;c++) {
 213       for(d=0;d<*k; d++) {
 214         tmp = 0;
 215         if(c!=(*k-1) && c!=i) {
 216           if(d!=(*k-1) && d!=i) {
 217             tmp += R[c][c]*((*J)[c][d])*R[d][d];
 218           }else {
 219             tmp += R[c][c]*((*J)[c][i])*R[d][i];
 220             tmp += R[c][c]*((*J)[c][*k-1])*R[d][*k-1];
 221           }
 222         } else {
 223           if(d!=(*k-1) && d !=i) {
 224             tmp += R[c][i]*((*J)[i][d])*R[d][d];
 225             tmp += R[c][*k-1]*((*J)[*k-1][d])*R[d][d];
 226           } else {
 227             tmp += R[c][i]*((*J)[i][*k-1])*R[d][*k-1];
 228             tmp += R[c][*k-1]*((*J)[*k-1][i])*R[d][i];
 229             tmp += R[c][i]*((*J)[i][i])*R[d][i];
 230             tmp += R[c][*k-1]*((*J)[*k-1][*k-1])*R[d][*k-1];
 231           }
 232         }
 233         tmpMatrix[c][d] = tmp%8;
 234       }
 235     }
 236     for(c=0; c<*k; c++)
 237       memcpy(((*J)[c]), tmpMatrix[c], sizeof(int)*(*k));
 238   /* multMatrixMod(*J,RT,tmpMatrix,*k,*k,*k,*k,8); */
 239   /* multMatrixMod(R,tmpMatrix,*J,*k,*k,*k,*k,8); */
 240   }
 241   // reset R & RT
 242   R[*k-1][*k-1] = 1; R[i][i] = 1;
 243   R[i][*k-1] = 0; R[*k-1][i] = 0;
 244   RT[*k-1][*k-1] = 1; RT[i][i] = 1;
 245   RT[i][*k-1] = 0; RT[*k-1][i] = 0;
 246
 247   for(a=0; a<n; a++)
 248     h[a] = (h[a] + beta*G[*k-1][a])%2;
 249
 250   // update Q & D using Eqs. 52 & 53 with y = beta*G[*k-1][:]
 251   // since h is expressed in F_2^n, in F_2^n, y_i = beta*delta_{i,*k-1}
 252   // (and so the second sum of Eq. 52 is always over zero terms)
 253   *Q = (*Q + ((*D)[*k-1])*beta)%8;
 254
 255   for(a=0; a<*k; a++) {
 256     (*D)[a] = ((*D)[a] + (*J)[a][*k-1]*beta)%8;
 257   }
 258
 259   // remove the kth row & column from J
 260   // remove the ith element from D
 261
 262   // we only need tmpMatrix to be (k-1)x(k-1) and it is kxk
 263   for(a=0; a<(*k-1); a++) {
 264     //tmpMatrix[a] = calloc(*k-1, sizeof(int));  // no need!
 265     for(b=0; b<(*k-1); b++)
 266       tmpMatrix[a][b] = (*J)[a][b];
 267   }
 268   deallocate_mem(J, *k);
 269   *J = calloc(*k-1, sizeof(int*));
 270   for(a=0; a<(*k-1); a++) {
 271     (*J)[a] = calloc(*k-1, sizeof(int));
 272     for(b=0; b<(*k-1); b++)
 273       (*J)[a][b] = tmpMatrix[a][b];
 274   }
 275   deallocate_mem(&tmpMatrix, *k);
 276
 277   tmpVec = calloc(*k-1, sizeof(int));
 278   for(a=0; a<(*k-1); a++)
 279     tmpVec[a] = (*D)[a];
 280   free(*D);
 281   *D = calloc(*k-1, sizeof(int));
 282   memcpy(*D, tmpVec, sizeof(int)*(*k-1));
 283   free(tmpVec);
 284
 285   deallocate_mem(&R, *k);
 286   deallocate_mem(&RT, *k);
 287
 288   // decrement 'k'
 289   *k = *k - 1;
 290
 291   freeList(&setS);
 292
 293   return 'U'; // SUCCESS
 294
 295 }
 296
 297
 298 void append(struct Node** head_ref, int new_data)
 299 {
 300     struct Node* new_node = (struct Node*) malloc(sizeof(struct Node));
 301     struct Node *last = *head_ref;
 302
 303     new_node->data = new_data;
 304     new_node->next = NULL;
 305
 306     // if the linked list is empty, then make the new node as head
 307     if (*head_ref == NULL) {
 308       *head_ref = new_node;
 309       return;
 310     }
 311
 312     // else traverse till the last node
 313     while (last->next != NULL) {
 314         last = last->next;
 315     }
 316     last->next = new_node;
 317
 318     return;
 319 }
 320
 321 /*void freeList(struct Node** head_ref)
 322 {
 323
 324   struct Node *second_last_walker = *head_ref;
 325   struct Node *walker = *head_ref;
 326
 327     // else traverse till the last node
 328     if(walker != NULL) {
 329       while(walker->next != NULL) {
 330         while(walker->next != NULL) {
 331           //printf("!%d\n", walker->data);
 332           second_last_walker = walker;
 333           walker = walker->next;
 334           //printf("%d\n", walker->data);
 335           //printf("!%d\n", second_last_walker->data);
 336         }
 337         free(walker);
 338         second_last_walker->next = NULL;
 339         walker = *head_ref;
 340         //printf("!!%d\n", second_last_walker->data);
 341       }
 342     }
 343     free(walker);
 344
 345     return;
 346     }*/
 347
 348 void freeList(struct Node** head_ref)
 349 {
 350
 351   struct Node *walker = *head_ref;
 352   struct Node *walker2;
 353
 354     // else traverse till the last node
 355     if(walker != NULL) {
 356       walker2 = walker->next;
 357       while(walker2 != NULL) {
 358         free(walker);
 359         walker = walker2;
 360         walker2 = walker->next;
 361       }
 362     }
 363     free(walker);
 364
 365     return;
 366 }
 367
 368 void printLinkedList(struct Node *node)
 369 {
 370   if(node == NULL)
 371     printf("NULL\n");
 372   else {
 373     while (node != NULL)
 374       {
 375         printf(" %d ", node->data);
 376         node = node->next;
 377       }
 378   }
 379   printf("\n");
 380 }