IT++ Logo
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 <string>
33 #include <iostream>
34 
35 
36 namespace itpp
37 {
38 
39 Array<Array<int> > GF::alphapow;
40 Array<Array<int> > GF::logalpha;
41 ivec GF::q = "1 2 4 8 16 32 64 128 256 512 1024 2048 4096 8192 16384 32768 65536";
42 
43 // set q=2^mvalue
44 void GF::set_size(int qvalue)
45 {
46  m = static_cast<char>(round_i(::log2(static_cast<double>(qvalue))));
47  it_assert((1 << m) == qvalue, "GF::setsize : q is not a power of 2");
48  it_assert((m > 0) && (m <= 16), "GF::setsize : q must be positive and "
49  "less than or equal to 2^16");
50 
51  /* Construct GF(q), q=2^m. From Wicker, "Error Control Systems
52  for digital communication and storage" pp. 463-465 */
53 
54  int reduce, temp, n;
55  const int reducetable[] = {3, 3, 3, 5, 3, 9, 29, 17, 9, 5, 83, 27, 43, 3, 4107}; // starts at m=2,..,16
56 
57  if (alphapow.size() < (m + 1)) {
58  alphapow.set_size(m + 1, true);
59  logalpha.set_size(m + 1, true);
60  }
61 
62  if (alphapow(m).size() == 0) {
63  alphapow(m).set_size(qvalue);
64  logalpha(m).set_size(qvalue);
65  alphapow(m) = 0;
66  logalpha(m) = 0;
67  if (m == 1) { // GF(2), special case
68  alphapow(1)(0) = 1;
69  logalpha(1)(0) = -1;
70  logalpha(1)(1) = 0;
71  }
72  else {
73  reduce = reducetable[m-2];
74  alphapow(m)(0) = 1; // alpha^0 = 1
75  for (n = 1; n < (1 << m) - 1; n++) {
76  temp = alphapow(m)(n - 1);
77  temp = (temp << 1); // multiply by alpha
78  if (temp & (1 << m)) // contains alpha**m term
79  alphapow(m)(n) = (temp & ~(1 << m)) ^ reduce;
80  else
81  alphapow(m)(n) = temp; // if no alpha**m term, store as is
82 
83  // create table to go in opposite direction
84  logalpha(m)(0) = -1; // special case, actually log(0)=-inf
85  }
86 
87  for (n = 0;n < (1 << m) - 1;n++)
88  logalpha(m)(alphapow(m)(n)) = n;
89  }
90  }
91 }
93 std::istream &operator>>(std::istream &is, GF &ingf)
94 {
95  int val; char c;
96  static const std::string prefix("alpha^");
97  c = is.get();
98  if(c == 'a') {
99  //read alpha^pow form from stream
100  std::string::const_iterator pr_it = prefix.begin(); pr_it++;
101  for(; pr_it < prefix.end(); ++pr_it) {
102  c = is.get();
103  if(*pr_it != c) {
104  is.setstate(std::ios_base::failbit);
105  return is;
106  }
107  }
108  is >> val;
109  if(is) ingf.set(ingf.get_size(),val);
110  }
111  else {
112  //try to read 0 from stream
113  is >> val;
114  if(is && (val==0)) {
115  ingf.set(ingf.get_size(),0);
116  }
117  else {
118  is.setstate(std::ios_base::failbit);
119  }
120  }
121  return is;
122 }
123 
125 std::ostream &operator<<(std::ostream &os, const GF &ingf)
126 {
127  if (ingf.value == -1)
128  os << "0";
129  else
130  os << "alpha^" << ingf.value;
131  return os;
132 }
133 
135 std::ostream &operator<<(std::ostream &os, const GFX &ingfx)
136 {
137  int terms = 0;
138  for (int i = 0; i < ingfx.degree + 1; i++) {
139  if (ingfx.coeffs(i) != GF(ingfx.q, -1)) {
140  if (terms != 0) os << " + ";
141  terms++;
142  if (ingfx.coeffs(i) == GF(ingfx.q, 0)) {// is the coefficient an one (=alpha^0=1)
143  os << "x^" << i;
144  }
145  else {
146  os << ingfx.coeffs(i) << "*x^" << i;
147  }
148  }
149  }
150  if (terms == 0) os << "0";
151  return os;
152 }
153 
154 //----------------- Help Functions -----------------
155 
157 GFX divgfx(const GFX &c, const GFX &g)
158 {
159  int q = c.get_size();
160  GFX temp = c;
161  int tempdegree = temp.get_true_degree();
162  int gdegree = g.get_true_degree();
163  int degreedif = tempdegree - gdegree;
164  if (degreedif < 0) return GFX(q, 0); // denominator larger than nominator. Return zero polynomial.
165  GFX m(q, degreedif), divisor(q);
166 
167  for (int i = 0; i < c.get_degree(); i++) {
168  m[degreedif] = temp[tempdegree] / g[gdegree];
169  divisor.set_degree(degreedif);
170  divisor.clear();
171  divisor[degreedif] = m[degreedif];
172  temp -= divisor * g;
173  tempdegree = temp.get_true_degree();
174  degreedif = tempdegree - gdegree;
175  if ((degreedif < 0) || (temp.get_true_degree() == 0 && temp[0] == GF(q, -1))) {
176  break;
177  }
178  }
179  return m;
180 }
181 
183 GFX modgfx(const GFX &a, const GFX &b)
184 {
185  int q = a.get_size();
186  GFX temp = a;
187  int tempdegree = temp.get_true_degree();
188  int bdegree = b.get_true_degree();
189  int degreedif = a.get_true_degree() - b.get_true_degree();
190  if (degreedif < 0) return temp; // Denominator larger than nominator. Return nominator.
191  GFX m(q, degreedif), divisor(q);
192 
193  for (int i = 0; i < a.get_degree(); i++) {
194  m[degreedif] = temp[tempdegree] / b[bdegree];
195  divisor.set_degree(degreedif);
196  divisor.clear();
197  divisor[degreedif] = m[degreedif];
198  temp -= divisor * b; // Bug-fixed. Used to be: temp -= divisor*a;
199  tempdegree = temp.get_true_degree();
200  degreedif = temp.get_true_degree() - bdegree;
201  if ((degreedif < 0) || (temp.get_true_degree() == 0 && temp[0] == GF(q, -1))) {
202  break;
203  }
204  }
205  return temp;
206 }
207 
208 } // namespace itpp
SourceForge Logo

Generated on Sat Jul 6 2013 10:54:22 for IT++ by Doxygen 1.8.2