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
Generated on Thu Apr 19 14:14:58 2007 for IT++ by Doxygen 1.5.1