IT++ Logo

galois.cpp

Go to the documentation of this file.
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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
SourceForge Logo

Generated on Sat Feb 26 2011 16:06:32 for IT++ by Doxygen 1.7.3