1 #include "sys.h"
2 #include "debug.h"
3 #include <cstring> // memset
4 #include <iostream>
5 #include <inttypes.h>
6 #include <vector>
7 #include <libecc/bitset.h>
8 #include <math.h>
9
10 // This program has been used to determine feedback points.
11 // It is multipurpose: you have to edit the source code a
12 // bit in order to make it do whatever you want.
13 // Compile as g++ -O3 -o rngTestMatrix -I../include rngTestMatrix.cc -L../.libs -Wl,-rpath -Wl,/home/carlo/c++/libecc/.libs -lecc
14
15 // Define this to use a faster algorithm that is MUCH harder to understand :p
16 #define COMPRESSED 1
17
18 // This program calculates the period of a
19 // shift register random number generator
20 // with arbitrary feedback points.
21 //
22 // Let the internal state of the RNG consist of
23 // 'b' buckets:
24 //
25 // bucket: 0 1 2 3 4 ... b-1
26 // value : n0 n1 n2 n3 n4 ... n(b-1)
27 //
28 // where b equals the number of buckets and each n represents
29 // an integer in the range 0 <= n < m, with m = 2^number_of_bits_per_bucket.
30 // Because of the memory constrains and because the period of
31 // the RNG doesn't increase much anyway, the rest of the program
32 // assumes that each bucket is only 1 bit.
33
34 // 2^number_of_bits - 1 must be a prime, which
35 // are called Mersenne primes.
36
37 // See http://www.utm.edu/research/primes/mersenne/
38 unsigned int const mersenne_prime_powers[] =
39 { 2, 3, 5, 7, 13, 17, 19, 31, 61, 89, 107, 127, 521, 607, 1279, 2203, 2281,
40 3217, 4253, 4423, 9689, 9941, 11213, 19937, 21701, 23209, 44497, 86243, 110503, 132049, 216091,
41 756839, 859433, 1257787, 1398269, 2976221, 3021377, /* ...???... */ 6972593, /* ...???...*/ 13466917 };
42
43 // which limits the possible values of number_of_bits.
44
45 #ifdef M
46 unsigned int const number_of_bits = M;
47 #else
48 unsigned int const number_of_bits = 63;
49 #endif
50
51 // The state of the RNG is changed so that we get
52 // the following state:
53 //
54 // bucket: 0 1 2 3 ... b-2 b-1
55 // value : n1 n2 n3 n4 ... n(b-1) x
56 //
57 // where x = n0 + n(b-f0) + n(b-f1) + n(b-f2) + ... + n(b-f(fp-1)) (MOD m)
58 // fp are the number of feedback points, and f0, f1, ... f(fp-1)
59 // are the feedback points.
60 //
61 // x is also passed as the output of the RNG.
62
63 unsigned int const number_of_feedback_points = 1;
64
65 // The feedback points must be small prime numbers.
66 std::vector<unsigned int> primes;
67 unsigned int const max_prime = number_of_bits;
68 void init_primes(void);
69
70 unsigned int feedback_points[number_of_feedback_points];
71 int feedback_point_indexes[number_of_feedback_points];
72 void init_feedback_points(void);
73
74 // If we consider the state of a RNG to be a vector.
75 // then we can calculate the next state of the RNG
76 // by multiplying the state with a matrix.
77
78 struct Matrix {
79 // columns rows
80 libecc::bitset<number_of_bits> M_matrix[number_of_bits];
clearMatrix81 void clear(void) { for (int row = 0; row < number_of_bits; ++row) M_matrix[row].reset(); }
MatrixMatrix82 Matrix(void) { clear(); }
83 void square(Matrix const& m);
84 void mult(Matrix const& m1, Matrix const& m2);
85 friend bool operator==(Matrix const& m1, Matrix const& m2);
86 friend std::ostream& operator<<(std::ostream& os, Matrix const& m);
87 };
88
init_m0(Matrix * m0)89 void init_m0(Matrix* m0)
90 {
91 // m0 is the matrix describing a shift of one.
92 // This means that:
93 //
94 // 0 1 0 0 0 0 0 ... 0 n0 n1
95 // 0 0 1 0 0 0 0 ... 0 n1 n2
96 // 0 0 0 1 0 0 0 ... 0 n2 n3
97 // 0 0 0 0 1 0 0 ... 0 n3 = .
98 // . . . .
99 // . . . .
100 // . . . n(b-1)
101 // 1 ... 1 . 1 ... n(b-1) n0 + n(b-f0) + n(b-f1) + n(b-f2) + ... + n(b-f(fp-1))
102 //
103 // Where the 1's on the bottom row are in column 0, b-f0, b-f1, etc.
104 m0->clear();
105 for (int row = 0; row < number_of_bits - 1; ++row)
106 m0->M_matrix[row].set(row + 1);
107 m0->M_matrix[number_of_bits - 1].set<0>();
108 for (int f = 0; f < number_of_feedback_points; ++f)
109 m0->M_matrix[number_of_bits - 1].set(number_of_bits - feedback_points[f]);
110 }
111
112 // For debugging purposes.
operator <<(std::ostream & os,libecc::bitset<number_of_bits> const * vector)113 std::ostream& operator<<(std::ostream& os, libecc::bitset<number_of_bits> const* vector)
114 {
115 for (int c = 0; c < number_of_bits; ++c)
116 {
117 if (vector->test(c))
118 os << "1 ";
119 else
120 os << "0 ";
121 }
122 return os;
123 }
124
main(void)125 int main(void)
126 {
127 Debug( check_configuration() );
128 // Debug( libcw_do.on() );
129 Debug( dc::notice.on() );
130
131 #if 0
132 // Make sure that 2^number_of_bits - 1 is a prime.
133 int i;
134 for (i = 0; i < sizeof(mersenne_prime_powers) / sizeof(unsigned int); ++i)
135 if (mersenne_prime_powers[i] == number_of_bits)
136 break;
137 assert( i != sizeof(mersenne_prime_powers) / sizeof(unsigned int) );
138 #endif
139
140 init_primes();
141
142 // for(std::vector<unsigned int>::iterator iter = primes.begin(); iter != primes.end(); ++iter)
143 // std::cout << (*iter) << std::endl;
144 // exit(0);
145
146 init_feedback_points();
147
148 #if !COMPRESSED
149 Matrix* m0 = new Matrix;
150 Matrix* m1 = new Matrix;
151 #endif
152
153 for(;;)
154 {
155 std::cout << "Trying " << number_of_bits;
156 for (int f = 0; f < number_of_feedback_points; ++f)
157 std::cout << '/' << feedback_points[f];
158 std::cout << "... " << std::flush;
159
160 #if !COMPRESSED
161 init_m0(m0);
162
163 // std::cout << "\nShift of 1:\n" << *m0 << std::endl;
164
165 #if 0
166 // Print out matrixes for arbitrary powers of M0.
167 Matrix* m2 = new Matrix;
168
169 int shift = 1;
170 int nextshift = 1;
171 *m1 = *m0;
172 for(;;)
173 {
174 if (shift == nextshift)
175 {
176 std::cout << "Shift = " << shift << '\n';
177 std::cout << *m1 << '\n';
178 nextshift *= 2; // Edit this to get the powers you want to see.
179 }
180 m2->mult(*m1, *m0);
181 ++shift;
182 *m1 = *m2;
183 if (shift == nextshift)
184 {
185 std::cout << "(hit return)\n";
186 char tmp[8];
187 std::cin.getline(tmp, 8);
188 }
189 }
190 #endif
191
192 // A shift of 2 is m[1] = m[0] * m[0]
193 // A shift of 4 is m[2] = m[1] * m[1]
194 // A shift of 8 is m[3] = m[2] * m[2]
195 // etc.
196 if (number_of_bits > 521)
197 std::cout << "Calculating... ";
198 for (int n = 1;;)
199 {
200 if (number_of_bits > 521)
201 std::cout << 'M' << n << ' ' << std::flush;
202 m1->square(*m0);
203
204 if (n++ == number_of_bits)
205 {
206 init_m0(m0);
207 if (*m0 == *m1)
208 std::cout << "successful, the period is 2^" << number_of_bits << " - 1.";
209 std::cout << std::endl;
210 break;
211 }
212
213 if (number_of_bits > 521)
214 std::cout << 'M' << n << ' ' << std::flush;
215 m0->square(*m1);
216
217 if (n++ == number_of_bits)
218 {
219 init_m0(m1);
220 if (*m0 == *m1)
221 std::cout << "successful, the period is 2^" << number_of_bits << " - 1.";
222 std::cout << std::endl;
223 break;
224 }
225 }
226 if (number_of_bits > 521)
227 std::cout << '\n';
228 #else // Using compressed matrices.
229
230 // Meaning of M0: See above.
231 //
232 // The period of n/f1/f2/f3.../fx is the same as the period of n/(n-f1)/(n-f2)/(n-f3).../(n-fx).
233 // All our feedback points are < n/2 but we calculate the period of a RNG with the mirrored
234 // values, thus a RNG for which all feedback points are > n/2.
235 //
236 // Example M0
237 //
238 // Consider the period of a 7/3 RNG (n == 7, feedback point at 3).
239 // Its period is the same as that of a 7/4 RNG.
240 //
241 // M0 is then:
242 //
243 // 0 1 0 0 0 0 0
244 // 0 0 1 0 0 0 0
245 // 0 0 0 1 0 0 0
246 // 0 0 0 0 1 0 0
247 // 0 0 0 0 0 1 0
248 // 0 0 0 0 0 0 1
249 // 1 0 0 \ 0 0 0
250 // |______ feedback point at '4' , '\' == '1' but represents a feedback point.
251 //
252 // M0^6 is
253 //
254 // 0 0 0 0 0 0 1
255 // 1 0 0 \ 0 0 0 <-- same feedback point, this row is equal to the last row of M0.
256 // 0 1 0 0 \ 0 0 <-- subsequential rows are shifted to the left.
257 // 0 0 1 0 0 \ 0
258 // 0 0 0 1 0 0 \
259 // \ 0 0 \ 1 0 0 \__ when a diagonal "feedback line" runs out of the matrix, it causes all feedback points to be toggled.
260 // 0 \ 0 0 \ 1 0
261
262 // The top row of M0^(number_of_bits-1) is (0 0 0 0 0 ... 0 0 1)
263 libecc::bitset<number_of_bits>* row0 = new libecc::bitset<number_of_bits>;
264 row0->reset();
265 row0->set<number_of_bits - 1>();
266
267 // The second row of M0^(number_of_bits-1) is equal to the bottom row of M0.
268 libecc::bitset<number_of_bits>* row1 = new libecc::bitset<number_of_bits>("1");
269 for (int f = 0; f < number_of_feedback_points; ++f)
270 row1->set(feedback_points[f]);
271
272 // Because all feedbacks are > n/2, each feedback point 'line' will reach the left
273 // edge at most once and the maximum number of 1's in any column will be
274 // number_of_feedback_points * (number_of_feedback_points + 1) + 1.
275 //
276 // Also, because every column is equal the previous column shifted one bit down
277 // except for the columns of feedback points because there the bit can be toggled
278 // when a "feedback line" (see graph above) runs left outside the matrix, we only
279 // need to store the first column and the columns that correspond to feedback points.
280 static unsigned short column[number_of_feedback_points + 1][number_of_feedback_points * (number_of_feedback_points + 1) + 2];
281 // Empty everything.
282 for (int i = 0; i <= number_of_feedback_points; ++i)
283 column[i][0] = 0; // No 1's in this column.
284
285 // Generation of the first column, and all columns corresponding to feedback points, of M0^(number_of_bits-1).
286 // Index of 0 corresponds to the first column, larger indices correspond to column 'feedback_points[index + 1]'.
287 for (unsigned short r = 0; r < number_of_bits; ++r) // r is the row of m0^(number_of_bits-1) that row0 refers to.
288 {
289 Dout(dc::notice, "row0 = " << row0);
290 // First column.
291 if (row0->test<0>())
292 column[0][++column[0][0]] = r;
293 // Columns corresponding to feedback points.
294 for (int f = 0; f < number_of_feedback_points; ++f)
295 {
296 if (row0->test(feedback_points[f]))
297 column[f + 1][++column[f + 1][0]] = r;
298 }
299 // Generate next `row'.
300 bool last_bit_set = row0->test<number_of_bits - 1>();
301 row0->shift_op<1, libecc::left, libecc::assign>(*row0);
302 if (last_bit_set)
303 *row0 ^= *row1;
304 }
305
306 // Abuse row0 and row1 alternating as input and output in the loop below.
307 libecc::bitset<number_of_bits>* out = row0;
308 libecc::bitset<number_of_bits>* in = row1;
309
310 // Calculate M0^(2^number_of_bits).
311 // Maxtrices are stored 'compressed' (just the top row).
312 // The top row of M0 is (0 1 0 0 0 ... 0).
313 out->reset();
314 out->set<1>();
315
316 if (number_of_bits > 2000)
317 std::cout << "Calculating... ";
318 // p == log2(current_power).
319 for (int p = 1;;)
320 {
321 Dout(dc::notice, "Calculating M0^(2^" << p << "))");
322
323 // Swap in and out.
324 libecc::bitset<number_of_bits>* tmp = in;
325 in = out;
326 out = tmp;
327
328 if (number_of_bits > 2000)
329 std::cout << 'M' << p << ' ' << std::flush;
330
331 // The square of a compressed matrix C (a vector) is given by
332 // the formula: bit n of C^2 is: CT Qn C, where CT is C transposed
333 // and Qn is a symmetric matrix where column x is the n-th column
334 // M0^x. This implies that the last column of every Qn is n-th
335 // column of M0^(number_of_bits-1) and because of the symmetry
336 // this is equal to the transpose of the bottom row.
337 // It turns out that Qn exists of a mirror-and-rotation part
338 // half of which exists in the left-top part of the matrix and
339 // half of which exists in the right-bottom part. The right
340 // bottom part also contains all bits that represent feedback
341 // points.
342 //
343 // A typical Qn looks as follows
344 // (this is Q1 of the previous example):
345 //
346 // Starts in
347 // col 1 Column 1 of M0^6 (both because this is Q1)
348 // | _______ |__
349 // V/ V \
350 // 0 1 0 0 0 0 0 _\ mirror/rotation part.
351 // 1 0 0 0 0 0 0 /
352 // 0 0 0 0 0 0 1
353 // 0 0 0 0 0 1 0
354 // 0 0 0 0 1 0 0
355 // 0 0 0 1 0 0 0
356 // 0 0 1 0 0 0 / <-- bottom row is transpose of last column.
357 //
358 // (this is Q4):
359 //
360 // Starts in
361 // col 4 Column 4 of M0^6
362 // \ |
363 // V V
364 // 0 0 0 0 1 0 0
365 // 0 0 0 1 0 0 0
366 // 0 0 1 0 0 0 /
367 // 0 1 0 0 0 / 0
368 // 1 0 0 0 / 0 0
369 // 0 0 0 / 0 0 1
370 // 0 0 / 0 0 1 /
371 //
372
373 // Instead of calculating Q_output_bit, we only determine its bottom row.
374 // The bottom row is, as said before, equal to the transpose of the last column.
375 // The last column is equal to column `output_bit' of M0^(number_of_bits-1).
376 // When output_bit == 0, this is equal to the column[0], when `output_bit'
377 // is equal to one of the feedback points, this last column is equal to
378 // column[feedback point + 1]. In all other cases, it is equal to the
379 // column[] that is on its left side, shifted down that distance.
380 // For example:
381 //
382 // Q6 is
383 // Column 6 of M0^6
384 // |
385 // V
386 // 0 0 0 0 0 0 1
387 // 0 0 0 0 0 1 0
388 // 0 0 0 0 1 0 0
389 // 0 0 0 1 0 0 0
390 // 0 0 1 0 0 0 /
391 // 0 1 0 0 0 / 0
392 // 1 0 0 0 / 0 0 <-- bottom row is transpose of column 6 of M0^6.
393 // |
394 // M0^6 is (see|above)
395 // V
396 // 0 0 0 0 0 0 1
397 // 1 0 0 \ 0 0 0
398 // 0 1 0 0 \ 0 0
399 // 0 0 1 0 0 \ 0
400 // 0 0 0 1 0 0 \
401 // \ 0 0 \ 1 0 0
402 // 0 \ 0 0 \ 1 0
403 // ^ ^
404 // | |
405 // | column[1] (first feedback point)
406 // column[0] (first column) Q6 feedback
407 // | |
408 // V V
409 // column[1], transposed and shifted (6 - 3) is: (0 0 0 0 / 0 0) which has
410 // all 1's that correspond to feedback points set as in the bottom row of Q6.
411 //
412 int next_feedback_index = 0;
413 int feedback_column_offset = 0;
414 unsigned short* feedback_column = column[0];
415
416 Dout(dc::notice, "Input is " << in);
417
418 // Calculate the square of 'in' and write the result to 'out'.
419 // We do this bit by bit.
420 for (int output_bit = 0; output_bit < number_of_bits;)
421 {
422 bool result = false;
423 // This bit is the inproduct of `*in' with 'Q_output_bit *in'.
424 // Loop over the bits of the left `*in'.
425 int left_in_bit_start = 0;
426 for (int left_in_bit_digit = 0;
427 left_in_bit_digit < libecc::bitset<number_of_bits>::digits;
428 ++left_in_bit_digit, left_in_bit_start += libecc::bitset_digit_bits)
429 {
430 libecc::bitset_digit_t in_digit = in->digit(left_in_bit_digit);
431 if (in_digit == 0)
432 continue;
433 libecc::bitset_digit_t mask = 1;
434 int left_in_bit = left_in_bit_start;
435 do
436 {
437 // If this bit is not set then it makes no sense to calculate the corresponding bit of Q_output_bit *in.
438 if (in_digit & mask)
439 {
440 // Now we need to toggle `result' if the inproduct
441 // of the right `*in' with the `left_in_bit's row of Q_output_bit is set.
442 // We achieve this by toggling result when the corresponding bit
443 // in `*in' is set for every set bit in `left_in_bit's row of Q_output_bit.
444
445 // First of all, the bits resulting from M0^0 in P (see Theory), the
446 // mirror/rotate part in the left-top part of the Qn matrix.
447 if (left_in_bit <= output_bit && in->test(output_bit - left_in_bit))
448 result = !result; // Toggle output bit.
449
450 // Next, all other bits.
451 for (int i = 1; i <= feedback_column[0]; ++i)
452 {
453 int b = feedback_column[i] + feedback_column_offset + number_of_bits - 1 - left_in_bit;
454 if (b >= number_of_bits)
455 break;
456 if (in->test(b))
457 result = !result; // Toggle output bit.
458 }
459 }
460 mask <<= 1;
461 ++left_in_bit;
462 }
463 while((left_in_bit & (libecc::bitset_digit_bits - 1)));
464 }
465 Dout(dc::notice, "Output bit " << output_bit << " is " << result);
466 if (result)
467 out->set(output_bit);
468 else
469 out->clear(output_bit);
470
471 if (++output_bit == feedback_points[next_feedback_index])
472 {
473 feedback_column = column[next_feedback_index + 1];
474 if (next_feedback_index < number_of_feedback_points - 1)
475 ++next_feedback_index;
476 feedback_column_offset = 0;
477 }
478 else
479 ++feedback_column_offset;
480 }
481
482 if (p++ == number_of_bits)
483 {
484 in->reset();
485 in->set<1>();
486 if (*in == *out)
487 std::cout << "successful, the period is 2^" << number_of_bits << " - 1.";
488 std::cout << std::endl;
489 break;
490 }
491 }
492 if (number_of_bits > 2000)
493 std::cout << '\n';
494 #endif
495
496 //feedback_point_indexes[8] = 54;
497 // Next feedback point
498 for (int f = number_of_feedback_points - 1;; --f)
499 {
500 feedback_points[f] = primes[++feedback_point_indexes[f]];
501 if (feedback_points[f] <= number_of_bits / 2)
502 {
503 while(++f < number_of_feedback_points)
504 {
505 feedback_point_indexes[f] = feedback_point_indexes[f - 1] + 1;
506 feedback_points[f] = primes[feedback_point_indexes[f]];
507 if (feedback_points[f] > number_of_bits / 2)
508 break;
509 }
510 if (feedback_points[f] <= number_of_bits / 2)
511 break;
512 }
513 if (f == 0)
514 return 0;
515 }
516 //feedback_point_indexes[number_of_feedback_points - 1] = 53;
517 //feedback_points[number_of_feedback_points - 1] = primes[53];
518 }
519
520 return 0;
521 }
522
init_primes(void)523 void init_primes(void)
524 {
525 #if 1
526 // Use any feedback point, not just primes.
527 for (int p = 1; p < number_of_bits; ++p)
528 primes.push_back(p);
529 return;
530 #endif
531 primes.push_back(2);
532 int const zs = max_prime / 64 + 1;
533 double mp = (32 * zs - 1) * 2 + 3;
534 unsigned int sr = (unsigned int)(sqrt(mp) + 0.5);
535 uint32_t* z = new uint32_t [zs];
536 std::memset(z, 0, zs * sizeof(uint32_t));
537 for (int i = 0; i < zs; ++i)
538 for (int s = 0; s < 32; ++s)
539 if ((z[i] & (1 << s)) == 0)
540 {
541 unsigned int p = (32 * i + s) * 2 + 3;
542 if (p > max_prime)
543 {
544 delete [] z;
545 return;
546 }
547 primes.push_back(p);
548 if (p >= sr)
549 continue;
550 for (int q = p; q < p * 33; q += p)
551 {
552 uint32_t m = 1 << ((s + q) % 32);
553 for (int j = i + (s + q) / 32; j < zs; j += p)
554 z[j] |= m;
555 }
556 }
557 delete [] z;
558 }
559
init_feedback_points(void)560 void init_feedback_points(void)
561 {
562 #if 1
563 // Use small feedback numbers.
564 for (int f = 0; f < number_of_feedback_points; ++f)
565 {
566 feedback_point_indexes[f] = f /* + 1 */; // +1 = skip feedback of 2.
567 }
568 #else
569 // Manual override.
570 feedback_point_indexes[0] = 0;
571 feedback_point_indexes[1] = 1;
572 feedback_point_indexes[2] = 3;
573 feedback_point_indexes[3] = 5;
574 feedback_point_indexes[4] = 10;
575 feedback_point_indexes[5] = 17;
576 feedback_point_indexes[6] = 31;
577 feedback_point_indexes[7] = 35;
578 feedback_point_indexes[8] = 53;
579 #endif
580
581 for (int f = 0; f < number_of_feedback_points; ++f)
582 {
583 feedback_points[f] = primes[feedback_point_indexes[f]];
584 }
585 }
586
587 struct BitIterator {
588 int cnt;
589 libecc::bitset<number_of_bits> const* bp;
590 libecc::bitset_digit_t bit_mask;
591 int digit;
BitIteratorBitIterator592 BitIterator(Matrix const& m) : cnt(number_of_bits), bp(&m.M_matrix[0]), bit_mask(1), digit(0) { }
next_bitBitIterator593 bool next_bit(void)
594 {
595 if (--cnt == 0)
596 return false;
597 ++bp;
598 bit_mask <<= 1;
599 if (bit_mask == 0)
600 {
601 bit_mask = 1;
602 ++digit;
603 }
604 return true;
605 }
606 };
607
square(Matrix const & m)608 void Matrix::square(Matrix const& m)
609 {
610 for (int row = 0; row < number_of_bits; ++row)
611 M_matrix[row].reset();
612 BitIterator col(m);
613 do
614 {
615 // Make a copy of column col.
616 libecc::bitset<number_of_bits> column;
617 libecc::bitset_digit_t res = 0;
618 BitIterator r(m);
619 do
620 {
621 if (r.bp->digit(col.digit) & col.bit_mask)
622 res |= r.bit_mask;
623 if (r.bit_mask == 0x80000000)
624 {
625 column.rawdigit(r.digit) = res;
626 res = 0;
627 }
628 }
629 while(r.next_bit());
630 column.rawdigit(libecc::bitset<number_of_bits>::digits - 1) = res;
631 // Calculate inproduct of this column with any of the rows
632 for (int row = 0; row < number_of_bits; ++row)
633 {
634 libecc::bitset<number_of_bits> tmp = column & m.M_matrix[row];
635 if (tmp.odd())
636 M_matrix[row].rawdigit(col.digit) |= col.bit_mask;
637 }
638 }
639 while(col.next_bit());
640 }
641
mult(Matrix const & m1,Matrix const & m2)642 void Matrix::mult(Matrix const& m1, Matrix const& m2)
643 {
644 for (int row = 0; row < number_of_bits; ++row)
645 M_matrix[row].reset();
646 BitIterator col(m2);
647 do
648 {
649 // Make a copy of column col.
650 libecc::bitset<number_of_bits> column;
651 libecc::bitset_digit_t res = 0;
652 BitIterator r(m1);
653 do
654 {
655 if (r.bp->digit(col.digit) & col.bit_mask)
656 res |= r.bit_mask;
657 if (r.bit_mask == 0x80000000)
658 {
659 column.rawdigit(r.digit) = res;
660 res = 0;
661 }
662 }
663 while(r.next_bit());
664 column.rawdigit(libecc::bitset<number_of_bits>::digits - 1) = res;
665 // Calculate inproduct of this column with any of the rows
666 for (int row = 0; row < number_of_bits; ++row)
667 {
668 libecc::bitset<number_of_bits> tmp = column & m2.M_matrix[row];
669 if (tmp.odd())
670 M_matrix[row].rawdigit(col.digit) |= col.bit_mask;
671 }
672 }
673 while(col.next_bit());
674 }
675
operator <<(std::ostream & os,Matrix const & m)676 std::ostream& operator<<(std::ostream& os, Matrix const& m)
677 {
678 for (int r = 0; r < number_of_bits; ++r)
679 {
680 for (int c = 0; c < number_of_bits; ++c)
681 {
682 if (m.M_matrix[r].test(c))
683 os << "1 ";
684 else
685 os << "0 ";
686 }
687 os << '\n';
688 }
689 return os;
690 }
691
operator ==(Matrix const & m1,Matrix const & m2)692 bool operator==(Matrix const& m1, Matrix const& m2)
693 {
694 for (int r = 0; r < number_of_bits; ++r)
695 if (m1.M_matrix[r] != m2.M_matrix[r])
696 return false;
697 return true;
698 }
699