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   Array<Array<int> > GF::alphapow;
00038   Array<Array<int> > GF::logalpha;
00039   ivec GF::q="1 2 4 8 16 32 64 128 256 512 1024 2048 4096 8192 16384 32768 65536";
00040 
00041   // set q=2^mvalue
00042   void GF::set_size(int qvalue)
00043   {
00044     int mtemp;
00045 
00046     mtemp = round_i(::log2(static_cast<double>(qvalue)));
00047     it_assert((1<<mtemp)==qvalue, "GF::setsize : q is not a power of 2");
00048     it_assert(mtemp<=16, "GF::setsize : q must be 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     it_error_if(mtemp < 1 || mtemp > 16, "createfield : m out of range");
00056     m=mtemp;
00057     if (alphapow.size() <  (m+1) ) {
00058       alphapow.set_size(m+1);
00059       logalpha.set_size(m+1);
00060     }
00061 
00062     if (alphapow(m).size() == 0) {
00063       alphapow(m).set_size(qvalue);
00064       logalpha(m).set_size(qvalue);
00065       alphapow(m) = 0;
00066       logalpha(m) = 0;
00067       if (m == 1) { // GF(2), special case
00068         alphapow(1)(0)=1;
00069         logalpha(1)(0)=-1; logalpha(1)(1)=0;
00070       } else {
00071         reduce=reducetable[m-2];
00072         alphapow(m)(0)=1; // alpha^0 = 1
00073         for (n=1; n<(1<<m)-1; n++) {
00074           temp=alphapow(m)(n-1);
00075           temp=(temp << 1); // multiply by alpha
00076           if (temp & (1<<m)) // contains alpha**m term
00077             alphapow(m)(n)=(temp & ~(1<<m))^reduce;
00078           else
00079             alphapow(m)(n)=temp; // if no alpha**m term, store as is
00080 
00081           // create table to go in opposite direction
00082           logalpha(m)(0)=-1; // special case, actually log(0)=-inf
00083         }
00084 
00085         for (n=0;n<(1<<m)-1;n++)
00086           logalpha(m)(alphapow(m)(n))=n;
00087       }
00088     }
00089   }
00090 
00092   std::ostream &operator<<(std::ostream &os, const GF &ingf)
00093   {
00094     if (ingf.value == -1)
00095       os << "0";
00096     else
00097       os << "alpha^" << ingf.value;
00098     return os;
00099   }
00100 
00102   std::ostream &operator<<(std::ostream &os, const GFX &ingfx)
00103   {
00104     int terms=0;
00105     for (int i=0; i<ingfx.degree+1; i++) {
00106       if (ingfx.coeffs(i) != GF(ingfx.q,-1) ) {
00107         if (terms != 0) os << " + ";
00108         terms++;
00109         if (ingfx.coeffs(i) == GF(ingfx.q,0) ) {// is the coefficient an one (=alpha^0=1)
00110           os  << "x^" << i;
00111         } else {
00112           os  << ingfx.coeffs(i) << "*x^" << i;
00113         }
00114       }
00115     }
00116     if (terms == 0) os << "0";
00117     return os;
00118   }
00119 
00120   //----------------- Help Functions -----------------
00121 
00123   GFX divgfx(const GFX &c, const GFX &g) {
00124     int q = c.get_size();
00125     GFX temp = c;
00126     int tempdegree = temp.get_true_degree();
00127     int gdegree = g.get_true_degree();
00128     int degreedif = tempdegree - gdegree;
00129     if (degreedif < 0) return GFX(q,0); // denominator larger than nominator. Return zero polynomial.
00130     GFX m(q,degreedif), divisor(q);
00131 
00132     for (int i=0; i<c.get_degree(); i++) {
00133       m[degreedif] = temp[tempdegree]/g[gdegree];
00134       divisor.set_degree(degreedif);
00135       divisor.clear();
00136       divisor[degreedif] = m[degreedif];
00137       temp -= divisor*g;
00138       tempdegree = temp.get_true_degree();
00139       degreedif = tempdegree - gdegree;
00140       if ( (degreedif<0) || (temp.get_true_degree()==0 && temp[0] == GF(q,-1) ) ) {
00141         break;
00142       }
00143     }
00144     return m;
00145   }
00146 
00148   GFX modgfx(const GFX &a, const GFX &b)
00149   {
00150     int q = a.get_size();
00151     GFX temp = a;
00152     int tempdegree = temp.get_true_degree();
00153     int bdegree = b.get_true_degree();
00154     int degreedif = a.get_true_degree() - b.get_true_degree();
00155     if (degreedif < 0) return temp; // Denominator larger than nominator. Return nominator.
00156     GFX m(q,degreedif), divisor(q);
00157 
00158     for (int i=0; i<a.get_degree(); i++) {
00159       m[degreedif] = temp[tempdegree]/b[bdegree];
00160       divisor.set_degree(degreedif);
00161       divisor.clear();
00162       divisor[degreedif] =  m[degreedif];
00163       temp -= divisor*b; // Bug-fixed. Used to be: temp -= divisor*a;
00164       tempdegree = temp.get_true_degree();
00165       degreedif = temp.get_true_degree() - bdegree;
00166       if ( (degreedif<0) || (temp.get_true_degree()==0 && temp[0] == GF(q,-1) ) ) {
00167         break;
00168       }
00169     }
00170     return temp;
00171   }
00172 
00173 } // namespace itpp
SourceForge Logo

Generated on Sun Dec 9 17:30:26 2007 for IT++ by Doxygen 1.5.4