IT++ Logo Newcom Logo

galois.cpp

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

Generated on Thu Apr 19 14:18:30 2007 for IT++ by Doxygen 1.5.1