--- /dev/null
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <complex.h>
+
+#include "matrix.h"
+#include "shrinkstar.h"
+
+// linked list
+struct Node {
+ int data;
+ struct Node* next;
+};
+
+void printLinkedList(struct Node *node);
+void append(struct Node** head_ref, int new_data);
+void freeList(struct Node** head_ref);
+
+/****************************************************************************
+ * Attempts to shrinks the dimension of the affine space by one dimension.
+ * The affine space is specified by
+ * the first 'k' columns of the 'n'x'n' matrix 'G', translated by 'h',
+ * that have an inner product with 'xi' of 'alpha'. 'GBar' is 'G^(-1)^T',
+ * and ('Q', 'D', 'J') specify the quadratic form on this affine space.
+ * ---
+ * If result is an EMPTY affine space then return value is 'E'.
+ * If result is the SAME affine space then return value is 'A'.
+ * If result is a SUCCESSfully shrunk affine space then return value is 'U'.
+ * ***
+ * Note that D, and J are assumed to be nx1 and nxn even though they are
+ * actually kx1 and kxk (extra elements are ignored).
+ ****************************************************************************/
+
+char shrinkstar(int n, int *k, int *h, int **G, int **GBar, int *xi, int alpha) {
+ // 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.
+
+ int a, b, c, d;
+ int beta;
+
+ struct Node* setS = NULL;
+ struct Node* setWalker = NULL;
+
+ for(a=0; a<*k; a++) {
+ if(dotProductMod(xi, G[a], n, 2) == 1)
+ append(&setS, a);
+ }
+
+ beta = (alpha + dotProductMod(xi, h, n, 2))%2;
+
+ if(setS == NULL) {
+ if(beta == 1) {
+ freeList(&setS);
+ return 'E'; // EMPTY
+ }
+ else {// beta = 0
+ freeList(&setS);
+ return 'A'; // SAME
+ }
+ }
+
+ // 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
+ int* tmpVec; // temporary vector used in swapping G and GBar rows and in Eq. 49 for updating D
+ /*int** tmpMatrix; // temporary matrix used to store intermediate value of J = R times J times R^T
+ tmpMatrix = calloc(*k, sizeof(int*));
+ for(a=0; a<*k; a++)
+ tmpMatrix[a] = calloc(*k, sizeof(int));*/
+
+ int** R; // basis change matrix R and R^T
+ int** RT;
+ R = calloc(*k, sizeof(int*));
+ RT = calloc(*k, sizeof(int*));
+ // initially set it to the k-dimensional identity matrix
+ for(a=0; a<*k; a++) {
+ R[a] = calloc(*k, sizeof(int));
+ R[a][a] = 1;
+ RT[a] = calloc(*k, sizeof(int));
+ RT[a][a] = 1;
+ }
+
+ int i = setS->data; // pick first element in setS to be added to the new space's shift vector
+ // remove from setS
+ if(setS->next != NULL) {
+ setWalker = setS->next;
+ setS->data = setS->next->data;
+ setS->next = setS->next->next;
+ free(setWalker);
+ } else
+ setS = NULL;
+
+ int tmp;
+ setWalker = setS;
+ while(setWalker != NULL) {
+ // update G rows, g
+ for(a=0; a<n; a++)
+ G[setWalker->data][a] = (G[setWalker->data][a] + G[i][a])%2;
+
+ // update D and J
+ /*R[setWalker->data][i] = 1;
+ RT[i][setWalker->data] = 1;
+
+ tmpVec = calloc(*k, sizeof(int));
+ for(c=0;c<*k;c++) {
+ tmp = 0;*/
+ /* tmpVec[c] = 0; */
+ /* for(a=0;a<*k;a++) */
+ /* tmpVec[c] = tmpVec[c] + R[c][a]*((*D)[a]); */
+ /* for(a=0;a<*k;a++) */
+ /* for(b=a+1;b<*k;b++) */
+ /* tmpVec[c] = tmpVec[c] + (*J)[a][b]*R[c][a]*R[c][b]; */
+ /* tmpVec[c] = tmpVec[c]%8; */
+ /* if(c!=i)
+ tmp += R[c][c]*((*D)[c]);
+ tmp += R[c][i]*((*D)[i]);
+ if(i>c)
+ tmp += ((*J)[c][i])*R[c][c]*R[c][i];
+ else if(i<c)
+ tmp += ((*J)[i][c])*R[c][i]*R[c][c];
+ tmpVec[c] = tmp%8;
+ }
+ memcpy(*D, tmpVec, sizeof(int)*(*k));
+ free(tmpVec);
+
+ // Eq. 50 update of J
+ for(c=0;c<*k;c++) {
+ for(d=0;d<*k; d++) {
+ tmp = 0;
+ if(c!=i) {
+ if(d!=i) {
+ tmp += R[c][c]*((*J)[c][d])*R[d][d];
+ tmp += R[c][c]*((*J)[c][i])*R[d][i];
+ tmp += R[c][i]*((*J)[i][d])*R[d][d];
+ }else {
+ tmp += R[c][c]*((*J)[c][i])*R[d][i];
+ }
+ } else {
+ if(d!=i) {
+ tmp += R[c][i]*((*J)[i][d])*R[d][d];
+ }
+ }
+ tmp += R[c][i]*((*J)[i][i])*R[d][i];
+ tmp %= 8;
+ tmpMatrix[c][d] = tmp;
+ }
+ }
+ for(c=0; c<*k; c++)
+ memcpy(((*J)[c]), tmpMatrix[c], sizeof(int)*(*k));*/
+ /* multMatrixMod(*J,RT,tmpMatrix,*k,*k,*k,*k,8); */
+ /* multMatrixMod(R,tmpMatrix,*J,*k,*k,*k,*k,8); */
+
+ if(setWalker->data != i) {
+ R[setWalker->data][i] = 0; // reset R
+ RT[i][setWalker->data] = 0;
+ }
+
+ // update GBar rows, gbar
+ for(a=0; a<n; a++)
+ GBar[i][a] = (GBar[i][a] + GBar[setWalker->data][a])%2;
+
+ setWalker = setWalker->next;
+ }
+
+ // Swap 'i' and 'k'
+ R[*k-1][*k-1] = 0; R[i][i] = 0;
+ R[i][*k-1] = 1; R[*k-1][i] = 1;
+ RT[*k-1][*k-1] = 0; RT[i][i] = 0;
+ RT[i][*k-1] = 1; RT[*k-1][i] = 1;
+
+ tmpVec = G[i];
+ G[i] = G[*k-1];
+ G[*k-1] = tmpVec;
+ tmpVec = GBar[i];
+ GBar[i] = GBar[*k-1];
+ GBar[*k-1] = tmpVec;
+
+ // update D and J
+ //(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.)
+ /* if(i!=(*k-1)) {
+ tmpVec = calloc(*k, sizeof(int));
+ for(c=0;c<*k;c++) {
+ tmp = 0;*/
+ /* tmpVec[c] = 0; */
+ /* for(a=0;a<*k;a++) */
+ /* tmpVec[c] = tmpVec[c] + R[c][a]*((*D)[a]); */
+ /* for(a=0;a<*k;a++) { */
+ /* for(b=a+1;b<*k;b++) */
+ /* tmpVec[c] = tmpVec[c] + (*J)[a][b]*R[c][a]*R[c][b]; */
+ /* } */
+ /* tmpVec[c] = tmpVec[c]%8; */
+ /* tmp += R[c][c]*((*D)[c]); // iff c!=i and c!=(*k-1) this can be non-zero
+ tmp += R[c][*k-1]*((*D)[*k-1]); // iff c=i this can be non-zero
+ tmp += R[c][i]*((*D)[i]); // iff c=*k-1 this can be non-zero
+ */
+ /* if(i>c) */
+ /* tmp += ((*J)[c][i])*R[c][c]*R[c][i]; */
+ /* else if(i<c) */
+ /* tmp += ((*J)[i][c])*R[c][i]*R[c][c]; */
+ /* if((*k-1)>c) */
+ /* tmp += ((*J)[c][*k-1])*R[c][c]*R[c][*k-1]; */
+ /* else if((*k-1)<c) // pretty sure this can't happen */
+ /* tmp += ((*J)[*k-1][c])*R[c][*k-1]*R[c][c]; */
+ /* if((*k-1)>i) */
+ /* tmp += ((*J)[i][*k-1])*R[c][i]*R[c][*k-1]; */
+ /* else if((*k-1)<i) */
+ /* tmp += ((*J)[*k-1][i])*R[c][*k-1]*R[c][i]; */
+ /* tmpVec[c] = tmp%8;
+ }
+ memcpy(*D, tmpVec, sizeof(int)*(*k));
+ free(tmpVec);
+
+ // Eq. 50 update of J
+ for(c=0;c<*k;c++) {
+ for(d=0;d<*k; d++) {
+ tmp = 0;
+ if(c!=(*k-1) && c!=i) {
+ if(d!=(*k-1) && d!=i) {
+ tmp += R[c][c]*((*J)[c][d])*R[d][d];
+ }else {
+ tmp += R[c][c]*((*J)[c][i])*R[d][i];
+ tmp += R[c][c]*((*J)[c][*k-1])*R[d][*k-1];
+ }
+ } else {
+ if(d!=(*k-1) && d !=i) {
+ tmp += R[c][i]*((*J)[i][d])*R[d][d];
+ tmp += R[c][*k-1]*((*J)[*k-1][d])*R[d][d];
+ } else {
+ tmp += R[c][i]*((*J)[i][*k-1])*R[d][*k-1];
+ tmp += R[c][*k-1]*((*J)[*k-1][i])*R[d][i];
+ tmp += R[c][i]*((*J)[i][i])*R[d][i];
+ tmp += R[c][*k-1]*((*J)[*k-1][*k-1])*R[d][*k-1];
+ }
+ }
+ tmpMatrix[c][d] = tmp%8;
+ }
+ }
+ for(c=0; c<*k; c++)
+ memcpy(((*J)[c]), tmpMatrix[c], sizeof(int)*(*k));*/
+ /* multMatrixMod(*J,RT,tmpMatrix,*k,*k,*k,*k,8); */
+ /* multMatrixMod(R,tmpMatrix,*J,*k,*k,*k,*k,8); */
+ /*}*/
+ // reset R & RT
+ R[*k-1][*k-1] = 1; R[i][i] = 1;
+ R[i][*k-1] = 0; R[*k-1][i] = 0;
+ RT[*k-1][*k-1] = 1; RT[i][i] = 1;
+ RT[i][*k-1] = 0; RT[*k-1][i] = 0;
+
+ for(a=0; a<n; a++)
+ h[a] = (h[a] + beta*G[*k-1][a])%2;
+
+ // update Q & D using Eqs. 52 & 53 with y = beta*G[*k-1][:]
+ // since h is expressed in F_2^n, in F_2^n, y_i = beta*delta_{i,*k-1}
+ // (and so the second sum of Eq. 52 is always over zero terms)
+ /**Q = (*Q + ((*D)[*k-1])*beta)%8;
+
+ for(a=0; a<*k; a++) {
+ (*D)[a] = ((*D)[a] + (*J)[a][*k-1]*beta)%8;
+ }*/
+
+ // remove the kth row & column from J
+ // remove the ith element from D
+
+ // we only need tmpMatrix to be (k-1)x(k-1) and it is kxk
+ /*for(a=0; a<(*k-1); a++) {
+ //tmpMatrix[a] = calloc(*k-1, sizeof(int)); // no need!
+ for(b=0; b<(*k-1); b++)
+ tmpMatrix[a][b] = (*J)[a][b];
+ }
+ deallocate_mem(J, *k);
+ *J = calloc(*k-1, sizeof(int*));
+ for(a=0; a<(*k-1); a++) {
+ (*J)[a] = calloc(*k-1, sizeof(int));
+ for(b=0; b<(*k-1); b++)
+ (*J)[a][b] = tmpMatrix[a][b];
+ }
+ deallocate_mem(&tmpMatrix, *k);*/
+
+ /*tmpVec = calloc(*k-1, sizeof(int));
+ for(a=0; a<(*k-1); a++)
+ tmpVec[a] = (*D)[a];
+ free(*D);
+ *D = calloc(*k-1, sizeof(int));
+ memcpy(*D, tmpVec, sizeof(int)*(*k-1));
+ free(tmpVec);*/
+
+ deallocate_mem(&R, *k);
+ deallocate_mem(&RT, *k);
+
+ // decrement 'k'
+ *k = *k - 1;
+
+ freeList(&setS);
+
+ return 'U'; // SUCCESS
+
+}
+
+
+void append(struct Node** head_ref, int new_data)
+{
+ struct Node* new_node = (struct Node*) malloc(sizeof(struct Node));
+ struct Node *last = *head_ref;
+
+ new_node->data = new_data;
+ new_node->next = NULL;
+
+ // if the linked list is empty, then make the new node as head
+ if (*head_ref == NULL) {
+ *head_ref = new_node;
+ return;
+ }
+
+ // else traverse till the last node
+ while (last->next != NULL) {
+ last = last->next;
+ }
+ last->next = new_node;
+
+ return;
+}
+
+/*void freeList(struct Node** head_ref)
+{
+
+ struct Node *second_last_walker = *head_ref;
+ struct Node *walker = *head_ref;
+
+ // else traverse till the last node
+ if(walker != NULL) {
+ while(walker->next != NULL) {
+ while(walker->next != NULL) {
+ //printf("!%d\n", walker->data);
+ second_last_walker = walker;
+ walker = walker->next;
+ //printf("%d\n", walker->data);
+ //printf("!%d\n", second_last_walker->data);
+ }
+ free(walker);
+ second_last_walker->next = NULL;
+ walker = *head_ref;
+ //printf("!!%d\n", second_last_walker->data);
+ }
+ }
+ free(walker);
+
+ return;
+ }*/
+
+void freeList(struct Node** head_ref)
+{
+
+ struct Node *walker = *head_ref;
+ struct Node *walker2;
+
+ // else traverse till the last node
+ if(walker != NULL) {
+ walker2 = walker->next;
+ while(walker2 != NULL) {
+ free(walker);
+ walker = walker2;
+ walker2 = walker->next;
+ }
+ }
+ free(walker);
+
+ return;
+}
+
+void printLinkedList(struct Node *node)
+{
+ if(node == NULL)
+ printf("NULL\n");
+ else {
+ while (node != NULL)
+ {
+ printf(" %d ", node->data);
+ node = node->next;
+ }
+ }
+ printf("\n");
+}
projects / weak_simulation_stab_extent.git / blobdiff