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
Generated on Sat Apr 19 10:59:23 2008 for IT++ by Doxygen 1.5.5