1 /*
2 * minimum-weight-gf3.c
3 *
4 * Minimum Hamming weight computation for codes over GF(3)
5 * The core of computation is in this src file
6 *
7 * NOTES:
8 * � The GF(3) vectors are represented as two packed long
9 * integers, i.e. GF3_VEC structure in minimum-weight-gf3.h.
10 * Each of the long integer is either 32-bit or 64-bit,
11 * depending on the machine architecture, see config.h
12 *
13 * � This packing of GF(3) vectors and their vector manipulations
14 * follow the description given in the paper by
15 *
16 * K. Harrison, D. Page and N. P. Smart, "Software implementation
17 * of finite fields of characteristic three, for use in
18 * pairing-based cryptosystems", London Mathematical Society
19 * Journal of Computation and Mathematics, vol. 5, pp.181--193,
20 * 2002
21 *
22 * � To efficiently generate combination patterns, revolving door
23 * algorithm is used here. This algorithm has the property that
24 * in going to the next combination pattern, there is only one
25 * position that is exchanged, i.e. one 'In' and one 'Out'. The
26 * implementation here follows the third algorithm in the paper
27 * by
28 *
29 * Clement W.H. Lam and Leonard H. Soicher, "Three new combination
30 * algorithms with minimal change property", Communications
31 * of the ACM, vol. 25, no. 8, August 1982
32 *
33 * In the current version of the code, the generation of combination
34 * patterns does not fully exploit the the revolving door property
35 * yet--further (speed) optimisation is possible.
36 *
37 * � To generate all n-tupple of non zero elements of GF(3), Bitner's
38 * algorithm is used. See Algorithm L (Loopless Gray binary
39 * generation) in Knuth's book.
40 *
41 * � The code here is not yet optmised, your contributions in this
42 * respect are greatly appreciated.
43 *
44 * Version Log:
45 * 0.1 18 March 2008 (first released to public -- GUAVA 3.3)
46 *
47 * CJ, Tjhai
48 * email: ctjhai@plymouth.ac.uk
49 * Homepage: www.plymouth.ac.uk/staff/ctjhai
50 *
51 */
52
53 #include <stdio.h>
54 #include <stdlib.h>
55 #include <string.h>
56 #include <math.h>
57 #include "minimum-weight-gf3.h"
58 #include "popcount.h"
59 #include "config.h"
60 #include "types.h"
61
62 /*--------- Some useful macros for updating upper-bound ---------*/
63 #define ROUND_0_MOD_3(d) (((int)ceil((double)d/3.0))*3) /* for self-orthogonal ternary code */
64
65 /*-------------------- Function prototypes ----------------------*/
66 void init_gf3(GF3 *gf);
67 void clear_gf3(GF3 *gf);
68 GF3_VEC **to_packed_integer_gf3(MATRIX *G, int nBlocks);
69 int gf3_gauss_jordan(GF3 *gf, MATRIX *M, int start);
70 void mat_packed_print_gf3(GF3_VEC **M, int nBlocks, int dim, int length);
71 void update_info_gf3(PACKED_MATRIX_GF3 *G, INFO *info);
72 void cyclic_update_info_gf3(PACKED_MATRIX_GF3 *G, INFO *info);
73 void __mindist_gf3(PACKED_MATRIX_GF3 *G, INFO *info);
74 void __cyclic_mindist_gf3(PACKED_MATRIX_GF3 *G, INFO *info);
75 /*---------------------------------------------------------------*/
76
77 /* Return the minimum weight of a general ternary linear code */
mindist_gf3(MATRIX G,mod_t m_mod,int lower_bound)78 int mindist_gf3(MATRIX G, mod_t m_mod, int lower_bound) {
79 int i, j, l, rank, pos;
80 double steps;
81 unsigned int sum=0;
82 PACKED_MATRIX_GF3 packedG;
83 INFO info;
84 GF3 gf;
85
86 /* Initialise GF(3) */
87 init_gf3(&gf);
88
89 /*
90 * Knowing G, now build an array of generator matrices with
91 * disjoint information sets. These matrices are stored in
92 * packed integer format
93 *
94 */
95 j = (G.cols - G.rows) + 1; /* maximum possible number of matrices */
96 packedG.mat = (GF3_VEC ***)malloc(j * sizeof(GF3_VEC **));
97 packedG.rank= (unsigned short *)malloc(j * sizeof(unsigned short));
98 packedG.dimension = G.rows;
99 packedG.block_length = G.cols;
100 packedG.nMatrices = packedG.nFullRankMatrices = 0;
101 packedG.nBlocks = (!(G.cols % BITSIZE)) ? (G.cols / BITSIZE) : (G.cols / BITSIZE)+1;
102 for (i=0,pos=0;;i++) {
103 if (pos >= G.cols) break;
104 rank = gf3_gauss_jordan(&gf, &G, pos); /* what is the rank at column 'pos'? */
105 if (!rank) {
106 pos++;
107 continue;
108 }
109
110 /* Convert the reduced echelon matrix into packed integer format */
111 packedG.mat[i] = to_packed_integer_gf3(&G, packedG.nBlocks);
112
113 packedG.rank[i] = rank;
114 pos += rank;
115 packedG.nMatrices++;
116 if (rank == G.rows) packedG.nFullRankMatrices++;
117 }
118
119 /*---------------------- core computation starts here -------------------------*/
120
121 printf("[%d,%d] linear code over GF(3) - minimum weight evaluation\n", G.cols, G.rows);
122 printf("Known lower-bound: %d\n", lower_bound);
123 printf("There are %d generator matrices, ranks : ", packedG.nMatrices);
124 for (i=0; i<packedG.nMatrices; i++) printf("%d ", packedG.rank[i]); printf("\n");
125 info.weight_constraint = m_mod;
126 if (m_mod == C_0MOD3)
127 printf("The weight of the minimum weight codeword satisfies 0 mod 3 congruence\n");
128 fflush(stdout);
129
130 /* Obtain the first lower-bound */
131 info.known_lower_bound = lower_bound;
132 info.lower_bound = 2 * packedG.nFullRankMatrices;
133 for (i=packedG.nFullRankMatrices; i<packedG.nMatrices; i++) {
134 j = 2 - (G.rows - packedG.rank[i]);
135 if (j>0) info.lower_bound += j;
136 }
137
138 /* Information weight 1 */
139 printf("Enumerating codewords with information weight 1 (w=1)\n"); fflush(stdout);
140 steps = 0; info.upper_bound = G.cols; i = packedG.nMatrices;
141 do { --i;
142 j = packedG.dimension;
143 do { --j;
144 steps++;
145 sum = 0; l = packedG.nBlocks;
146 do { --l;
147 sum += (popcount(packedG.mat[i][j][l].v1) + popcount(packedG.mat[i][j][l].v2));
148 } while (l);
149 if (sum < info.upper_bound) {
150 info.upper_bound = sum;
151 printf(" Found new minimum weight %d\n", info.upper_bound); fflush(stdout);
152 if (info.upper_bound <= info.known_lower_bound) break;
153 }
154 } while (j);
155 if (sum <= info.known_lower_bound) {
156 printf("The known lower-bound is met, enumeration completed.\n");
157 info.upper_bound = sum;
158 goto completed;
159 }
160 } while(i);
161
162 /* Update progress information */
163 info.info_weight_completed = 1;
164 update_info_gf3(&packedG, &info);
165
166 printf("Number of matrices required for codeword enumeration %d\n", info.matrices_req);
167 printf("Completed w= 1, %.0lf codewords enumerated, lower-bound %d, upper-bound %d\n",
168 steps, info.displayed_lower_bound, info.upper_bound); fflush(stdout);
169 if (info.displayed_lower_bound >= info.upper_bound) goto completed;
170 printf("Termination expected with information weight %d at matrix %d\n",
171 info.max_info_weight_req, info.ith_matrix);
172 printf("-----------------------------------------------------------------------------\n");
173
174 /* The following function is executed for higher minimum weight */
175 __mindist_gf3(&packedG, &info);
176
177 /*------------------------ end of core computation ----------------------------*/
178
179 completed:
180 /* Deallocating memory */
181 for (i=0; i<packedG.nMatrices; i++) {
182 for (j=0; j<G.rows; j++) free(packedG.mat[i][j]);
183 free(packedG.mat[i]);
184 }
185 free(packedG.mat);
186 clear_gf3(&gf);
187
188 return info.upper_bound;
189 }
190
191 /* Return the minimum weight of a ternary cyclic code */
cyclic_mindist_gf3(MATRIX G,mod_t m_mod,int lower_bound)192 int cyclic_mindist_gf3(MATRIX G, mod_t m_mod, int lower_bound) {
193 int j, l, rank;
194 unsigned int sum=0;
195 PACKED_MATRIX_GF3 packedG;
196 INFO info;
197 GF3 gf;
198
199 /* Initialise GF(3) */
200 init_gf3(&gf);
201
202 /*
203 * For cyclic code, only one information set is required, so
204 * only one generator matrix (in packed integer format) is
205 * needed.
206 *
207 */
208 j = 1; /* only one generator matrix */
209 packedG.mat = (GF3_VEC ***)malloc(j * sizeof(GF3_VEC **));
210 packedG.rank= (unsigned short *)malloc(j * sizeof(unsigned short));
211 packedG.dimension = G.rows;
212 packedG.block_length = G.cols;
213 packedG.nMatrices = j;
214 packedG.nFullRankMatrices = 1 + ((G.cols - G.rows)/G.rows);
215 packedG.nBlocks = (!(G.cols % BITSIZE)) ? (G.cols / BITSIZE) : (G.cols / BITSIZE)+1;
216 rank = gf3_gauss_jordan(&gf, &G, 0);
217 if (!rank) {
218 fprintf(stderr, "ERROR: rank is 0\n\n"); return -1;
219 }
220
221 /* Convert the reduced echelon matrix into packed integer format */
222 packedG.mat[0] = to_packed_integer_gf3(&G, packedG.nBlocks);
223
224 packedG.rank[0] = rank;
225 if (rank != G.rows) {
226 fprintf(stderr, "ERROR: full rank generator matrix cannot be obtained.\n");
227 fprintf(stderr, " Please provide a full-rank generator matrix as an input.\n\n");
228 return -1;
229 }
230
231 /*---------------------- core computation starts here -------------------------*/
232
233 printf("[%d,%d] cyclic code over GF(3) - minimum weight evaluation\n", G.cols, G.rows);
234 printf("Known lower-bound: %d\n", lower_bound);
235 info.weight_constraint = m_mod;
236 if (m_mod == C_0MOD3)
237 printf("The weight of the minimum weight codeword satisfies 0 mod 3 congruence\n");
238 fflush(stdout);
239 info.matrices_req = 1;
240
241 /* Obtain the first lower-bound */
242 info.known_lower_bound = lower_bound;
243 info.lower_bound = (unsigned int)ceil(((double)G.cols/(double)G.rows)*2.0);
244
245 /* Information weight 1 */
246 printf("Enumerating codewords with information weight 1 (w=1)\n"); fflush(stdout);
247 info.upper_bound = G.cols;
248 j = packedG.dimension;
249 do { --j;
250 sum = 0; l = packedG.nBlocks;
251 do { --l;
252 sum += (popcount(packedG.mat[0][j][l].v1) + popcount(packedG.mat[0][j][l].v2));
253 } while (l);
254 if (sum < info.upper_bound) {
255 info.upper_bound = sum;
256 printf(" Found new minimum weight %d\n", info.upper_bound); fflush(stdout);
257 if (info.upper_bound <= info.known_lower_bound) break;
258 }
259 } while (j);
260 if (sum <= info.known_lower_bound) {
261 printf("The known lower-bound is met, enumeration completed.\n");
262 info.upper_bound = sum;
263 goto completed;
264 }
265
266 /* Update progress information */
267 info.info_weight_completed = 1;
268 cyclic_update_info_gf3(&packedG, &info);
269
270 printf("Number of matrices required for codeword enumeration %d\n", info.matrices_req);
271 printf("Completed w= 1, %d codewords enumerated, lower-bound %d, upper-bound %d\n",
272 packedG.dimension, info.displayed_lower_bound, info.upper_bound); fflush(stdout);
273 if (info.displayed_lower_bound >= info.upper_bound) goto completed;
274 printf("Termination expected with information weight %d\n", info.max_info_weight_req);
275 printf("-----------------------------------------------------------------------------\n");
276
277 /* The following function is executed for higher minimum weight */
278 __cyclic_mindist_gf3(&packedG, &info);
279
280 /*------------------------ end of core computation ----------------------------*/
281
282 completed:
283 /* Deallocating memory */
284 for (j=0; j<G.rows; j++) free(packedG.mat[0][j]);
285 free(packedG.mat[0]);
286 free(packedG.mat);
287 clear_gf3(&gf);
288
289 return info.upper_bound;
290 }
291
292 /* Minimum weight enumeration function for information weight >= 2 */
__mindist_gf3(PACKED_MATRIX_GF3 * G,INFO * info)293 void __mindist_gf3(PACKED_MATRIX_GF3 *G, INFO *info) {
294 #define EVALUATE_MINIMUM_WEIGHT_GF3 {\
295 memset(a, 0, m * sizeof(short));\
296 f[m] = m; j = m; do { --j; f[j] = j; } while (j);\
297 while (1) {\
298 l=info->matrices_req; do { --l;\
299 steps++;\
300 i = G->nBlocks; w = 0;\
301 do { --i;\
302 /* FIXME: (optimisation) */\
303 /* There must a clever (faster) way of doing this */\
304 /* by exploiting the revolving door property */\
305 c[l][i] = G2[0][l][v[1]][i];\
306 j = m; do { --j;\
307 add_gf3_vec(c[l][i], G2[a[j]][l][v[j+2]][i], c[l][i]);\
308 } while (j);\
309 w += (popcount(c[l][i].v1) + popcount(c[l][i].v2));\
310 } while (i);\
311 if (w < info->upper_bound) {\
312 info->upper_bound = w;\
313 printf(" Found new minimum weight %d\n", w); fflush(stdout);\
314 if (w <= info->known_lower_bound) goto lower_bound_met;\
315 }\
316 } while (l);\
317 j = f[0]; f[0] = 0;\
318 if (!(j-m)) break;\
319 else { f[j] = f[j+1]; f[j+1] = j+1; }\
320 a[j] = 1 - a[j];\
321 }\
322 }
323 double steps;
324 short m, iw, top, nMat;
325 register short i, j, l, w;
326 register GF3_VEC **c;
327 short *v, *s, *a, *f;
328 packed_t t1, t2;
329
330 /* G2 contains a set of generator matrices of G->nFullRankMatrices disjoint */
331 /* information sets. For each information set, G2 contains (q-1) matrices, */
332 /* where the identity element of each is non zeros of GF(q) */
333 GF3_VEC ****G2;
334 nMat = info->matrices_req;
335 G2 = (GF3_VEC ****) malloc(2 * sizeof(GF3_VEC ***));
336 G2[0] = G->mat; /* the first one (identity = 1) points to 'G->mat' */
337 G2[1] = (GF3_VEC ***) malloc(nMat * sizeof(GF3_VEC **));
338 for (i=0; i<nMat; i++) {
339 G2[1][i] = (GF3_VEC **) malloc(G->dimension * sizeof(GF3_VEC *));
340 for (j=0; j<G->dimension; j++) {
341 G2[1][i][j] = (GF3_VEC *) malloc(G->nBlocks * sizeof(GF3_VEC));
342 for (l=0; l<G->nBlocks; l++) {
343 G2[1][i][j][l].v1 = G2[0][i][j][l].v2;
344 G2[1][i][j][l].v2 = G2[0][i][j][l].v1;
345 }
346 }
347 }
348
349 w = top = 0;
350 c = (GF3_VEC **)malloc(nMat * sizeof(GF3_VEC*));
351 l = nMat; do { --l;
352 c[l] = (GF3_VEC *)malloc(G->nBlocks * sizeof(GF3_VEC));
353 } while (l);
354 for (iw=2;;iw++) { /* starting from information weight 2 */
355 if (iw == info->max_info_weight_req) info->matrices_req = info->ith_matrix;
356 printf("Enumerating codewords with information weight %d (w=%d) using %d matrices\n",
357 iw, iw, info->matrices_req);
358 fflush(stdout);
359 steps = 0;
360
361 m = iw - 1;
362 v = (short *) malloc((iw+2)*sizeof(short));
363 s = (short *) malloc((iw+2)*sizeof(short));
364 a = (short *) malloc(iw *sizeof(short));
365 f = (short *) malloc((iw+1)*sizeof(short));
366
367 if ( !(iw & 1) ) {
368 v[iw + 1] = G->dimension; v[iw] = iw - 1;
369 if (iw < G->dimension) top = iw;
370 } else {
371 v[iw] = G->dimension - 1;
372 if (iw < G->dimension) top = iw - 1;
373 }
374 v[1] = s[iw] = 0;
375 for (i=2; i<iw; ++i) { v[i] = i-1; s[i] = i+1; }
376
377 /* first combination */
378 EVALUATE_MINIMUM_WEIGHT_GF3;
379
380 /* main loop to generate all other combinations */
381 while (top != 0) {
382 if (top == 2) { /* special handling far v[1] and v[2] */
383 top = s[2]; s[2] = 3;
384 do {
385 v[1] = v[2]; ++v[2];
386 EVALUATE_MINIMUM_WEIGHT_GF3;
387 do {
388 --v[1];
389 EVALUATE_MINIMUM_WEIGHT_GF3;
390 } while (v[1]);
391 } while ( !(v[2] == v[3] - 1) );
392 } else {
393 if (top & 0x1) {
394 --v[top];
395 if (v[top] > top - 1) {
396 top = top - 1; v[top] = top - 1;
397 } else {
398 v[top-1] = top - 2; i = top;
399 top = s[top]; s[i] = i + 1;
400 }
401 } else {
402 v[top-1] = v[top]; ++v[top];
403 if (v[top] == v[top+1] - 1) {
404 s[top-1] = s[top]; s[top] = top + 1;
405 }
406 top -= 2;
407 }
408 EVALUATE_MINIMUM_WEIGHT_GF3;
409 }
410 }
411
412 /* Update progress information */
413 info->info_weight_completed = iw;
414 update_info_gf3(G, info);
415
416 printf("Completed w=%2d, %.0lf codewords enumerated, lower-bound %d, upper-bound %d\n",
417 iw, steps, info->displayed_lower_bound, info->upper_bound); fflush(stdout);
418 lower_bound_met:
419 free(a);
420 free(f);
421 free(v);
422 free(s);
423 if (w <= info->known_lower_bound) {
424 printf("The known lower-bound is met, enumeration completed.\n");
425 printf("-----------------------------------------------------------------------------\n");
426 info->upper_bound = w;
427 for (l=0; l<nMat; l++) {
428 for (i=0; i<G->dimension; i++) free(G2[1][l][i]);
429 free(G2[1][l]);
430 free(c[l]);
431 }
432 free(G2[1]); free(G2);
433 return;
434 }
435 if (info->displayed_lower_bound >= info->upper_bound) break;
436 printf("Termination expected with information weight %d at matrix %d\n",
437 info->max_info_weight_req, info->ith_matrix);
438 printf("-----------------------------------------------------------------------------\n");
439 }
440 printf("-----------------------------------------------------------------------------\n");
441 for (l=0; l<nMat; l++) {
442 for (i=0; i<G->dimension; i++) free(G2[1][l][i]);
443 free(c[l]); free(G2[1][l]);
444 }
445 free(G2[1]); free(G2);
446 free(c);
447 }
448
449 /* Cyclic code: minimum weight enumeration function for information weight >= 2 */
__cyclic_mindist_gf3(PACKED_MATRIX_GF3 * G,INFO * info)450 void __cyclic_mindist_gf3(PACKED_MATRIX_GF3 *G, INFO *info) {
451 #define EVALUATE_CYCLIC_MINIMUM_WEIGHT_GF3 {\
452 memset(a, 0, m * sizeof(short));\
453 f[m] = m; j = m; do { --j; f[j] = j; } while (j);\
454 while (1) {\
455 steps++;\
456 i = G->nBlocks; w = 0;\
457 do { --i;\
458 /* FIXME: (optimisation) */\
459 /* There must a clever (faster) way of doing this */\
460 /* by exploiting the revolving door property */\
461 c[i] = G2[0][0][v[1]][i];\
462 j = m; do { --j; add_gf3_vec(c[i], G2[a[j]][0][v[j+2]][i], c[i]); } while (j);\
463 w += (popcount(c[i].v1) + popcount(c[i].v2));\
464 } while (i);\
465 if (w < info->upper_bound) {\
466 info->upper_bound = w;\
467 printf(" Found new minimum weight %d\n", w); fflush(stdout);\
468 if (w <= info->known_lower_bound) goto lower_bound_met;\
469 }\
470 j = f[0]; f[0] = 0;\
471 if (!(j-m)) break;\
472 else { f[j] = f[j+1]; f[j+1] = j+1; }\
473 a[j] = 1 - a[j];\
474 }\
475 }
476 double steps;
477 short m, iw, top;
478 register short i, j, l, w;
479 register GF3_VEC *c;
480 packed_t t1, t2;
481 short *v, *s, *a, *f;
482
483 /* G2 contains a set of generator matrices of G->nFullRankMatrices disjoint */
484 /* information sets. For each information set, G2 contains (q-1) matrices, */
485 /* where the identity element of each is non zeros of GF(q) */
486 GF3_VEC ****G2;
487 G2 = (GF3_VEC ****) malloc(2 * sizeof(GF3_VEC ***));
488 G2[0] = G->mat; /* the first one (identity = 1) points to 'G->mat' */
489 G2[1] = (GF3_VEC ***) malloc(1 * sizeof(GF3_VEC **));
490 G2[1][0] = (GF3_VEC **) malloc(G->dimension * sizeof(GF3_VEC *));
491 for (j=0; j<G->dimension; j++) {
492 G2[1][0][j] = (GF3_VEC *) malloc(G->nBlocks * sizeof(GF3_VEC));
493 for (l=0; l<G->nBlocks; l++) {
494 G2[1][0][j][l].v1 = G2[0][0][j][l].v2;
495 G2[1][0][j][l].v2 = G2[0][0][j][l].v1;
496 }
497 }
498
499 w = top = 0;
500 c = (GF3_VEC *)malloc(G->nBlocks * sizeof(GF3_VEC));
501 for (iw=2;;iw++) { /* starting from information weight 2 */
502 printf("Enumerating codewords with information weight %d (w=%d) using 1 matrix\n",
503 iw, iw);
504 fflush(stdout);
505 steps = 0;
506
507 m = iw - 1;
508 v = (short *) malloc((iw+2)*sizeof(short));
509 s = (short *) malloc((iw+2)*sizeof(short));
510 a = (short *) malloc(iw *sizeof(short));
511 f = (short *) malloc((iw+1)*sizeof(short));
512
513 if ( !(iw & 1) ) {
514 v[iw + 1] = G->dimension; v[iw] = iw - 1;
515 if (iw < G->dimension) top = iw;
516 } else {
517 v[iw] = G->dimension - 1;
518 if (iw < G->dimension) top = iw - 1;
519 }
520 v[1] = s[iw] = 0;
521 for (i=2; i<iw; ++i) { v[i] = i-1; s[i] = i+1; }
522
523 /* first combination */
524 EVALUATE_CYCLIC_MINIMUM_WEIGHT_GF3;
525
526 /* main loop to generate all other combinations */
527 while (top != 0) {
528 if (top == 2) { /* special handling far v[1] and v[2] */
529 top = s[2]; s[2] = 3;
530 do {
531 v[1] = v[2]; ++v[2];
532 EVALUATE_CYCLIC_MINIMUM_WEIGHT_GF3;
533 do {
534 --v[1];
535 EVALUATE_CYCLIC_MINIMUM_WEIGHT_GF3;
536 } while (v[1]);
537 } while ( !(v[2] == v[3] - 1) );
538 } else {
539 if (top & 0x1) {
540 --v[top];
541 if (v[top] > top - 1) {
542 top = top - 1; v[top] = top - 1;
543 } else {
544 v[top-1] = top - 2; i = top;
545 top = s[top]; s[i] = i + 1;
546 }
547 } else {
548 v[top-1] = v[top]; ++v[top];
549 if (v[top] == v[top+1] - 1) {
550 s[top-1] = s[top]; s[top] = top + 1;
551 }
552 top -= 2;
553 }
554 EVALUATE_CYCLIC_MINIMUM_WEIGHT_GF3;
555 }
556 }
557
558 /* Update progress information */
559 info->info_weight_completed = iw;
560 cyclic_update_info_gf3(G, info);
561
562 printf("Completed w=%2d, %.0lf codewords enumerated, lower-bound %d, upper-bound %d\n",
563 iw, steps, info->displayed_lower_bound, info->upper_bound); fflush(stdout);
564 lower_bound_met:
565 free(v);
566 free(s);
567 free(a);
568 free(f);
569 if (w <= info->known_lower_bound) {
570 printf("The known lower-bound is met, enumeration completed.\n");
571 printf("-----------------------------------------------------------------------------\n");
572 info->upper_bound = w;
573 free(c);
574 for (i=0; i<G->dimension; i++) free(G2[1][0][i]);
575 free(G2[1][0]); free(G2[1]); free(G2);
576 return;
577 }
578 if (info->displayed_lower_bound >= info->upper_bound) break;
579 printf("Termination expected with information weight %d\n",
580 info->max_info_weight_req);
581 printf("-----------------------------------------------------------------------------\n");
582 }
583 printf("-----------------------------------------------------------------------------\n");
584 free(c);
585 for (i=0; i<G->dimension; i++) free(G2[1][0][i]);
586 free(G2[1][0]); free(G2[1]); free(G2);
587 }
588
589 /* Update INFO structure, which contains various information on the */
590 /* current progress of enumeration. */
update_info_gf3(PACKED_MATRIX_GF3 * G,INFO * info)591 void update_info_gf3(PACKED_MATRIX_GF3 *G, INFO *info) {
592 bool end;
593 short i, j, l, d, w;
594
595 /* Estimate the number of matrices required to complete enumeration */
596 info->matrices_req = 0;
597 for (w=0,i=1;;i++) {
598 w = (i+1) * G->nFullRankMatrices;
599 info->matrices_req = G->nFullRankMatrices;
600 for (j=G->nFullRankMatrices; j<G->nMatrices; j++) {
601 l = i - (G->dimension - G->rank[j]) + 1;
602 if (l > 0) { w += l; info->matrices_req++; }
603 }
604 end = false;
605 switch (info->weight_constraint) {
606 case C_0MOD3: if (ROUND_0_MOD_3(w) >= info->upper_bound) end = true; break;
607 default: if (w >= info->upper_bound) end = true;
608 }
609 if (end) break;
610 }
611
612 /* Estimate the minimum distance lower bound after current enumeration */
613 if (info->info_weight_completed == 1) /* required for safety */
614 info->max_info_weight_req = G->dimension;
615 if (info->info_weight_completed < info->max_info_weight_req) {
616 info->lower_bound = (info->info_weight_completed + 1) * G->nFullRankMatrices;
617 for (j=G->nFullRankMatrices; j<G->nMatrices; j++) {
618 l = info->info_weight_completed - (G->dimension - G->rank[j]) + 1;
619 if (l > 0) info->lower_bound += l;
620 }
621 } else { /* last information weight in enumeration */
622 for (j=0; j<info->ith_matrix; j++) {
623 if (j < G->nFullRankMatrices) info->lower_bound++;
624 else {
625 l = i - (G->dimension - G->rank[j]) + 1;
626 if (l > 0) info->lower_bound++;
627 }
628 }
629 }
630
631 /* lower-bound information to display, this value could be different */
632 /* from the actual lower-bound if the code satisfies some weight */
633 /* constraint, i.e. 0 mod 3 (self-orthogonal ternay code). */
634 switch (info->weight_constraint) {
635 case C_0MOD3: info->displayed_lower_bound = ROUND_0_MOD_3(info->lower_bound); break;
636 default: info->displayed_lower_bound = info->lower_bound;
637 }
638
639 /* Estimate the maximum information weight required */
640 d = 0;
641 w = info->lower_bound;
642 for (i=info->info_weight_completed+1;;i++) {
643 for (j=0; j<G->nMatrices; j++) {
644 if (j < G->nFullRankMatrices)
645 l = 1; /* this is a trick to reduce the number of lines of codes */
646 else
647 l = i - (G->dimension - G->rank[j]) + 1;
648 if (l > 0) {
649 w++;
650 switch (info->weight_constraint) {
651 case C_0MOD3: d = ROUND_0_MOD_3(w); break;
652 default: d = w;
653 }
654 if (d >= info->upper_bound) { info->ith_matrix = j + 1; break; }
655 }
656 }
657 if (d >= info->upper_bound) { info->max_info_weight_req = i; break; }
658 }
659 }
660
661 /* Update INFO structure, which contains various information on the */
662 /* current progress of enumeration. */
cyclic_update_info_gf3(PACKED_MATRIX_GF3 * G,INFO * info)663 void cyclic_update_info_gf3(PACKED_MATRIX_GF3 *G, INFO *info) {
664 short i, d, w;
665
666 /* Estimate the minimum distance lower bound after current enumeration */
667 info->lower_bound = (info->info_weight_completed + 1) * G->nFullRankMatrices;
668 info->lower_bound = (unsigned int)
669 ceil(((double)G->block_length/(double)G->dimension)*((double)info->info_weight_completed+1.0));
670
671 /* lower-bound information to display, this value could be different */
672 /* from the actual lower-bound if the code satisfies some weight */
673 /* constraint, i.e. 0 mod 3 (self-orthogonal ternary code). */
674 switch (info->weight_constraint) {
675 case C_0MOD3: info->displayed_lower_bound = ROUND_0_MOD_3(info->lower_bound); break;
676 default: info->displayed_lower_bound = info->lower_bound;
677 }
678
679 /* Estimate the maximum information weight required */
680 d = 0;
681 for (i=info->info_weight_completed+1;;i++) {
682 w = (unsigned int)
683 ceil(((double)G->block_length/(double)G->dimension)*((double)i+1.0));
684 switch (info->weight_constraint) {
685 case C_0MOD3: d = ROUND_0_MOD_3(w); break;
686 default: d = w;
687 }
688 if (d >= info->upper_bound) { info->max_info_weight_req = i; break; }
689 }
690 }
691
692 /* Initialise GF subtraction and division tables */
init_gf3(GF3 * gf)693 void init_gf3(GF3 *gf) {
694 int i;
695
696 gf->sub = (short **)malloc(3 * sizeof(short *));
697 for (i=0; i<3; i++) gf->sub[i] = (short *)malloc(3*sizeof(short));
698 gf->sub[0][0] = 0; gf->sub[0][1] = 2; gf->sub[0][2] = 1;
699 gf->sub[1][0] = 1; gf->sub[1][1] = 0; gf->sub[1][2] = 2;
700 gf->sub[2][0] = 2; gf->sub[2][1] = 1; gf->sub[2][2] = 0;
701
702 gf->div = (short **)malloc(3 * sizeof(short int *));
703 for (i=0; i<3; i++) gf->div[i] = (short *)malloc(3*sizeof(short));
704 gf->div[0][0] = 3; gf->div[0][1] = 0; gf->div[0][2] = 0;
705 gf->div[1][0] = 3; gf->div[1][1] = 1; gf->div[1][2] = 2;
706 gf->div[2][0] = 3; gf->div[2][1] = 2; gf->div[2][2] = 1;
707 }
708
709 /* Deallocation */
clear_gf3(GF3 * gf)710 void clear_gf3(GF3 *gf) {
711 int i;
712 for (i=0; i<3; i++) {
713 free(gf->sub[i]);
714 free(gf->div[i]);
715 }
716 free(gf->sub);
717 free(gf->div);
718 }
719
720 /* Convert a ternary matrix into packed integer format */
to_packed_integer_gf3(MATRIX * G,int nBlocks)721 GF3_VEC **to_packed_integer_gf3(MATRIX *G, int nBlocks) {
722 int i, j;
723 GF3_VEC **M;
724
725 M = (GF3_VEC **)malloc(G->rows * sizeof(GF3_VEC *));
726 for (i=0; i<G->rows; i++) {
727 M[i] = (GF3_VEC *)malloc(nBlocks * sizeof(GF3_VEC));
728 memset(M[i], 0, nBlocks * sizeof(GF3_VEC));
729 for (j=0; j<nBlocks; j++) M[i][j].v1 = M[i][j].v2 = ZERO;
730 for (j=0; j<G->cols; j++) {
731 switch (G->m[i][j]) {
732 case 1 : M[i][j >> LOG2].v2 |= (ONE << (j & MOD)); break;
733 case 2 : M[i][j >> LOG2].v1 |= (ONE << (j & MOD)); break;
734 default: break;
735 }
736 }
737 }
738 return M;
739 }
740
741 /* Print out a packed GF3 matrix in unpacked format */
mat_packed_print_gf3(GF3_VEC ** M,int nBlocks,int dim,int length)742 void mat_packed_print_gf3(GF3_VEC **M, int nBlocks, int dim, int length) {
743 unsigned int r, i, j;
744 packed_t u1, u2;
745 for (r=0; r<dim; r++) {
746 for (i=0; i<nBlocks-1; i++) {
747 for (j=0; j<BITSIZE; j++) {
748 u1 = M[r][i].v1 & (ONE << j);
749 u2 = M[r][i].v2 & (ONE << j);
750 if ((!u1) && (!u2)) printf("0");
751 else if ((!u1) && u2) printf("1");
752 else if (u1 && !u2) printf("2");
753 else printf("?");
754 }
755 }
756 for (i=nBlocks-1,j=0; j<length-(BITSIZE*(nBlocks-1)); j++) {
757 u1 = M[r][i].v1 & (ONE << j);
758 u2 = M[r][i].v2 & (ONE << j);
759 if ((!u1) && (!u2)) printf("0");
760 else if ((!u1) && u2) printf("1");
761 else if (u1 && !u2) printf("2");
762 else printf("?");
763 }
764 printf("\n");
765 }
766 }
767
768 /* gauss jordan - turn matrix M into reduced-echelon form at position 'start' */
gf3_gauss_jordan(GF3 * gf,MATRIX * M,int start)769 int gf3_gauss_jordan(GF3 *gf, MATRIX *M, int start) {
770 bool found;
771 int i, r, c, pivot, swap;
772 for (c=0; c<M->rows; ++c) { /* diagonalise the element from col=1 to rows */
773 if (c+start >= M->cols) return c;
774 if (!M->m[c][c+start]) { /* pivot is zero, do row swapping */
775 for (found=false,r=c+1; r<M->rows; ++r) { if (M->m[r][c+start]) { found=true; break; } }
776 if (found) {
777 for (i=0; i<M->cols; ++i) { swap=M->m[r][i]; M->m[r][i]=M->m[c][i]; M->m[c][i]=swap; }
778 } else { /* row swapping failed, do column swapping */
779 for (found=false,i=start+c+1; i<M->cols; ++i) {
780 for (r=c; r<M->rows; ++r) { if (M->m[r][i]) { found=true; break; } }
781 if (found) break;
782 }
783 if (!found) return c;
784 for (pivot=i,i=0; i<M->cols; ++i) { swap=M->m[r][i]; M->m[r][i]=M->m[c][i]; M->m[c][i]=swap; }
785 for (r=0; r<M->rows; ++r) {
786 swap=M->m[r][pivot]; M->m[r][pivot]=M->m[r][c+start]; M->m[r][c+start]=swap;
787 }
788 }
789 }
790 if (M->m[c][c+start]>1) { /* ensure the diagonal element is 1 */
791 for (pivot=M->m[c][c+start],i=0; i<M->cols; ++i)
792 M->m[c][i] = divide_gf3(gf, M->m[c][i], pivot);
793 }
794 for (r=0; r<M->rows; ++r) { /* zero the whole rows at column c, except at row c (diagonalisation) */
795 if (r==c) continue;
796 if (M->m[r][c+start]) {
797 pivot = M->m[r][c+start];
798 for (i=0; i<M->cols; ++i)
799 M->m[r][i] = subtract_gf3(gf, divide_gf3(gf, M->m[r][i], pivot), M->m[c][i]);
800 }
801 }
802 }
803 /* ensure that the diagonal element is 1 */
804 for (c=0; c<M->rows; ++c) {
805 if (M->m[c][c+start]>1) {
806 for (pivot=M->m[c][c+start],i=0; i<M->cols; ++i)
807 M->m[c][i] = divide_gf3(gf, M->m[c][i], pivot);
808 }
809 }
810 return c;
811 }
812