IT++ Logo
galois.cpp
Go to the documentation of this file.
1
29#include <itpp/comm/galois.h>
31#include <itpp/base/itcompat.h>
32#include <string>
33#include <iostream>
34
35
36namespace itpp
37{
38
39Array<Array<int> > GF::alphapow;
40Array<Array<int> > GF::logalpha;
41ivec 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
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}
93std::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
125std::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
135std::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
157GFX 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();
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();
172 temp -= divisor * g;
173 tempdegree = temp.get_true_degree();
175 if ((degreedif < 0) || (temp.get_true_degree() == 0 && temp[0] == GF(q, -1))) {
176 break;
177 }
178 }
179 return m;
180}
181
183GFX 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();
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
General array class.
Definition array.h:105
int size() const
Returns the number of data elements in the array object.
Definition array.h:155
void set_size(int n, bool copy=false)
Resizing an Array<T>.
Definition array.h:257
Polynomials over GF(q)[x], where q=2^m, m=1,...,16.
Definition galois.h:176
int get_true_degree() const
Return true degree of GF(q)[x].
Definition galois.h:466
Galois Field GF(q).
Definition galois.h:75
void set_size(int qvalue)
set q=2^mvalue
Definition galois.cpp:44
Definitions of Galois Field algebra classes and functions.
#define it_assert(t, s)
Abort if t is not true.
Definition itassert.h:94
vec log2(const vec &x)
log-2 of the elements
Definition log_exp.cpp:36
int size(const Vec< T > &v)
Length of vector.
Definition matfunc.h:55
IT++ compatibility types and functions.
Logarithmic and exponenential functions - header file.
itpp namespace
Definition itmex.h:37
ITPP_EXPORT int round_i(double x)
Round to nearest integer.
std::ostream & operator<<(std::ostream &output, const bin &inbin)
Output stream of bin.
Definition binary.cpp:36
GFX divgfx(const GFX &c, const GFX &g)
Division of two GFX (local help function)
Definition galois.cpp:157
std::istream & operator>>(std::istream &input, bin &outbin)
Input stream of bin.
Definition binary.cpp:42
GFX modgfx(const GFX &a, const GFX &b)
Modulo function of two GFX (local help function)
Definition galois.cpp:183

Generated on Tue Mar 26 2024 19:08:31 for IT++ by Doxygen 1.9.8