00001 00030 #include <itpp/comm/galois.h> 00031 #include <itpp/base/math/log_exp.h> 00032 #include <iostream> 00033 00034 00035 namespace itpp 00036 { 00037 00038 Array<Array<int> > GF::alphapow; 00039 Array<Array<int> > GF::logalpha; 00040 ivec GF::q = "1 2 4 8 16 32 64 128 256 512 1024 2048 4096 8192 16384 32768 65536"; 00041 00042 // set q=2^mvalue 00043 void GF::set_size(int qvalue) 00044 { 00045 m = static_cast<char>(round_i(::log2(static_cast<double>(qvalue)))); 00046 it_assert((1 << m) == qvalue, "GF::setsize : q is not a power of 2"); 00047 it_assert((m > 0) && (m <= 16), "GF::setsize : q must be positive and " 00048 "less than or equal to 2^16"); 00049 00050 /* Construct GF(q), q=2^m. From Wicker, "Error Control Systems 00051 for digital communication and storage" pp. 463-465 */ 00052 00053 int reduce, temp, n; 00054 const int reducetable[] = {3, 3, 3, 5, 3, 9, 29, 17, 9, 5, 83, 27, 43, 3, 4107}; // starts at m=2,..,16 00055 00056 if (alphapow.size() < (m + 1)) { 00057 alphapow.set_size(m + 1, true); 00058 logalpha.set_size(m + 1, true); 00059 } 00060 00061 if (alphapow(m).size() == 0) { 00062 alphapow(m).set_size(qvalue); 00063 logalpha(m).set_size(qvalue); 00064 alphapow(m) = 0; 00065 logalpha(m) = 0; 00066 if (m == 1) { // GF(2), special case 00067 alphapow(1)(0) = 1; 00068 logalpha(1)(0) = -1; 00069 logalpha(1)(1) = 0; 00070 } 00071 else { 00072 reduce = reducetable[m-2]; 00073 alphapow(m)(0) = 1; // alpha^0 = 1 00074 for (n = 1; n < (1 << m) - 1; n++) { 00075 temp = alphapow(m)(n - 1); 00076 temp = (temp << 1); // multiply by alpha 00077 if (temp & (1 << m)) // contains alpha**m term 00078 alphapow(m)(n) = (temp & ~(1 << m)) ^ reduce; 00079 else 00080 alphapow(m)(n) = temp; // if no alpha**m term, store as is 00081 00082 // create table to go in opposite direction 00083 logalpha(m)(0) = -1; // special case, actually log(0)=-inf 00084 } 00085 00086 for (n = 0;n < (1 << m) - 1;n++) 00087 logalpha(m)(alphapow(m)(n)) = n; 00088 } 00089 } 00090 } 00091 00093 std::ostream &operator<<(std::ostream &os, const GF &ingf) 00094 { 00095 if (ingf.value == -1) 00096 os << "0"; 00097 else 00098 os << "alpha^" << ingf.value; 00099 return os; 00100 } 00101 00103 std::ostream &operator<<(std::ostream &os, const GFX &ingfx) 00104 { 00105 int terms = 0; 00106 for (int i = 0; i < ingfx.degree + 1; i++) { 00107 if (ingfx.coeffs(i) != GF(ingfx.q, -1)) { 00108 if (terms != 0) os << " + "; 00109 terms++; 00110 if (ingfx.coeffs(i) == GF(ingfx.q, 0)) {// is the coefficient an one (=alpha^0=1) 00111 os << "x^" << i; 00112 } 00113 else { 00114 os << ingfx.coeffs(i) << "*x^" << i; 00115 } 00116 } 00117 } 00118 if (terms == 0) os << "0"; 00119 return os; 00120 } 00121 00122 //----------------- Help Functions ----------------- 00123 00125 GFX divgfx(const GFX &c, const GFX &g) 00126 { 00127 int q = c.get_size(); 00128 GFX temp = c; 00129 int tempdegree = temp.get_true_degree(); 00130 int gdegree = g.get_true_degree(); 00131 int degreedif = tempdegree - gdegree; 00132 if (degreedif < 0) return GFX(q, 0); // denominator larger than nominator. Return zero polynomial. 00133 GFX m(q, degreedif), divisor(q); 00134 00135 for (int i = 0; i < c.get_degree(); i++) { 00136 m[degreedif] = temp[tempdegree] / g[gdegree]; 00137 divisor.set_degree(degreedif); 00138 divisor.clear(); 00139 divisor[degreedif] = m[degreedif]; 00140 temp -= divisor * g; 00141 tempdegree = temp.get_true_degree(); 00142 degreedif = tempdegree - gdegree; 00143 if ((degreedif < 0) || (temp.get_true_degree() == 0 && temp[0] == GF(q, -1))) { 00144 break; 00145 } 00146 } 00147 return m; 00148 } 00149 00151 GFX modgfx(const GFX &a, const GFX &b) 00152 { 00153 int q = a.get_size(); 00154 GFX temp = a; 00155 int tempdegree = temp.get_true_degree(); 00156 int bdegree = b.get_true_degree(); 00157 int degreedif = a.get_true_degree() - b.get_true_degree(); 00158 if (degreedif < 0) return temp; // Denominator larger than nominator. Return nominator. 00159 GFX m(q, degreedif), divisor(q); 00160 00161 for (int i = 0; i < a.get_degree(); i++) { 00162 m[degreedif] = temp[tempdegree] / b[bdegree]; 00163 divisor.set_degree(degreedif); 00164 divisor.clear(); 00165 divisor[degreedif] = m[degreedif]; 00166 temp -= divisor * b; // Bug-fixed. Used to be: temp -= divisor*a; 00167 tempdegree = temp.get_true_degree(); 00168 degreedif = temp.get_true_degree() - bdegree; 00169 if ((degreedif < 0) || (temp.get_true_degree() == 0 && temp[0] == GF(q, -1))) { 00170 break; 00171 } 00172 } 00173 return temp; 00174 } 00175 00176 } // namespace itpp
Generated on Sun Jul 26 08:36:49 2009 for IT++ by Doxygen 1.5.9