IT++ Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
galois.cpp
Go to the documentation of this file.
1 
29 #include <itpp/comm/galois.h>
30 #include <itpp/base/math/log_exp.h>
31 #include <itpp/base/itcompat.h>
32 #include <iostream>
33 
34 
35 namespace itpp
36 {
37 
38 Array<Array<int> > GF::alphapow;
39 Array<Array<int> > GF::logalpha;
40 ivec GF::q = "1 2 4 8 16 32 64 128 256 512 1024 2048 4096 8192 16384 32768 65536";
41 
42 // set q=2^mvalue
43 void GF::set_size(int qvalue)
44 {
45  m = static_cast<char>(round_i(::log2(static_cast<double>(qvalue))));
46  it_assert((1 << m) == qvalue, "GF::setsize : q is not a power of 2");
47  it_assert((m > 0) && (m <= 16), "GF::setsize : q must be positive and "
48  "less than or equal to 2^16");
49 
50  /* Construct GF(q), q=2^m. From Wicker, "Error Control Systems
51  for digital communication and storage" pp. 463-465 */
52 
53  int reduce, temp, n;
54  const int reducetable[] = {3, 3, 3, 5, 3, 9, 29, 17, 9, 5, 83, 27, 43, 3, 4107}; // starts at m=2,..,16
55 
56  if (alphapow.size() < (m + 1)) {
57  alphapow.set_size(m + 1, true);
58  logalpha.set_size(m + 1, true);
59  }
60 
61  if (alphapow(m).size() == 0) {
62  alphapow(m).set_size(qvalue);
63  logalpha(m).set_size(qvalue);
64  alphapow(m) = 0;
65  logalpha(m) = 0;
66  if (m == 1) { // GF(2), special case
67  alphapow(1)(0) = 1;
68  logalpha(1)(0) = -1;
69  logalpha(1)(1) = 0;
70  }
71  else {
72  reduce = reducetable[m-2];
73  alphapow(m)(0) = 1; // alpha^0 = 1
74  for (n = 1; n < (1 << m) - 1; n++) {
75  temp = alphapow(m)(n - 1);
76  temp = (temp << 1); // multiply by alpha
77  if (temp & (1 << m)) // contains alpha**m term
78  alphapow(m)(n) = (temp & ~(1 << m)) ^ reduce;
79  else
80  alphapow(m)(n) = temp; // if no alpha**m term, store as is
81 
82  // create table to go in opposite direction
83  logalpha(m)(0) = -1; // special case, actually log(0)=-inf
84  }
85 
86  for (n = 0;n < (1 << m) - 1;n++)
87  logalpha(m)(alphapow(m)(n)) = n;
88  }
89  }
90 }
91 
93 std::ostream &operator<<(std::ostream &os, const GF &ingf)
94 {
95  if (ingf.value == -1)
96  os << "0";
97  else
98  os << "alpha^" << ingf.value;
99  return os;
100 }
101 
103 std::ostream &operator<<(std::ostream &os, const GFX &ingfx)
104 {
105  int terms = 0;
106  for (int i = 0; i < ingfx.degree + 1; i++) {
107  if (ingfx.coeffs(i) != GF(ingfx.q, -1)) {
108  if (terms != 0) os << " + ";
109  terms++;
110  if (ingfx.coeffs(i) == GF(ingfx.q, 0)) {// is the coefficient an one (=alpha^0=1)
111  os << "x^" << i;
112  }
113  else {
114  os << ingfx.coeffs(i) << "*x^" << i;
115  }
116  }
117  }
118  if (terms == 0) os << "0";
119  return os;
120 }
121 
122 //----------------- Help Functions -----------------
123 
125 GFX divgfx(const GFX &c, const GFX &g)
126 {
127  int q = c.get_size();
128  GFX temp = c;
129  int tempdegree = temp.get_true_degree();
130  int gdegree = g.get_true_degree();
131  int degreedif = tempdegree - gdegree;
132  if (degreedif < 0) return GFX(q, 0); // denominator larger than nominator. Return zero polynomial.
133  GFX m(q, degreedif), divisor(q);
134 
135  for (int i = 0; i < c.get_degree(); i++) {
136  m[degreedif] = temp[tempdegree] / g[gdegree];
137  divisor.set_degree(degreedif);
138  divisor.clear();
139  divisor[degreedif] = m[degreedif];
140  temp -= divisor * g;
141  tempdegree = temp.get_true_degree();
142  degreedif = tempdegree - gdegree;
143  if ((degreedif < 0) || (temp.get_true_degree() == 0 && temp[0] == GF(q, -1))) {
144  break;
145  }
146  }
147  return m;
148 }
149 
151 GFX modgfx(const GFX &a, const GFX &b)
152 {
153  int q = a.get_size();
154  GFX temp = a;
155  int tempdegree = temp.get_true_degree();
156  int bdegree = b.get_true_degree();
157  int degreedif = a.get_true_degree() - b.get_true_degree();
158  if (degreedif < 0) return temp; // Denominator larger than nominator. Return nominator.
159  GFX m(q, degreedif), divisor(q);
160 
161  for (int i = 0; i < a.get_degree(); i++) {
162  m[degreedif] = temp[tempdegree] / b[bdegree];
163  divisor.set_degree(degreedif);
164  divisor.clear();
165  divisor[degreedif] = m[degreedif];
166  temp -= divisor * b; // Bug-fixed. Used to be: temp -= divisor*a;
167  tempdegree = temp.get_true_degree();
168  degreedif = temp.get_true_degree() - bdegree;
169  if ((degreedif < 0) || (temp.get_true_degree() == 0 && temp[0] == GF(q, -1))) {
170  break;
171  }
172  }
173  return temp;
174 }
175 
176 } // namespace itpp
SourceForge Logo

Generated on Fri Mar 21 2014 17:14:13 for IT++ by Doxygen 1.8.1.2