[288382] | 1 | /////////////////////////////////////////////////////////////////////////////// |
---|
[7f3ad4] | 2 | version="$Id: decodegb.lib,v 1.9 2009-01-14 16:07:03 Singular Exp $"; |
---|
[288382] | 3 | category="Coding theory"; |
---|
| 4 | info=" |
---|
[fba86f8] | 5 | LIBRARY: decodegb.lib Decoding and min distance of linear codes with GB |
---|
[7f3ad4] | 6 | AUTHOR: Stanislav Bulygin, bulygin@mathematik.uni-kl.de |
---|
[288382] | 7 | |
---|
| 8 | OVERVIEW: |
---|
[fba86f8] | 9 | In this library we generate several systems used for decoding cyclic codes and |
---|
[7f3ad4] | 10 | finding their minimum distance. Namely, we work with the Cooper's philosophy |
---|
| 11 | and generalized Newton identities. The original method of quadratic equations |
---|
| 12 | is worked out here as well. We also (for comparison) enable to work with the |
---|
| 13 | system of Fitzgerald-Lax. We provide some auxiliary functions for further |
---|
| 14 | manipulations and decoding. For an overview of the methods mentioned above @ref{Decoding codes with GB} |
---|
| 15 | For the vanishing ideal computation the algorithm of Farr and Gao is |
---|
[fba86f8] | 16 | implemented. |
---|
[288382] | 17 | |
---|
[7f3ad4] | 18 | MAIN PROCEDURES: |
---|
[fba86f8] | 19 | sysCRHT(..); generates the CRHT-ideal as in Cooper's philosophy |
---|
| 20 | sysCRHTMindist(..); CRHT-ideal to find the minimum distance in the binary case |
---|
| 21 | sysNewton(..); generates the ideal with the generalized Newton identities |
---|
| 22 | sysBin(..); generates Bin system using Waring function |
---|
| 23 | encode(x,g); encodes given message x with the given generator matrix g |
---|
| 24 | syndrome(h,c); computes a syndrome w.r.t. the given check matrix |
---|
| 25 | sysQE(..); generates the system of quadratic equations for decoding |
---|
| 26 | errorInsert(..); inserts errors in a word |
---|
| 27 | errorRand(y,num,e); inserts random errors in a word |
---|
| 28 | randomCheck(m,n,e); generates a random check matrix |
---|
| 29 | genMDSMat(n,a); generates an MDS (actually an RS) matrix |
---|
| 30 | mindist(check); computes the minimum distance of a code |
---|
[7f3ad4] | 31 | decode(rec); decoding of a word using the system of quadratic equations |
---|
[fba86f8] | 32 | decodeRandom(..); a procedure for manipulation with random codes |
---|
| 33 | decodeCode(..); a procedure for manipulation with the given code |
---|
| 34 | vanishId(points); computes the vanishing ideal for the given set of points |
---|
| 35 | sysFL(..); generates the Fitzgerald-Lax system |
---|
| 36 | decodeRandomFL(..); manipulation with random codes via Fitzgerald-Lax |
---|
[288382] | 37 | |
---|
[7f3ad4] | 38 | |
---|
| 39 | KEYWORDS: Cyclic code; Linear code; Decoding; |
---|
[288382] | 40 | Minimum distance; Groebner bases |
---|
| 41 | "; |
---|
| 42 | |
---|
| 43 | LIB "linalg.lib"; |
---|
| 44 | LIB "brnoeth.lib"; |
---|
| 45 | |
---|
| 46 | /////////////////////////////////////////////////////////////////////////////// |
---|
[6e118cf] | 47 | // creates a list result, where result[i]=i, 1<=i<=n |
---|
[288382] | 48 | static proc lis (int n) |
---|
| 49 | { |
---|
| 50 | list result; |
---|
| 51 | if (n<=0) {print("ERRORlis");} |
---|
| 52 | for (int i=1; i<=n; i++) |
---|
| 53 | { |
---|
| 54 | result=result+list(i); |
---|
| 55 | } |
---|
| 56 | return(result); |
---|
| 57 | } |
---|
| 58 | |
---|
[6e118cf] | 59 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 60 | // creates a list of all combinations without repititions of m objects out of n |
---|
[288382] | 61 | static proc combinations (int m, int n) |
---|
| 62 | { |
---|
| 63 | list result; |
---|
| 64 | if (m>n) {print("ERRORcombinations");} |
---|
| 65 | if (m==n) {result[size(result)+1]=lis(m);return(result);} |
---|
| 66 | if (m==0) {result[size(result)+1]=list();return(result);} |
---|
| 67 | list temp=combinations(m-1,n-1); |
---|
| 68 | for (int i=1; i<=size(temp); i++) |
---|
| 69 | { |
---|
| 70 | temp[i]=temp[i]+list(n); |
---|
| 71 | } |
---|
| 72 | result=combinations(m,n-1)+temp; |
---|
| 73 | return(result); |
---|
| 74 | } |
---|
| 75 | |
---|
| 76 | |
---|
[6e118cf] | 77 | /////////////////////////////////////////////////////////////////////////////// |
---|
[7f3ad4] | 78 | // the polynomial for Sala's restrictions |
---|
[288382] | 79 | static proc p_poly(int n, int a, int b) |
---|
| 80 | { |
---|
| 81 | poly f; |
---|
| 82 | for (int i=0; i<=n-1; i++) |
---|
| 83 | { |
---|
| 84 | f=f+Z(a)^i*Z(b)^(n-1-i); |
---|
| 85 | } |
---|
| 86 | return(f); |
---|
| 87 | } |
---|
| 88 | |
---|
[6e118cf] | 89 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 90 | |
---|
[fba86f8] | 91 | proc sysCRHT (int n, list defset, int e, int q, int m, list #) |
---|
| 92 | "USAGE: sysCRHT(n,defset,e,q,m,[k]); n,e,q,m,k int, defset list of int's |
---|
| 93 | @format |
---|
[7f3ad4] | 94 | - n length of the cyclic code, |
---|
[fba86f8] | 95 | - defset is a list representing the defining set, |
---|
[7f3ad4] | 96 | - e the error-correcting capacity, |
---|
[fba86f8] | 97 | - q field size |
---|
[7f3ad4] | 98 | - m degree extension of the splitting field, |
---|
| 99 | - if k>0 additional equations representing the fact that every two |
---|
[fba86f8] | 100 | error positions are either different or at least one of them is zero |
---|
| 101 | @end format |
---|
[7f3ad4] | 102 | RETURN: the ring to work with the CRHT-ideal (with Sala's additions), |
---|
[fba86f8] | 103 | containig an ideal with name 'crht' |
---|
[7f3ad4] | 104 | THEORY: Based on 'defset' of the given cyclic code, the procedure constructs |
---|
| 105 | the corresponding Cooper-Reed-Heleseth-Truong ideal 'crht'. With its |
---|
[24f6cd9] | 106 | help one can solve the decoding problem. For basics of the method @ref{Cooper philosophy}. |
---|
[3599e9] | 107 | SEE ALSO: sysNewton, sysBin |
---|
[fba86f8] | 108 | EXAMPLE: example sysCRHT; shows an example |
---|
[288382] | 109 | " |
---|
| 110 | { |
---|
[7f3ad4] | 111 | int r=size(defset); |
---|
[288382] | 112 | ring @crht=(q,a),(Y(e..1),Z(1..e),X(r..1)),lp; |
---|
| 113 | ideal crht; |
---|
| 114 | int i,j; |
---|
| 115 | poly sum; |
---|
[fba86f8] | 116 | int k; |
---|
[7f3ad4] | 117 | if ( size(#) > 0) |
---|
| 118 | { |
---|
| 119 | k = #[1]; |
---|
[fba86f8] | 120 | } |
---|
[7f3ad4] | 121 | |
---|
[fba86f8] | 122 | //---------------------- add check equations -------------------------- |
---|
[288382] | 123 | for (i=1; i<=r; i++) |
---|
| 124 | { |
---|
| 125 | sum=0; |
---|
| 126 | for (j=1; j<=e; j++) |
---|
| 127 | { |
---|
| 128 | sum=sum+Y(j)*Z(j)^defset[i]; |
---|
| 129 | } |
---|
| 130 | crht[i]=sum-X(i); |
---|
| 131 | } |
---|
[7f3ad4] | 132 | |
---|
[fba86f8] | 133 | //--------------------- field equations on syndromes ------------------ |
---|
[288382] | 134 | for (i=1; i<=r; i++) |
---|
| 135 | { |
---|
| 136 | crht=crht,X(i)^(q^m)-X(i); |
---|
| 137 | } |
---|
[7f3ad4] | 138 | |
---|
[fba86f8] | 139 | //------ restrictions on error-locations: n-th roots of unity ---------- |
---|
[288382] | 140 | for (i=1; i<=e; i++) |
---|
| 141 | { |
---|
| 142 | crht=crht,Z(i)^(n+1)-Z(i); |
---|
| 143 | } |
---|
[7f3ad4] | 144 | |
---|
[288382] | 145 | for (i=1; i<=e; i++) |
---|
| 146 | { |
---|
| 147 | crht=crht,Y(i)^(q-1)-1; |
---|
[7f3ad4] | 148 | } |
---|
| 149 | |
---|
[fba86f8] | 150 | //--------- add Sala's additional conditions if necessary -------------- |
---|
| 151 | if ( k > 0 ) |
---|
[7f3ad4] | 152 | |
---|
[288382] | 153 | { |
---|
| 154 | for (i=1; i<=e; i++) |
---|
| 155 | { |
---|
| 156 | for (j=i+1; j<=e; j++) |
---|
| 157 | { |
---|
| 158 | crht=crht,Z(i)*Z(j)*p_poly(n,i,j); |
---|
| 159 | } |
---|
| 160 | } |
---|
[7f3ad4] | 161 | } |
---|
| 162 | export crht; |
---|
| 163 | return(@crht); |
---|
| 164 | } |
---|
[fba86f8] | 165 | example |
---|
| 166 | { "EXAMPLE:"; echo=2; |
---|
[288382] | 167 | // binary cyclic [15,7,5] code with defining set (1,3) |
---|
[fba86f8] | 168 | intvec v = option(get); |
---|
| 169 | |
---|
[7f3ad4] | 170 | list defset=1,3; // defining set |
---|
[fba86f8] | 171 | int n=15; // length |
---|
| 172 | int e=2; // error-correcting capacity |
---|
| 173 | int q=2; // basefield size |
---|
| 174 | int m=4; // degree extension of the splitting field |
---|
| 175 | int sala=1; // indicator to add additional equations |
---|
[7f3ad4] | 176 | |
---|
| 177 | def A=sysCRHT(n,defset,e,q,m); |
---|
[288382] | 178 | setring A; |
---|
[fba86f8] | 179 | A; // shows the ring we are working in |
---|
[7f3ad4] | 180 | print(crht); // the CRHT-ideal |
---|
[288382] | 181 | option(redSB); |
---|
[fba86f8] | 182 | ideal red_crht=std(crht); // reduced Groebner basis |
---|
[288382] | 183 | print(red_crht); |
---|
[7f3ad4] | 184 | |
---|
[288382] | 185 | //============================ |
---|
[7f3ad4] | 186 | A=sysCRHT(n,defset,e,q,m,sala); |
---|
| 187 | setring A; |
---|
[fba86f8] | 188 | print(crht); // CRHT-ideal with additional equations from Sala |
---|
[288382] | 189 | option(redSB); |
---|
[fba86f8] | 190 | ideal red_crht=std(crht); // reduced Groebner basis |
---|
[7f3ad4] | 191 | print(red_crht); |
---|
[fba86f8] | 192 | red_crht[5]; // general error-locator polynomial for this code |
---|
| 193 | option(set,v); |
---|
[288382] | 194 | } |
---|
| 195 | |
---|
[6e118cf] | 196 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 197 | |
---|
[288382] | 198 | |
---|
[fba86f8] | 199 | proc sysCRHTMindist (int n, list defset, int w) |
---|
[7f3ad4] | 200 | "USAGE: sysCRHTMindist(n,defset,w); n,w are int, defset is list of int's |
---|
[fba86f8] | 201 | @format |
---|
[7f3ad4] | 202 | - n length of the cyclic code, |
---|
[fba86f8] | 203 | - defset is a list representing the defining set, |
---|
| 204 | - w is a candidate for the minimum distance |
---|
| 205 | @end format |
---|
[7f3ad4] | 206 | RETURN: the ring to work with the Sala's ideal for the minimum distance |
---|
[fba86f8] | 207 | containing the ideal with name 'crht_md' |
---|
[7f3ad4] | 208 | THEORY: Based on 'defset' of the given cyclic code, the procedure constructs |
---|
| 209 | the corresponding Cooper-Reed-Heleseth-Truong ideal 'crht_md'. With |
---|
| 210 | its help one can find minimum distance of the code in the binary |
---|
[24f6cd9] | 211 | case. For basics of the method @ref{Cooper philosophy}. |
---|
[fba86f8] | 212 | EXAMPLE: example sysCRHTMindist; shows an example |
---|
[288382] | 213 | " |
---|
| 214 | { |
---|
[7f3ad4] | 215 | int r=size(defset); |
---|
[288382] | 216 | ring @crht_md=2,Z(1..w),lp; |
---|
| 217 | ideal crht_md; |
---|
| 218 | int i,j; |
---|
| 219 | poly sum; |
---|
[7f3ad4] | 220 | |
---|
[6e118cf] | 221 | //------------ add check equations -------------- |
---|
[288382] | 222 | for (i=1; i<=r; i++) |
---|
| 223 | { |
---|
| 224 | sum=0; |
---|
| 225 | for (j=1; j<=w; j++) |
---|
| 226 | { |
---|
| 227 | sum=sum+Z(j)^defset[i]; |
---|
| 228 | } |
---|
| 229 | crht_md[i]=sum; |
---|
[7f3ad4] | 230 | } |
---|
| 231 | |
---|
| 232 | |
---|
[6e118cf] | 233 | //----------- locations are n-th roots of unity ------------ |
---|
[288382] | 234 | for (i=1; i<=w; i++) |
---|
| 235 | { |
---|
| 236 | crht_md=crht_md,Z(i)^n-1; |
---|
[7f3ad4] | 237 | } |
---|
| 238 | |
---|
[6e118cf] | 239 | //------------ adding conditions on locations being different ------------ |
---|
[288382] | 240 | for (i=1; i<=w; i++) |
---|
| 241 | { |
---|
| 242 | for (j=i+1; j<=w; j++) |
---|
| 243 | { |
---|
| 244 | crht_md=crht_md,Z(i)*Z(j)*p_poly(n,i,j); |
---|
| 245 | } |
---|
| 246 | } |
---|
[7f3ad4] | 247 | |
---|
| 248 | export crht_md; |
---|
| 249 | return(@crht_md); |
---|
| 250 | } |
---|
[fba86f8] | 251 | example |
---|
[288382] | 252 | { |
---|
| 253 | "EXAMPLE:"; echo=2; |
---|
[fba86f8] | 254 | intvec v = option(get); |
---|
[288382] | 255 | // binary cyclic [15,7,5] code with defining set (1,3) |
---|
[7f3ad4] | 256 | |
---|
| 257 | list defset=1,3; // defining set |
---|
| 258 | int n=15; // length |
---|
| 259 | int d=5; // candidate for the minimum distance |
---|
| 260 | |
---|
| 261 | def A=sysCRHTMindist(n,defset,d); |
---|
[288382] | 262 | setring A; |
---|
[fba86f8] | 263 | A; // shows the ring we are working in |
---|
| 264 | print(crht_md); // the Sala's ideal for mindist |
---|
[288382] | 265 | option(redSB); |
---|
[7f3ad4] | 266 | ideal red_crht_md=std(crht_md); |
---|
[fba86f8] | 267 | print(red_crht_md); // reduced Groebner basis |
---|
[7f3ad4] | 268 | |
---|
[fba86f8] | 269 | option(set,v); |
---|
[288382] | 270 | } |
---|
| 271 | |
---|
[6e118cf] | 272 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 273 | // slightly modified mod function |
---|
[288382] | 274 | static proc mod_ (int n, int m) |
---|
| 275 | { |
---|
| 276 | if (n mod m==0) {return(m);} |
---|
| 277 | if (n mod m!=0) {return(n mod m);} |
---|
| 278 | } |
---|
| 279 | |
---|
[6e118cf] | 280 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 281 | |
---|
[fba86f8] | 282 | proc sysNewton (int n, list defset, int t, int q, int m, list #) |
---|
[7f3ad4] | 283 | "USAGE: sysNewton (n,defset,t,q,m,[tr]); n,t,q,m,tr int, defset list int's |
---|
[fba86f8] | 284 | @format |
---|
[7f3ad4] | 285 | - n is length, |
---|
[fba86f8] | 286 | - defset is the defining set, |
---|
[7f3ad4] | 287 | - t is the number of errors, |
---|
| 288 | - q is basefield size, |
---|
[fba86f8] | 289 | - m is degree extension of the splitting field, |
---|
[7f3ad4] | 290 | - if tr>0 it indicates that Newton identities in triangular |
---|
[fba86f8] | 291 | form should be constructed |
---|
| 292 | @end format |
---|
[7f3ad4] | 293 | RETURN: the ring to work with the generalized Newton identities (in |
---|
[3599e9] | 294 | triangular form if applicable) containing the ideal with name 'newton' |
---|
[7f3ad4] | 295 | THEORY: Based on 'defset' of the given cyclic code, the procedure constructs |
---|
| 296 | the corresponding ideal 'newton' with the generalized Newton |
---|
| 297 | identities. With its help one can solve the decoding problem. For |
---|
[24f6cd9] | 298 | basics of the method @ref{Cooper philosophy}. |
---|
[3599e9] | 299 | SEE ALSO: sysCRHT, sysBin |
---|
[288382] | 300 | EXAMPLE: example sysNewton; shows an example |
---|
| 301 | " |
---|
[7f3ad4] | 302 | { |
---|
[288382] | 303 | string s="ring @newton=("+string(q)+",a),("; |
---|
| 304 | int i,j; |
---|
| 305 | int flag; |
---|
[fba86f8] | 306 | int tr; |
---|
[7f3ad4] | 307 | |
---|
[fba86f8] | 308 | if (size(#)>0) |
---|
| 309 | { |
---|
| 310 | tr=#[1]; |
---|
| 311 | } |
---|
[7f3ad4] | 312 | |
---|
[288382] | 313 | for (i=n; i>=1; i--) |
---|
| 314 | { |
---|
| 315 | for (j=1; j<=size(defset); j++) |
---|
| 316 | { |
---|
| 317 | flag=1; |
---|
| 318 | if (i==defset[j]) |
---|
| 319 | { |
---|
| 320 | flag=0; |
---|
[7f3ad4] | 321 | break; |
---|
[288382] | 322 | } |
---|
| 323 | } |
---|
| 324 | if (flag) |
---|
| 325 | { |
---|
| 326 | s=s+"S("+string(i)+"),"; |
---|
| 327 | } |
---|
| 328 | } |
---|
| 329 | s=s+"sigma(1.."+string(t)+"),"; |
---|
| 330 | for (i=size(defset); i>=2; i--) |
---|
| 331 | { |
---|
| 332 | s=s+"S("+string(defset[i])+"),"; |
---|
| 333 | } |
---|
[7f3ad4] | 334 | s=s+"S("+string(defset[1])+")),lp;"; |
---|
| 335 | |
---|
[288382] | 336 | execute(s); |
---|
[7f3ad4] | 337 | |
---|
| 338 | ideal newton; |
---|
[288382] | 339 | poly sum; |
---|
[7f3ad4] | 340 | |
---|
| 341 | |
---|
[6e118cf] | 342 | //------------ generate generalized Newton identities ---------- |
---|
[fba86f8] | 343 | if (tr) |
---|
[288382] | 344 | { |
---|
| 345 | for (i=1; i<=t; i++) |
---|
| 346 | { |
---|
| 347 | sum=0; |
---|
| 348 | for (j=1; j<=i-1; j++) |
---|
| 349 | { |
---|
| 350 | sum=sum+sigma(j)*S(i-j); |
---|
| 351 | } |
---|
| 352 | newton=newton,S(i)+sum+number(i)*sigma(i); |
---|
| 353 | } |
---|
| 354 | } else |
---|
[7f3ad4] | 355 | { |
---|
[288382] | 356 | for (i=1; i<=t; i++) |
---|
| 357 | { |
---|
| 358 | sum=0; |
---|
| 359 | for (j=1; j<=t; j++) |
---|
| 360 | { |
---|
| 361 | sum=sum+sigma(j)*S(mod_(i-j,n)); |
---|
| 362 | } |
---|
| 363 | newton=newton,S(i)+sum; |
---|
| 364 | } |
---|
| 365 | } |
---|
| 366 | for (i=1; i<=n-t; i++) |
---|
| 367 | { |
---|
| 368 | sum=0; |
---|
| 369 | for (j=1; j<=t; j++) |
---|
| 370 | { |
---|
| 371 | sum=sum+sigma(j)*S(t+i-j); |
---|
| 372 | } |
---|
| 373 | newton=newton,S(t+i)+sum; |
---|
[7f3ad4] | 374 | } |
---|
| 375 | |
---|
[6e118cf] | 376 | //----------- add field equations on sigma's -------------- |
---|
[288382] | 377 | for (i=1; i<=t; i++) |
---|
| 378 | { |
---|
| 379 | newton=newton,sigma(i)^(q^m)-sigma(i); |
---|
| 380 | } |
---|
[7f3ad4] | 381 | |
---|
[6e118cf] | 382 | //----------- add conjugacy relations ------------------ |
---|
[288382] | 383 | for (i=1; i<=n; i++) |
---|
| 384 | { |
---|
| 385 | newton=newton,S(i)^q-S(mod_(q*i,n)); |
---|
| 386 | } |
---|
| 387 | newton=simplify(newton,2); |
---|
| 388 | export newton; |
---|
[7f3ad4] | 389 | return(@newton); |
---|
| 390 | } |
---|
| 391 | example |
---|
[288382] | 392 | { |
---|
| 393 | "EXAMPLE:"; echo = 2; |
---|
[7f3ad4] | 394 | // Newton identities for a binary 3-error-correcting cyclic code of |
---|
[fba86f8] | 395 | //length 31 with defining set (1,5,7) |
---|
[7f3ad4] | 396 | |
---|
[fba86f8] | 397 | int n=31; // length |
---|
[288382] | 398 | list defset=1,5,7; //defining set |
---|
[fba86f8] | 399 | int t=3; // number of errors |
---|
| 400 | int q=2; // basefield size |
---|
| 401 | int m=5; // degree extension of the splitting field |
---|
| 402 | int tr=1; // indicator of triangular form of Newton identities |
---|
[7f3ad4] | 403 | |
---|
| 404 | |
---|
[288382] | 405 | def A=sysNewton(n,defset,t,q,m); |
---|
| 406 | setring A; |
---|
[fba86f8] | 407 | A; // shows the ring we are working in |
---|
| 408 | print(newton); // generalized Newton identities |
---|
[7f3ad4] | 409 | |
---|
[288382] | 410 | //=============================== |
---|
[fba86f8] | 411 | A=sysNewton(n,defset,t,q,m,tr); |
---|
[7f3ad4] | 412 | setring A; |
---|
[fba86f8] | 413 | print(newton); // generalized Newton identities in triangular form |
---|
[288382] | 414 | } |
---|
| 415 | |
---|
[6e118cf] | 416 | /////////////////////////////////////////////////////////////////////////////// |
---|
[7f3ad4] | 417 | // forms a list of special combinations needed for computation of Waring's |
---|
[fba86f8] | 418 | //function |
---|
[288382] | 419 | static proc combinations_sum (int m, int n) |
---|
| 420 | { |
---|
| 421 | list result; |
---|
| 422 | list comb=combinations(m-1,n+m-1); |
---|
| 423 | int i,j,flag,count; |
---|
| 424 | list interm=comb; |
---|
| 425 | for (i=1; i<=size(comb); i++) |
---|
| 426 | { |
---|
| 427 | interm[i][1]=comb[i][1]-1; |
---|
| 428 | for (j=2; j<=m-1; j++) |
---|
| 429 | { |
---|
| 430 | interm[i][j]=comb[i][j]-comb[i][j-1]-1; |
---|
| 431 | } |
---|
| 432 | interm[i][m]=n+m-comb[i][m-1]-1; |
---|
| 433 | flag=1; |
---|
| 434 | count=2; |
---|
| 435 | while ((flag)&&(count<=m)) |
---|
| 436 | { |
---|
| 437 | if (interm[i][count] mod count != 0) {flag=0;} |
---|
| 438 | count++; |
---|
| 439 | } |
---|
[7f3ad4] | 440 | if (flag) |
---|
[288382] | 441 | { |
---|
| 442 | for (j=2; j<=m; j++) |
---|
| 443 | { |
---|
| 444 | interm[i][j]=interm[i][j] div j; |
---|
| 445 | } |
---|
[7f3ad4] | 446 | result[size(result)+1]=interm[i]; |
---|
[288382] | 447 | } |
---|
| 448 | } |
---|
| 449 | return(result); |
---|
| 450 | } |
---|
| 451 | |
---|
[6e118cf] | 452 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 453 | //if n=q^e*m, m and n are coprime, then return e |
---|
[288382] | 454 | static proc exp_count (int n, int q) |
---|
| 455 | { |
---|
| 456 | int flag=1; |
---|
| 457 | int result=0; |
---|
| 458 | while(flag) |
---|
| 459 | { |
---|
| 460 | if (n mod q != 0) {flag=0;} |
---|
| 461 | else {n=n div q; result++;} |
---|
| 462 | } |
---|
| 463 | return(result); |
---|
| 464 | } |
---|
| 465 | |
---|
[6e118cf] | 466 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 467 | |
---|
[288382] | 468 | |
---|
[fba86f8] | 469 | proc sysBin (int v, list Q, int n, list #) |
---|
| 470 | "USAGE: sysBin (v,Q,n,[odd]); v,n,odd are int, Q is list of int's |
---|
| 471 | @format |
---|
[7f3ad4] | 472 | - v a number if errors, |
---|
| 473 | - Q is a generating set of the code, |
---|
| 474 | - n the length, |
---|
| 475 | - odd is an additional parameter: if |
---|
| 476 | set to 1, then the generating set is enlarged by odd elements, |
---|
[fba86f8] | 477 | which are 2^(some power)*(some elment in the gen.set) mod n |
---|
| 478 | @end format |
---|
| 479 | RETURN: the ring with the resulting system called 'bin' |
---|
[7f3ad4] | 480 | THEORY: Based on Q of the given cyclic code, the procedure constructs |
---|
| 481 | the corresponding ideal 'bin' with the use of Waring function. |
---|
| 482 | With its help one can solve the decoding problem. |
---|
[3599e9] | 483 | For basics of the method @ref{Generalized Newton identities}. |
---|
| 484 | SEE ALSO: sysNewton, sysCRHT |
---|
[288382] | 485 | EXAMPLE: example sysBin; shows an example |
---|
| 486 | " |
---|
| 487 | { |
---|
[fba86f8] | 488 | int odd; |
---|
| 489 | if (size(#)>0) |
---|
| 490 | { |
---|
| 491 | odd=#[1]; |
---|
| 492 | } |
---|
| 493 | |
---|
[288382] | 494 | //ring r=2,(sigma(1..v),S(1..n)),(lp(v),dp(n)); |
---|
| 495 | ring r=2,(S(1..n),sigma(1..v)),lp; |
---|
| 496 | list cyclot; |
---|
| 497 | ideal result; |
---|
| 498 | int i,j,k,s; |
---|
| 499 | list comb; |
---|
| 500 | poly sum_, mon; |
---|
| 501 | int count1, count2, upper, coef_, flag, gener; |
---|
| 502 | list Q_update; |
---|
[fba86f8] | 503 | if (odd==1) |
---|
[288382] | 504 | { |
---|
| 505 | for (i=1; i<=n; i++) |
---|
| 506 | { |
---|
| 507 | cyclot[i]=0; |
---|
| 508 | } |
---|
| 509 | for (i=1; i<=size(Q); i++) |
---|
| 510 | { |
---|
| 511 | flag=1; |
---|
| 512 | gener=Q[i]; |
---|
| 513 | while(flag) |
---|
| 514 | { |
---|
| 515 | cyclot[gener]=1; |
---|
| 516 | gener=2*gener mod n; |
---|
| 517 | if (gener == Q[i]) {flag=0;} |
---|
| 518 | } |
---|
| 519 | } |
---|
| 520 | for (i=1; i<=n; i++) |
---|
| 521 | { |
---|
| 522 | if ((cyclot[i] == 1)&&(i mod 2 == 1)) {Q_update[size(Q_update)+1]=i;} |
---|
| 523 | } |
---|
| 524 | } |
---|
| 525 | else |
---|
| 526 | { |
---|
| 527 | Q_update=Q; |
---|
| 528 | } |
---|
[7f3ad4] | 529 | |
---|
| 530 | //---- form polynomials for the Bin system via Waring's function --------- |
---|
[288382] | 531 | for (i=1; i<=size(Q_update); i++) |
---|
| 532 | { |
---|
| 533 | comb=combinations_sum(v,Q_update[i]); |
---|
| 534 | sum_=0; |
---|
| 535 | for (k=1; k<=size(comb); k++) |
---|
| 536 | { |
---|
| 537 | upper=0; |
---|
| 538 | for (j=1; j<=v; j++) |
---|
| 539 | { |
---|
| 540 | upper=upper+comb[k][j]; |
---|
| 541 | } |
---|
| 542 | count1=0; |
---|
[7f3ad4] | 543 | for (j=2; j<=upper-1; j++) |
---|
[288382] | 544 | { |
---|
| 545 | count1=count1+exp_count(j,2); |
---|
| 546 | } |
---|
| 547 | count1=count1+exp_count(Q_update[i],2); |
---|
| 548 | count2=0; |
---|
| 549 | for (j=1; j<=v; j++) |
---|
| 550 | { |
---|
| 551 | for (s=2; s<=comb[k][j]; s++) |
---|
| 552 | { |
---|
| 553 | count2=count2+exp_count(s,2); |
---|
| 554 | } |
---|
| 555 | } |
---|
| 556 | if (count1<count2) {print("ERRORsysBin");} |
---|
| 557 | if (count1>count2) {coef_=0;} |
---|
| 558 | if (count1 == count2) {coef_=1;} |
---|
| 559 | mon=1; |
---|
| 560 | for (j=1; j<=v; j++) |
---|
| 561 | { |
---|
| 562 | mon=mon*sigma(j)^(comb[k][j]); |
---|
| 563 | } |
---|
| 564 | sum_=sum_+coef_*mon; |
---|
| 565 | } |
---|
[7f3ad4] | 566 | result=result,S(Q_update[i])-sum_; |
---|
[288382] | 567 | } |
---|
| 568 | ideal bin=simplify(result,2); |
---|
| 569 | export bin; |
---|
| 570 | return(r); |
---|
[7f3ad4] | 571 | } |
---|
| 572 | example |
---|
[288382] | 573 | { |
---|
| 574 | "EXAMPLE:"; echo = 2; |
---|
| 575 | // [31,16,7] quadratic residue code |
---|
| 576 | list l=1,5,7,9,19,25; |
---|
| 577 | // we do not need even synromes here |
---|
| 578 | def A=sysBin(3,l,31); |
---|
| 579 | setring A; |
---|
| 580 | print(bin); |
---|
| 581 | } |
---|
| 582 | |
---|
[6e118cf] | 583 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 584 | |
---|
[288382] | 585 | proc encode (matrix x, matrix g) |
---|
| 586 | "USAGE: encode (x, g); x a row vector (message), and g a generator matrix |
---|
| 587 | RETURN: corresponding codeword |
---|
| 588 | EXAMPLE: example encode; shows an example |
---|
| 589 | " |
---|
| 590 | { |
---|
| 591 | if (nrows(x)>1) {print("ERRORencode1!");} |
---|
| 592 | if (ncols(x)!=nrows(g)) {print("ERRORencode2!");} |
---|
| 593 | return(x*g); |
---|
[7f3ad4] | 594 | } |
---|
| 595 | example |
---|
[288382] | 596 | { |
---|
| 597 | "EXAMPLE:"; echo = 2; |
---|
| 598 | ring r=2,x,dp; |
---|
| 599 | matrix x[1][4]=1,0,1,0; |
---|
| 600 | matrix g[4][7]=1,0,0,0,0,1,1, |
---|
| 601 | 0,1,0,0,1,0,1, |
---|
| 602 | 0,0,1,0,1,1,1, |
---|
| 603 | 0,0,0,1,1,1,0; |
---|
| 604 | //encode x with the generator matrix g |
---|
[7f3ad4] | 605 | print(encode(x,g)); |
---|
[288382] | 606 | } |
---|
| 607 | |
---|
[6e118cf] | 608 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 609 | |
---|
[288382] | 610 | proc syndrome (matrix h, matrix c) |
---|
| 611 | "USAGE: syndrome (h, c); h a check matrix, c a row vector (codeword) |
---|
| 612 | RETURN: corresponding syndrome |
---|
| 613 | EXAMPLE: example syndrome; shows an example |
---|
| 614 | " |
---|
| 615 | { |
---|
| 616 | if (nrows(c)>1) {print("ERRORsyndrome1!");} |
---|
| 617 | if (ncols(c)!=ncols(h)) {print("ERRORsyndrome2!");} |
---|
[7f3ad4] | 618 | return(h*transpose(c)); |
---|
| 619 | } |
---|
| 620 | example |
---|
[288382] | 621 | { |
---|
| 622 | "EXAMPLE:"; echo = 2; |
---|
| 623 | ring r=2,x,dp; |
---|
| 624 | matrix x[1][4]=1,0,1,0; |
---|
| 625 | matrix g[4][7]=1,0,0,0,0,1,1, |
---|
| 626 | 0,1,0,0,1,0,1, |
---|
| 627 | 0,0,1,0,1,1,1, |
---|
| 628 | 0,0,0,1,1,1,0; |
---|
| 629 | //encode x with the generator matrix g |
---|
| 630 | matrix c=encode(x,g); |
---|
| 631 | // disturb |
---|
| 632 | c[1,3]=0; |
---|
| 633 | //compute syndrome |
---|
| 634 | //corresponding check matrix |
---|
| 635 | matrix check[3][7]=1,0,0,1,1,0,1,0,1,0,1,0,1,1,0,0,1,0,1,1,1; |
---|
| 636 | print(syndrome(check,c)); |
---|
| 637 | c[1,3]=1; |
---|
[7f3ad4] | 638 | //now c is a codeword |
---|
[288382] | 639 | print(syndrome(check,c)); |
---|
| 640 | } |
---|
| 641 | |
---|
[6e118cf] | 642 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 643 | // (coordinatewise) star product of two vectors |
---|
[288382] | 644 | static proc star(matrix m, int i, int j) |
---|
| 645 | { |
---|
| 646 | matrix result[ncols(m)][1]; |
---|
| 647 | for (int k=1; k<=ncols(m); k++) |
---|
| 648 | { |
---|
| 649 | result[k,1]=m[i,k]*m[j,k]; |
---|
| 650 | } |
---|
| 651 | return(result); |
---|
| 652 | } |
---|
| 653 | |
---|
[6e118cf] | 654 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 655 | |
---|
[fba86f8] | 656 | proc sysQE(matrix check, matrix y, int t, list #) |
---|
| 657 | "USAGE: sysQE(check,y,t,[fieldeq,formal]);check,y matrix;t,fieldeq,formal int |
---|
| 658 | @format |
---|
[7f3ad4] | 659 | - check is the check matrix of the code |
---|
| 660 | - y is a received word, |
---|
| 661 | - t the number of errors to be corrected, |
---|
| 662 | - if fieldeq=1, then field equations are added, |
---|
| 663 | - if formal=0, field equations on (known) syndrome variables |
---|
| 664 | are not added, in order to add them (note that the exponent should |
---|
[fba86f8] | 665 | be as a number of elements in the INITIAL alphabet) one |
---|
[288382] | 666 | needs to set formal>0 for the exponent |
---|
[fba86f8] | 667 | @end format |
---|
| 668 | RETURN: the ring to work with together with the resulting system called 'qe' |
---|
[7f3ad4] | 669 | THEORY: Based on 'check' of the given linear code, the procedure constructs |
---|
[3599e9] | 670 | the corresponding ideal that gives an opportunity to compute |
---|
| 671 | unknown syndrome of the received word y. Further |
---|
[7f3ad4] | 672 | one is able to solve the decoding problem. |
---|
[3599e9] | 673 | For basics of the method @ref{Decoding method based on quadratic equations}. |
---|
| 674 | SEE ALSO: sysFL |
---|
[288382] | 675 | EXAMPLE: example sysQE; shows an example |
---|
| 676 | " |
---|
| 677 | { |
---|
[fba86f8] | 678 | int fieldeq; |
---|
| 679 | int formal; |
---|
| 680 | if (size(#)>0) |
---|
| 681 | { |
---|
| 682 | fieldeq=#[1]; |
---|
| 683 | } |
---|
| 684 | if (size(#)>1) |
---|
| 685 | { |
---|
| 686 | formal=#[2]; |
---|
| 687 | } |
---|
[7f3ad4] | 688 | |
---|
| 689 | def br=basering; |
---|
[288382] | 690 | list rl=ringlist(br); |
---|
[7f3ad4] | 691 | |
---|
[288382] | 692 | int red=nrows(check); |
---|
| 693 | int n=ncols(check); |
---|
| 694 | int q=rl[1][1]; |
---|
[7f3ad4] | 695 | |
---|
[288382] | 696 | if (formal==0) |
---|
[7f3ad4] | 697 | { |
---|
[288382] | 698 | ring work=(q,a),(V(1..t),U(1..n)),dp; |
---|
| 699 | } else |
---|
| 700 | { |
---|
| 701 | ring work=(q,a),(V(1..t),U(1..n),s(1..red)),(dp(t),lp(n),dp(red)); |
---|
| 702 | } |
---|
[7f3ad4] | 703 | |
---|
[288382] | 704 | matrix check=imap(br,check); |
---|
| 705 | matrix y=imap(br,y); |
---|
[7f3ad4] | 706 | |
---|
[288382] | 707 | matrix h_full=genMDSMat(n,a); |
---|
| 708 | matrix h=submat(h_full,1..red,1..n); |
---|
| 709 | if (nrows(y)!=1) {print("ERROR1Pell");} |
---|
| 710 | if (ncols(y)!=n) {print("ERROR2Pell");} |
---|
[7f3ad4] | 711 | |
---|
[288382] | 712 | ideal result; |
---|
[7f3ad4] | 713 | |
---|
[288382] | 714 | list c; |
---|
| 715 | list a; |
---|
| 716 | list tmp,tmp2; |
---|
| 717 | int i,j,l,k; |
---|
| 718 | number sum,prod,sig; |
---|
| 719 | poly sum1,sum2,sum3; |
---|
| 720 | for (i=1; i<=n; i++) |
---|
| 721 | { |
---|
| 722 | c[i]=tmp; |
---|
| 723 | } |
---|
[7f3ad4] | 724 | |
---|
[288382] | 725 | int tim=rtimer; |
---|
| 726 | matrix transf=inverse(transpose(h_full)); |
---|
[7f3ad4] | 727 | |
---|
[fba86f8] | 728 | //------ expression matrix of check vectors w.r.t. the MDS basis ----------- |
---|
[288382] | 729 | tim=rtimer; |
---|
| 730 | for (i=1; i<=red ; i++) |
---|
| 731 | { |
---|
| 732 | a[i]=transpose(submat(check,i..i,1..n)); |
---|
| 733 | a[i]=transf*a[i]; |
---|
| 734 | } |
---|
[7f3ad4] | 735 | |
---|
[6e118cf] | 736 | //----------- compute the structure constants ------------------------ |
---|
[288382] | 737 | tim=rtimer; |
---|
[7f3ad4] | 738 | matrix te[n][1]; |
---|
[288382] | 739 | for (i=1; i<=n; i++) |
---|
| 740 | { |
---|
| 741 | for (j=1; j<=t+1; j++) |
---|
| 742 | { |
---|
| 743 | if ((j<i)&&(i<=t+1)) {c[i][j]=c[j][i];} |
---|
[7f3ad4] | 744 | else |
---|
[288382] | 745 | { |
---|
[7f3ad4] | 746 | if (i+j<=n+1) |
---|
[288382] | 747 | { |
---|
| 748 | c[i][j]=te; |
---|
[7f3ad4] | 749 | c[i][j][i+j-1,1]=1; |
---|
[288382] | 750 | } |
---|
| 751 | else |
---|
| 752 | { |
---|
| 753 | c[i][j]=star(h_full,i,j); |
---|
| 754 | c[i][j]=transf*c[i][j]; |
---|
[7f3ad4] | 755 | } |
---|
[288382] | 756 | } |
---|
| 757 | } |
---|
| 758 | } |
---|
[7f3ad4] | 759 | |
---|
| 760 | |
---|
| 761 | tim=rtimer; |
---|
[288382] | 762 | if (formal==0) |
---|
| 763 | { |
---|
| 764 | matrix s[red][1]=syndrome(check,y); |
---|
| 765 | for (j=1; j<=red; j++) |
---|
| 766 | { |
---|
| 767 | sum1=0; |
---|
| 768 | for (l=1; l<=n; l++) |
---|
| 769 | { |
---|
| 770 | sum1=sum1+a[j][l,1]*U(l); |
---|
| 771 | } |
---|
| 772 | result=result,sum1-s[j,1]; |
---|
| 773 | } |
---|
| 774 | } else |
---|
| 775 | { |
---|
| 776 | for (j=1; j<=red; j++) |
---|
| 777 | { |
---|
| 778 | sum1=0; |
---|
| 779 | for (l=1; l<=n; l++) |
---|
| 780 | { |
---|
| 781 | sum1=sum1+a[j][l,1]*U(l); |
---|
| 782 | } |
---|
| 783 | result=result,sum1-s(j); |
---|
| 784 | } |
---|
| 785 | for (j=1; j<=red; j++) |
---|
| 786 | { |
---|
| 787 | result=result,s(j)^(formal)-s(j); |
---|
| 788 | } |
---|
| 789 | } |
---|
| 790 | if (fieldeq) |
---|
| 791 | { |
---|
| 792 | for (i=1; i<=n; i++) |
---|
| 793 | { |
---|
| 794 | result=result,U(i)^q-U(i); |
---|
| 795 | } |
---|
| 796 | for (j=1; j<=t; j++) |
---|
| 797 | { |
---|
| 798 | result=result,V(j)^q-V(j); |
---|
| 799 | } |
---|
[7f3ad4] | 800 | } |
---|
| 801 | |
---|
[fba86f8] | 802 | //----- form the quadratic equations according to the theory ----------- |
---|
[288382] | 803 | for (i=1; i<=n; i++) |
---|
| 804 | { |
---|
| 805 | sum1=0; |
---|
| 806 | for (j=1; j<=t; j++) |
---|
| 807 | { |
---|
| 808 | sum2=0; |
---|
| 809 | for (l=1; l<=n; l++) |
---|
| 810 | { |
---|
| 811 | sum2=sum2+c[i][j][l,1]*U(l); |
---|
| 812 | } |
---|
| 813 | sum1=sum1+sum2*V(j); |
---|
| 814 | } |
---|
| 815 | sum3=0; |
---|
| 816 | for (l=1; l<=n; l++) |
---|
| 817 | { |
---|
| 818 | sum3=sum3+c[i][t+1][l,1]*U(l); |
---|
| 819 | } |
---|
| 820 | result=result,sum1-sum3; |
---|
[7f3ad4] | 821 | } |
---|
| 822 | |
---|
[288382] | 823 | result=simplify(result,2); |
---|
[7f3ad4] | 824 | |
---|
[288382] | 825 | ideal qe=result; |
---|
| 826 | export qe; |
---|
[7f3ad4] | 827 | return(work); |
---|
| 828 | } |
---|
| 829 | example |
---|
[288382] | 830 | { |
---|
| 831 | "EXAMPLE:"; echo = 2; |
---|
[fba86f8] | 832 | intvec v = option(get); |
---|
[7f3ad4] | 833 | |
---|
[288382] | 834 | //correct 2 errors in [7,3] 8-ary code RS code |
---|
| 835 | int t=2; int q=8; int n=7; int redun=4; |
---|
| 836 | ring r=(q,a),x,dp; |
---|
| 837 | matrix h_full=genMDSMat(n,a); |
---|
| 838 | matrix h=submat(h_full,1..redun,1..n); |
---|
| 839 | matrix g=dual_code(h); |
---|
| 840 | matrix x[1][3]=0,0,1,0; |
---|
| 841 | matrix y[1][7]=encode(x,g); |
---|
[7f3ad4] | 842 | |
---|
[288382] | 843 | //disturb with 2 errors |
---|
[7f3ad4] | 844 | matrix rec[1][7]=errorInsert(y,list(2,4),list(1,a)); |
---|
| 845 | |
---|
[288382] | 846 | //generate the system |
---|
[fba86f8] | 847 | def A=sysQE(h,rec,t); |
---|
[288382] | 848 | setring A; |
---|
| 849 | print(qe); |
---|
[7f3ad4] | 850 | |
---|
[288382] | 851 | //let us decode |
---|
| 852 | option(redSB); |
---|
| 853 | ideal sys_qe=std(qe); |
---|
[7f3ad4] | 854 | print(sys_qe); |
---|
| 855 | |
---|
[fba86f8] | 856 | option(set,v); |
---|
[288382] | 857 | } |
---|
| 858 | |
---|
[6e118cf] | 859 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 860 | |
---|
[00142a] | 861 | proc errorInsert(matrix y, list pos, list val) |
---|
[fba86f8] | 862 | "USAGE: errorInsert(y,pos,val); y is matrix, pos,val list of int's |
---|
| 863 | @format |
---|
[7f3ad4] | 864 | - y is a (code) word, |
---|
| 865 | - pos = positions where errors occured, |
---|
[fba86f8] | 866 | - val = their corresponding values |
---|
| 867 | @end format |
---|
[288382] | 868 | RETURN: corresponding received word |
---|
[fba86f8] | 869 | EXAMPLE: example errorInsert; shows an example |
---|
[288382] | 870 | " |
---|
| 871 | { |
---|
| 872 | matrix result[1][ncols(y)]=y; |
---|
| 873 | if (size(pos)!=size(val)) {print("ERRORerror");} |
---|
| 874 | for (int i=1; i<=size(pos); i++) |
---|
| 875 | { |
---|
| 876 | result[1,pos[i]]=y[1,pos[i]]+val[i]; |
---|
| 877 | } |
---|
| 878 | return(result); |
---|
[7f3ad4] | 879 | } |
---|
| 880 | example |
---|
[288382] | 881 | { |
---|
| 882 | "EXAMPLE:"; echo = 2; |
---|
| 883 | //correct 2 errors in [7,3] 8-ary code RS code |
---|
| 884 | int t=2; int q=8; int n=7; int redun=4; |
---|
| 885 | ring r=(q,a),x,dp; |
---|
| 886 | matrix h_full=genMDSMat(n,a); |
---|
| 887 | matrix h=submat(h_full,1..redun,1..n); |
---|
| 888 | matrix g=dual_code(h); |
---|
| 889 | matrix x[1][3]=0,0,1,0; |
---|
| 890 | matrix y[1][7]=encode(x,g); |
---|
| 891 | print(y); |
---|
[7f3ad4] | 892 | |
---|
[288382] | 893 | //disturb with 2 errors |
---|
[00142a] | 894 | matrix rec[1][7]=errorInsert(y,list(2,4),list(1,a)); |
---|
[288382] | 895 | print(rec); |
---|
| 896 | print(rec-y); |
---|
| 897 | } |
---|
| 898 | |
---|
[6e118cf] | 899 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 900 | |
---|
[288382] | 901 | proc errorRand(matrix y, int num, int e) |
---|
[fba86f8] | 902 | "USAGE: errorRand(y, num, e); y matrix, num,e int |
---|
| 903 | @format |
---|
[7f3ad4] | 904 | - y is a (code) word, |
---|
| 905 | - num is the number of errors, |
---|
[fba86f8] | 906 | - e is an extension degree (if one wants values to be from GF(p^e)) |
---|
| 907 | @end format |
---|
[288382] | 908 | RETURN: corresponding received word |
---|
| 909 | EXAMPLE: example errorRand; shows an example |
---|
| 910 | " |
---|
| 911 | { |
---|
| 912 | matrix result[1][ncols(y)]=y; |
---|
| 913 | int i,j, flag, temp; |
---|
| 914 | list pos, val; |
---|
| 915 | matrix tempnum; |
---|
| 916 | |
---|
| 917 | for (i=1; i<=num; i++) |
---|
| 918 | { |
---|
| 919 | while(1) |
---|
| 920 | { |
---|
| 921 | temp=random(1,ncols(y)); |
---|
| 922 | flag=1; |
---|
| 923 | for (j=1; j<=size(pos); j++) |
---|
| 924 | { |
---|
| 925 | if (temp==pos[j]) {flag=0;} |
---|
| 926 | } |
---|
| 927 | if (flag) {pos[i]=temp;break;} |
---|
[7f3ad4] | 928 | } |
---|
[288382] | 929 | } |
---|
[7f3ad4] | 930 | |
---|
[288382] | 931 | for (i=1; i<=num; i++) |
---|
| 932 | { |
---|
| 933 | flag=1; |
---|
| 934 | while(flag) |
---|
[7f3ad4] | 935 | { |
---|
| 936 | tempnum=randomvector(1,e); |
---|
| 937 | if (tempnum!=0) {flag=0;} |
---|
[288382] | 938 | } |
---|
[7f3ad4] | 939 | val[i]=tempnum; |
---|
[288382] | 940 | } |
---|
[7f3ad4] | 941 | |
---|
[288382] | 942 | for (i=1; i<=size(pos); i++) |
---|
| 943 | { |
---|
| 944 | result[1,pos[i]]=y[1,pos[i]]+val[i]; |
---|
| 945 | } |
---|
| 946 | return(result); |
---|
[7f3ad4] | 947 | } |
---|
| 948 | example |
---|
[288382] | 949 | { |
---|
| 950 | "EXAMPLE:"; echo = 2; |
---|
| 951 | //correct 2 errors in [7,3] 8-ary code RS code |
---|
| 952 | int t=2; int q=8; int n=7; int redun=4; |
---|
| 953 | ring r=(q,a),x,dp; |
---|
| 954 | matrix h_full=genMDSMat(n,a); |
---|
| 955 | matrix h=submat(h_full,1..redun,1..n); |
---|
| 956 | matrix g=dual_code(h); |
---|
| 957 | matrix x[1][3]=0,0,1,0; |
---|
| 958 | matrix y[1][7]=encode(x,g); |
---|
[7f3ad4] | 959 | |
---|
[288382] | 960 | //disturb with 2 random errors |
---|
| 961 | matrix rec[1][7]=errorRand(y,2,3); |
---|
| 962 | print(rec); |
---|
[7f3ad4] | 963 | print(rec-y); |
---|
[288382] | 964 | } |
---|
| 965 | |
---|
[6e118cf] | 966 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 967 | |
---|
[fba86f8] | 968 | proc randomCheck(int m, int n, int e) |
---|
| 969 | "USAGE: randomCheck(m, n, e); m,n,e are int |
---|
| 970 | @format |
---|
[7f3ad4] | 971 | - m x n are dimensions of the matrix, |
---|
[fba86f8] | 972 | - e is an extension degree (if one wants values to be from GF(p^e)) |
---|
| 973 | @end format |
---|
[288382] | 974 | RETURN: random check matrix |
---|
| 975 | EXAMPLE: example randomCheck; shows an example |
---|
| 976 | " |
---|
| 977 | { |
---|
| 978 | matrix result[m][n]; |
---|
| 979 | matrix rand[m][n-m]; |
---|
| 980 | int i,j; |
---|
| 981 | matrix temp; |
---|
| 982 | for (i=1; i<=m; i++) |
---|
| 983 | { |
---|
[fba86f8] | 984 | temp=randomvector(n-m,e); |
---|
[288382] | 985 | for (j=1; j<=n-m; j++) |
---|
| 986 | { |
---|
| 987 | rand[i,j]=temp[j,1]; |
---|
[7f3ad4] | 988 | } |
---|
[288382] | 989 | } |
---|
| 990 | result=concat(rand,unitmat(m)); |
---|
| 991 | return(result); |
---|
[7f3ad4] | 992 | } |
---|
| 993 | example |
---|
[288382] | 994 | { |
---|
[7f3ad4] | 995 | "EXAMPLE:"; echo = 2; |
---|
[288382] | 996 | int redun=5; int n=15; |
---|
| 997 | ring r=2,x,dp; |
---|
[7f3ad4] | 998 | |
---|
[288382] | 999 | //generate random check matrix for a [15,5] binary code |
---|
| 1000 | matrix h=randomCheck(redun,n,1); |
---|
| 1001 | print(h); |
---|
[7f3ad4] | 1002 | |
---|
[288382] | 1003 | //corresponding generator matrix |
---|
| 1004 | matrix g=dual_code(h); |
---|
[7f3ad4] | 1005 | print(g); |
---|
[288382] | 1006 | } |
---|
| 1007 | |
---|
[6e118cf] | 1008 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1009 | |
---|
[288382] | 1010 | proc genMDSMat(int n, number a) |
---|
[fba86f8] | 1011 | "USAGE: genMDSMat(n, a); n is int, a is number |
---|
| 1012 | @format |
---|
[7f3ad4] | 1013 | - n x n are dimensions of the MDS matrix, |
---|
[fba86f8] | 1014 | - a is a primitive element of the field. |
---|
[7f3ad4] | 1015 | @end format |
---|
| 1016 | NOTE: An MDS matrix is constructed in the following way. We take a to be a |
---|
| 1017 | generator of the multiplicative group of the field. Then we construct |
---|
[fba86f8] | 1018 | the Vandermonde matrix with this a. |
---|
[7f3ad4] | 1019 | ASSUME: extension field should already be defined |
---|
[288382] | 1020 | RETURN: a matrix with the MDS property |
---|
[7f3ad4] | 1021 | EXAMPLE: example genMDSMat; shows an example |
---|
[288382] | 1022 | " |
---|
| 1023 | { |
---|
| 1024 | int i,j; |
---|
| 1025 | matrix result[n][n]; |
---|
| 1026 | for (i=0; i<=n-1; i++) |
---|
| 1027 | { |
---|
| 1028 | for (j=0; j<=n-1; j++) |
---|
| 1029 | { |
---|
| 1030 | result[j+1,i+1]=(a^i)^j; |
---|
[7f3ad4] | 1031 | } |
---|
[288382] | 1032 | } |
---|
| 1033 | return(result); |
---|
[7f3ad4] | 1034 | } |
---|
| 1035 | example |
---|
[288382] | 1036 | { |
---|
| 1037 | "EXAMPLE:"; echo = 2; |
---|
| 1038 | int q=16; int n=15; |
---|
| 1039 | ring r=(q,a),x,dp; |
---|
[7f3ad4] | 1040 | |
---|
[288382] | 1041 | //generate an MDS (Vandermonde) matrix |
---|
| 1042 | matrix h_full=genMDSMat(n,a); |
---|
[7f3ad4] | 1043 | print(h_full); |
---|
[288382] | 1044 | } |
---|
| 1045 | |
---|
[6e118cf] | 1046 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1047 | |
---|
[288382] | 1048 | |
---|
| 1049 | proc mindist (matrix check) |
---|
[fba86f8] | 1050 | "USAGE: mindist (check, q); check matrix, q int |
---|
| 1051 | @format |
---|
[7f3ad4] | 1052 | - check is a check matrix, |
---|
[fba86f8] | 1053 | - q is the field size |
---|
| 1054 | @end format |
---|
[7f3ad4] | 1055 | RETURN: minimum distance of the code together with the time needed for its |
---|
[3599e9] | 1056 | calculation |
---|
[7f3ad4] | 1057 | EXAMPLE: example mindist; shows an example |
---|
[288382] | 1058 | " |
---|
| 1059 | { |
---|
[fba86f8] | 1060 | intvec vopt = option(get); |
---|
| 1061 | |
---|
[288382] | 1062 | int n=ncols(check); int redun=nrows(check); int t=redun+1; |
---|
[7f3ad4] | 1063 | |
---|
| 1064 | def br=basering; |
---|
| 1065 | list rl=ringlist(br); |
---|
[288382] | 1066 | int q=rl[1][1]; |
---|
[7f3ad4] | 1067 | |
---|
| 1068 | ring work=(q,a),(V(1..t),U(1..n)),dp; |
---|
| 1069 | matrix check=imap(br,check); |
---|
| 1070 | |
---|
[288382] | 1071 | ideal temp; |
---|
| 1072 | int count=1; |
---|
| 1073 | int flag=1; |
---|
| 1074 | int flag2; |
---|
| 1075 | int i, tim, timsolve; |
---|
[7f3ad4] | 1076 | matrix z[1][n]; |
---|
[288382] | 1077 | option(redSB); |
---|
[fba86f8] | 1078 | def A=sysQE(check,z,count); |
---|
[6e118cf] | 1079 | |
---|
[7f3ad4] | 1080 | //proceed with solving the system w.r.t zero vector until some solutions |
---|
[fba86f8] | 1081 | //are found |
---|
[288382] | 1082 | while (flag) |
---|
| 1083 | { |
---|
[fba86f8] | 1084 | A=sysQE(check,z,count); |
---|
[288382] | 1085 | setring A; |
---|
| 1086 | ideal temp=qe; |
---|
| 1087 | tim=rtimer; |
---|
| 1088 | temp=std(temp); |
---|
[7f3ad4] | 1089 | timsolve=timsolve+rtimer-tim; |
---|
[288382] | 1090 | flag2=1; |
---|
| 1091 | setring work; |
---|
[7f3ad4] | 1092 | temp=imap(A,temp); |
---|
[288382] | 1093 | for (i=1; i<=n; i++) |
---|
| 1094 | { |
---|
[7f3ad4] | 1095 | if |
---|
| 1096 | (temp[i]!=U(n-i+1)) |
---|
[288382] | 1097 | { |
---|
| 1098 | flag2=0; |
---|
| 1099 | } |
---|
| 1100 | } |
---|
[7f3ad4] | 1101 | if (!flag2) |
---|
[288382] | 1102 | { |
---|
| 1103 | flag=0; |
---|
| 1104 | } |
---|
[7f3ad4] | 1105 | else |
---|
[288382] | 1106 | { |
---|
| 1107 | count++; |
---|
| 1108 | } |
---|
| 1109 | } |
---|
[7f3ad4] | 1110 | list result=list(count,timsolve); |
---|
| 1111 | |
---|
[fba86f8] | 1112 | option(set,vopt); |
---|
[7f3ad4] | 1113 | return(result); |
---|
| 1114 | } |
---|
| 1115 | example |
---|
[288382] | 1116 | { |
---|
| 1117 | "EXAMPLE:"; echo = 2; |
---|
[7f3ad4] | 1118 | //determine a minimum distance for a [7,3] binary code |
---|
| 1119 | int q=8; int n=7; int redun=4; int t=redun+1; |
---|
| 1120 | ring r=(q,a),x,dp; |
---|
| 1121 | |
---|
[288382] | 1122 | //generate random check matrix |
---|
| 1123 | matrix h=randomCheck(redun,n,1); |
---|
| 1124 | print(h); |
---|
[7f3ad4] | 1125 | list l=mindist(h); |
---|
| 1126 | print(l[1]); |
---|
| 1127 | //time for the comutation in secs |
---|
| 1128 | print(l[2]); |
---|
[288382] | 1129 | } |
---|
| 1130 | |
---|
[6e118cf] | 1131 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1132 | |
---|
[288382] | 1133 | proc decode(matrix check, matrix rec) |
---|
[fba86f8] | 1134 | "USAGE: decode(check, rec, t); check, rec matrix, t int |
---|
| 1135 | @format |
---|
| 1136 | - check is the check matrix of the code, |
---|
[7f3ad4] | 1137 | - rec is a received word, |
---|
[fba86f8] | 1138 | - t is an upper bound for the number of errors one wants to correct |
---|
| 1139 | @end format |
---|
[7f3ad4] | 1140 | ASSUME: Errors in rec should be correctable, otherwise the output is |
---|
[fba86f8] | 1141 | unpredictable |
---|
[3599e9] | 1142 | RETURN: a codeword that is closest to rec |
---|
[7f3ad4] | 1143 | EXAMPLE: example decode; shows an example |
---|
[288382] | 1144 | " |
---|
[fba86f8] | 1145 | { |
---|
[7f3ad4] | 1146 | intvec vopt = option(get); |
---|
| 1147 | |
---|
[288382] | 1148 | def br=basering; |
---|
[7f3ad4] | 1149 | int n=ncols(check); |
---|
| 1150 | |
---|
| 1151 | int count=1; |
---|
[fba86f8] | 1152 | def A=sysQE(check,rec,count); |
---|
[288382] | 1153 | while(1) |
---|
| 1154 | { |
---|
[fba86f8] | 1155 | A=sysQE(check,rec,count); |
---|
[7f3ad4] | 1156 | setring A; |
---|
[288382] | 1157 | matrix h_full=genMDSMat(n,a); |
---|
| 1158 | matrix rec=imap(br,rec); |
---|
| 1159 | option(redSB); |
---|
[7f3ad4] | 1160 | ideal qe_red=std(qe); |
---|
[288382] | 1161 | if (qe_red[1]!=1) |
---|
| 1162 | { |
---|
| 1163 | break; |
---|
| 1164 | } |
---|
[7f3ad4] | 1165 | else |
---|
[288382] | 1166 | { |
---|
| 1167 | count++; |
---|
| 1168 | } |
---|
| 1169 | setring br; |
---|
| 1170 | } |
---|
[7f3ad4] | 1171 | |
---|
| 1172 | setring A; |
---|
| 1173 | |
---|
[288382] | 1174 | //obtain a codeword |
---|
| 1175 | //this works only if our code is indeed can correct these errors |
---|
| 1176 | matrix syn[n][1]; |
---|
| 1177 | for (int i=1; i<=n; i++) |
---|
| 1178 | { |
---|
| 1179 | syn[i,1]=-qe_red[n-i+1]+lead(qe_red[n-i+1]); |
---|
[7f3ad4] | 1180 | } |
---|
| 1181 | |
---|
| 1182 | matrix real_syn=inverse(h_full)*syn; |
---|
[288382] | 1183 | setring br; |
---|
| 1184 | matrix real_syn=imap(A,real_syn); |
---|
[7f3ad4] | 1185 | |
---|
[fba86f8] | 1186 | option(set,vopt); |
---|
[7f3ad4] | 1187 | return(rec-transpose(real_syn)); |
---|
| 1188 | } |
---|
| 1189 | example |
---|
[288382] | 1190 | { |
---|
[7f3ad4] | 1191 | "EXAMPLE:"; echo = 2; |
---|
| 1192 | //correct 1 error in [15,7] binary code |
---|
[288382] | 1193 | int t=1; int q=16; int n=15; int redun=10; |
---|
[7f3ad4] | 1194 | ring r=(q,a),x,dp; |
---|
| 1195 | |
---|
[288382] | 1196 | //generate random check matrix |
---|
[7f3ad4] | 1197 | matrix h=randomCheck(redun,n,1); |
---|
[288382] | 1198 | matrix g=dual_code(h); |
---|
| 1199 | matrix x[1][n-redun]=0,0,1,0,1,0,1; |
---|
| 1200 | matrix y[1][n]=encode(x,g); |
---|
| 1201 | print(y); |
---|
[7f3ad4] | 1202 | |
---|
[288382] | 1203 | // find out the minimum distance of the code |
---|
| 1204 | list l=mindist(h); |
---|
[7f3ad4] | 1205 | |
---|
[288382] | 1206 | //disturb with errors |
---|
| 1207 | "Correct ",(l[1]-1) div 2," errors"; |
---|
| 1208 | matrix rec[1][n]=errorRand(y,(l[1]-1) div 2,1); |
---|
[7f3ad4] | 1209 | print(rec); |
---|
| 1210 | |
---|
[288382] | 1211 | //let us decode |
---|
| 1212 | matrix dec_word=decode(h,rec); |
---|
[7f3ad4] | 1213 | print(dec_word); |
---|
[288382] | 1214 | } |
---|
| 1215 | |
---|
[6e118cf] | 1216 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1217 | |
---|
[288382] | 1218 | |
---|
[fba86f8] | 1219 | proc decodeRandom(int n, int redun, int ncodes, int ntrials, list #) |
---|
| 1220 | "USAGE: decodeRandom(redun,q,ncodes,ntrials,[e]); all parameters int |
---|
| 1221 | @format |
---|
| 1222 | - redun is a redundabcy of a (random) code, |
---|
[7f3ad4] | 1223 | - q is the field size, |
---|
[fba86f8] | 1224 | - ncodes is the number of random codes to be processed, |
---|
| 1225 | - ntrials is the number of received vectors per code to be corrected |
---|
[7f3ad4] | 1226 | - If e is given it sets the correction capacity explicitly. It |
---|
[fba86f8] | 1227 | should be used in case one expects some lower bound, |
---|
[7f3ad4] | 1228 | otherwise the procedure tries to compute the real minimum distance |
---|
[fba86f8] | 1229 | to find out the error-correction capacity |
---|
| 1230 | @end format |
---|
[7f3ad4] | 1231 | RETURN: nothing; |
---|
| 1232 | EXAMPLE: example decodeRandom; shows an example |
---|
[288382] | 1233 | " |
---|
| 1234 | { |
---|
[fba86f8] | 1235 | intvec vopt = option(get); |
---|
| 1236 | |
---|
[288382] | 1237 | int i,j; |
---|
| 1238 | matrix h; |
---|
| 1239 | int dist, t, tim, tim2, tim3, timdist, timdec, timdist2, timdec2, timdec3; |
---|
| 1240 | ideal sys; |
---|
| 1241 | list tmp; |
---|
[fba86f8] | 1242 | int e; |
---|
| 1243 | if (size(#)>0) |
---|
| 1244 | { |
---|
| 1245 | e=#[1]; |
---|
| 1246 | } |
---|
[288382] | 1247 | |
---|
| 1248 | option(redSB); |
---|
| 1249 | def br=basering; |
---|
[7f3ad4] | 1250 | matrix h_full=genMDSMat(n,a); |
---|
[288382] | 1251 | matrix z[1][ncols(h_full)]; |
---|
[6e118cf] | 1252 | |
---|
| 1253 | //------------------ determine error-correction capacity ------------------- |
---|
[288382] | 1254 | for (i=1; i<=ncodes; i++) |
---|
| 1255 | { |
---|
| 1256 | setring br; |
---|
| 1257 | h=randomCheck(redun,n,1); |
---|
[bcb97dd] | 1258 | "check matrix:"; |
---|
[288382] | 1259 | print(h); |
---|
[fba86f8] | 1260 | if (e>0) |
---|
[288382] | 1261 | { |
---|
[fba86f8] | 1262 | t=e; |
---|
[288382] | 1263 | } else { |
---|
| 1264 | tim=rtimer; |
---|
| 1265 | tmp=mindist(h); |
---|
| 1266 | timdist=timdist+rtimer-tim; |
---|
| 1267 | timdist2=timdist2+tmp[2]; |
---|
| 1268 | dist=tmp[1]; |
---|
| 1269 | printf("d= %p",dist); |
---|
[7f3ad4] | 1270 | t=(dist-1) div 2; |
---|
| 1271 | } |
---|
| 1272 | tim2=rtimer; |
---|
[288382] | 1273 | tim3=rtimer; |
---|
[6e118cf] | 1274 | |
---|
| 1275 | //------------- generate the template system ---------------------- |
---|
[fba86f8] | 1276 | def A=sysQE(h,z,t); |
---|
[288382] | 1277 | setring A; |
---|
| 1278 | matrix word,y,rec; |
---|
| 1279 | ideal sys2,sys3; |
---|
| 1280 | matrix h=imap(br,h); |
---|
[7f3ad4] | 1281 | matrix g=dual_code(h); |
---|
[288382] | 1282 | ideal sys=qe; |
---|
| 1283 | print("The system is generated"); |
---|
[7f3ad4] | 1284 | timdec3=timdec3+rtimer-tim3; |
---|
[6e118cf] | 1285 | |
---|
[fba86f8] | 1286 | //------ modify the template according to every received word -------------- |
---|
[288382] | 1287 | for (j=1; j<=ntrials; j++) |
---|
| 1288 | { |
---|
[7f3ad4] | 1289 | word=randomvector(n-redun,1); |
---|
[288382] | 1290 | y=encode(transpose(word),g); |
---|
[7f3ad4] | 1291 | rec=errorRand(y,t,1); |
---|
| 1292 | sys2=add_synd(rec,h,redun,sys); |
---|
| 1293 | |
---|
[288382] | 1294 | tim=rtimer; |
---|
[7f3ad4] | 1295 | sys3=std(sys2); |
---|
| 1296 | timdec=timdec+rtimer-tim; |
---|
| 1297 | } |
---|
[c51e4a] | 1298 | timdec2=timdec2+rtimer-tim2; |
---|
| 1299 | kill A; |
---|
[fba86f8] | 1300 | option(set,vopt); |
---|
[7f3ad4] | 1301 | } |
---|
[288382] | 1302 | printf("Time for mindist: %p", timdist); |
---|
| 1303 | printf("Time for GB in mindist: %p", timdist); |
---|
| 1304 | printf("Time for decoding: %p", timdec2); |
---|
| 1305 | printf("Time for GB in decoding: %p", timdec); |
---|
[7f3ad4] | 1306 | printf("Time for sysQE in decoding: %p", timdec3); |
---|
| 1307 | } |
---|
| 1308 | example |
---|
[288382] | 1309 | { |
---|
| 1310 | "EXAMPLE:"; echo = 2; |
---|
| 1311 | int q=32; int n=25; int redun=n-11; int t=redun+1; |
---|
| 1312 | ring r=(q,a),x,dp; |
---|
[7f3ad4] | 1313 | |
---|
[288382] | 1314 | // correct 2 errors in 5 random binary codes, 50 trials each |
---|
[7f3ad4] | 1315 | decodeRandom(n,redun,5,50,2); |
---|
[288382] | 1316 | } |
---|
| 1317 | |
---|
[6e118cf] | 1318 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1319 | |
---|
[288382] | 1320 | |
---|
[fba86f8] | 1321 | proc decodeCode(matrix check, int ntrials, list #) |
---|
[7f3ad4] | 1322 | "USAGE: decodeCode(check, ntrials, [e]); check matrix, ntrials,e int |
---|
| 1323 | @format |
---|
| 1324 | - check is a check matrix for the code, |
---|
| 1325 | - ntrials is the number of received vectors per code to be |
---|
[fba86f8] | 1326 | corrected. |
---|
[7f3ad4] | 1327 | - If e is given it sets the correction capacity explicitly. It |
---|
[fba86f8] | 1328 | should be used in case one expects some lower bound, |
---|
[7f3ad4] | 1329 | otherwise the procedure tries to compute the real minimum distance |
---|
[fba86f8] | 1330 | to find out the error-correction capacity |
---|
| 1331 | @end format |
---|
[7f3ad4] | 1332 | RETURN: nothing; |
---|
| 1333 | EXAMPLE: example decodeCode; shows an example |
---|
[288382] | 1334 | " |
---|
| 1335 | { |
---|
[fba86f8] | 1336 | intvec vopt = option(get); |
---|
| 1337 | |
---|
[288382] | 1338 | int n=ncols(check); |
---|
| 1339 | int redun=nrows(check); |
---|
| 1340 | int i,j; |
---|
| 1341 | matrix h; |
---|
| 1342 | int dist, t, tim, tim2, tim3, timdist, timdec, timdist2, timdec2, timdec3; |
---|
| 1343 | ideal sys; |
---|
| 1344 | list tmp; |
---|
[fba86f8] | 1345 | int e; |
---|
| 1346 | if (size(#)>0) |
---|
| 1347 | { |
---|
| 1348 | e=#[1]; |
---|
| 1349 | } |
---|
[288382] | 1350 | |
---|
| 1351 | option(redSB); |
---|
| 1352 | def br=basering; |
---|
[7f3ad4] | 1353 | matrix h_full=genMDSMat(n,a); |
---|
[288382] | 1354 | matrix z[1][ncols(h_full)]; |
---|
| 1355 | setring br; |
---|
| 1356 | h=check; |
---|
[bcb97dd] | 1357 | "check matrix:"; |
---|
[288382] | 1358 | print(h); |
---|
[6e118cf] | 1359 | |
---|
| 1360 | //------------------ determine error-correction capacity ------------------- |
---|
[fba86f8] | 1361 | if (e>0) |
---|
[288382] | 1362 | { |
---|
[fba86f8] | 1363 | t=e; |
---|
[288382] | 1364 | } else { |
---|
| 1365 | tim=rtimer; |
---|
| 1366 | tmp=mindist(h); |
---|
| 1367 | timdist=timdist+rtimer-tim; |
---|
| 1368 | timdist2=timdist2+tmp[2]; |
---|
| 1369 | dist=tmp[1]; |
---|
| 1370 | printf("d= %p",dist); |
---|
[7f3ad4] | 1371 | t=(dist-1) div 2; |
---|
| 1372 | } |
---|
| 1373 | tim2=rtimer; |
---|
[288382] | 1374 | tim3=rtimer; |
---|
[6e118cf] | 1375 | |
---|
| 1376 | //------------- generate the template system ---------------------- |
---|
[fba86f8] | 1377 | def A=sysQE(h,z,t); |
---|
[288382] | 1378 | setring A; |
---|
| 1379 | matrix word,y,rec; |
---|
| 1380 | ideal sys2,sys3; |
---|
| 1381 | matrix h=imap(br,h); |
---|
[7f3ad4] | 1382 | matrix g=dual_code(h); |
---|
[288382] | 1383 | ideal sys=qe; |
---|
| 1384 | print("The system is generated"); |
---|
[7f3ad4] | 1385 | timdec3=timdec3+rtimer-tim3; |
---|
[6e118cf] | 1386 | |
---|
[7f3ad4] | 1387 | //--- modify the template according to every received word --------------- |
---|
[288382] | 1388 | for (j=1; j<=ntrials; j++) |
---|
| 1389 | { |
---|
[7f3ad4] | 1390 | word=randomvector(n-redun,1); |
---|
[288382] | 1391 | y=encode(transpose(word),g); |
---|
[7f3ad4] | 1392 | rec=errorRand(y,t,1); |
---|
| 1393 | sys2=add_synd(rec,h,redun,sys); |
---|
| 1394 | |
---|
[288382] | 1395 | tim=rtimer; |
---|
[7f3ad4] | 1396 | sys3=std(sys2); |
---|
| 1397 | timdec=timdec+rtimer-tim; |
---|
| 1398 | } |
---|
| 1399 | timdec2=timdec2+rtimer-tim2; |
---|
| 1400 | |
---|
[288382] | 1401 | printf("Time for mindist: %p", timdist); |
---|
| 1402 | printf("Time for GB in mindist: %p", timdist); |
---|
| 1403 | printf("Time for decoding: %p", timdec2); |
---|
| 1404 | printf("Time for GB in decoding: %p", timdec); |
---|
[7f3ad4] | 1405 | printf("Time for sysQE in decoding: %p", timdec3); |
---|
| 1406 | |
---|
[fba86f8] | 1407 | option(set,vopt); |
---|
[7f3ad4] | 1408 | } |
---|
| 1409 | example |
---|
[288382] | 1410 | { |
---|
| 1411 | "EXAMPLE:"; echo = 2; |
---|
| 1412 | int q=32; int n=25; int redun=n-11; int t=redun+1; |
---|
| 1413 | ring r=(q,a),x,dp; |
---|
| 1414 | matrix check=randomCheck(redun,n,1); |
---|
[7f3ad4] | 1415 | |
---|
| 1416 | // correct 2 errors in using the code above, 50 trials |
---|
| 1417 | decodeCode(check,50,2); |
---|
[288382] | 1418 | } |
---|
| 1419 | |
---|
| 1420 | |
---|
[6e118cf] | 1421 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1422 | // adding syndrome values to the template system |
---|
[288382] | 1423 | static proc add_synd (matrix rec, matrix check, int redun, ideal sys) |
---|
| 1424 | { |
---|
| 1425 | ideal result=sys; |
---|
| 1426 | matrix s[redun][1]=syndrome(check,rec); |
---|
| 1427 | for (int i=1; i<=redun; i++) |
---|
[7f3ad4] | 1428 | |
---|
[288382] | 1429 | { |
---|
[7f3ad4] | 1430 | result[i]=result[i]-s[i,1]; |
---|
[288382] | 1431 | } |
---|
| 1432 | return(result); |
---|
| 1433 | } |
---|
| 1434 | |
---|
[6e118cf] | 1435 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1436 | // evaluate a polynomial at a given point |
---|
[288382] | 1437 | static proc ev (poly f, matrix p) |
---|
| 1438 | { |
---|
| 1439 | if (ncols(p)>1) {ERROR("not a column vector");}; |
---|
| 1440 | int m=size(p); |
---|
| 1441 | poly temp=f; |
---|
| 1442 | for (int i=1; i<=m; i++) |
---|
| 1443 | { |
---|
| 1444 | temp=subst(temp,x(i),p[i,1]); |
---|
| 1445 | } |
---|
| 1446 | return(number(temp)); |
---|
| 1447 | } |
---|
| 1448 | |
---|
[fba86f8] | 1449 | /////////////////////////////////////////////////////////////////////////////// |
---|
[7f3ad4] | 1450 | // return index of an element in the ideal where it does not vanish at the |
---|
[fba86f8] | 1451 | //given point |
---|
[288382] | 1452 | static proc find_index (ideal G, matrix p) |
---|
| 1453 | { |
---|
| 1454 | if (ncols(p)>1) {ERROR("not a column vector");}; |
---|
| 1455 | int i=1; |
---|
| 1456 | int n=size(G); |
---|
| 1457 | while(i<=n) |
---|
| 1458 | { |
---|
| 1459 | if (ev(G[i],p)!=0) {return(i);} |
---|
| 1460 | i++; |
---|
| 1461 | } |
---|
| 1462 | return(-1); |
---|
| 1463 | } |
---|
| 1464 | |
---|
[6e118cf] | 1465 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1466 | // convert ideal to list |
---|
[288382] | 1467 | static proc ideal2list (ideal id) |
---|
| 1468 | { |
---|
| 1469 | list l; |
---|
| 1470 | for (int i=1; i<=size(id); i++) |
---|
| 1471 | { |
---|
| 1472 | l[i]=id[i]; |
---|
| 1473 | } |
---|
| 1474 | return(l); |
---|
| 1475 | } |
---|
| 1476 | |
---|
[6e118cf] | 1477 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1478 | // convert list to ideal |
---|
[288382] | 1479 | static proc list2ideal (list l) |
---|
| 1480 | { |
---|
| 1481 | ideal id; |
---|
| 1482 | for (int i=1; i<=size(l); i++) |
---|
| 1483 | { |
---|
| 1484 | id[i]=l[i]; |
---|
| 1485 | } |
---|
| 1486 | return(id); |
---|
| 1487 | } |
---|
| 1488 | |
---|
[fba86f8] | 1489 | /////////////////////////////////////////////////////////////////////////////// |
---|
[7f3ad4] | 1490 | // check whether given polynomial is divisible by some leading monomial of the |
---|
[fba86f8] | 1491 | //ideal |
---|
[288382] | 1492 | static proc divisible (poly m, ideal G) |
---|
| 1493 | { |
---|
| 1494 | for (int i=1; i<=size(G); i++) |
---|
| 1495 | { |
---|
[7f3ad4] | 1496 | if (m/leadmonom(G[i])!=0) {return(1);} |
---|
[288382] | 1497 | } |
---|
| 1498 | return(0); |
---|
| 1499 | } |
---|
| 1500 | |
---|
[6e118cf] | 1501 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1502 | |
---|
[288382] | 1503 | proc vanishId (list points) |
---|
[fba86f8] | 1504 | "USAGE: vanishId (points); point is a list of matrices |
---|
[7f3ad4] | 1505 | 'points' is a list of points for which the vanishing ideal is to be |
---|
| 1506 | constructed |
---|
[288382] | 1507 | RETURN: Vanishing ideal corresponding to the given set of points |
---|
[7f3ad4] | 1508 | EXAMPLE: example vanishId; shows an example |
---|
[288382] | 1509 | " |
---|
| 1510 | { |
---|
| 1511 | int m=size(points[1]); |
---|
| 1512 | int n=size(points); |
---|
[7f3ad4] | 1513 | |
---|
[288382] | 1514 | ideal G=1; |
---|
| 1515 | int i,k,j; |
---|
| 1516 | list temp; |
---|
| 1517 | poly h,cur; |
---|
[6e118cf] | 1518 | |
---|
| 1519 | //------------- proceed according to Farr-Gao algorithm ---------------- |
---|
[288382] | 1520 | for (k=1; k<=n; k++) |
---|
[7f3ad4] | 1521 | { |
---|
[288382] | 1522 | i=find_index(G,points[k]); |
---|
[7f3ad4] | 1523 | cur=G[i]; |
---|
[288382] | 1524 | for(j=i+1; j<=size(G); j++) |
---|
| 1525 | { |
---|
| 1526 | G[j]=G[j]-ev(G[j],points[k])/ev(G[i],points[k])*G[i]; |
---|
| 1527 | } |
---|
| 1528 | G=simplify(G,2); |
---|
| 1529 | temp=ideal2list(G); |
---|
| 1530 | temp=delete(temp,i); |
---|
[7f3ad4] | 1531 | G=list2ideal(temp); |
---|
[288382] | 1532 | for (j=1; j<=m; j++) |
---|
| 1533 | { |
---|
| 1534 | if (!divisible(x(j)*leadmonom(cur),G)) |
---|
| 1535 | { |
---|
| 1536 | attrib(G,"isSB",1); |
---|
[7f3ad4] | 1537 | h=NF((x(j)-points[k][j,1])*cur,G); |
---|
[288382] | 1538 | temp=ideal2list(G); |
---|
| 1539 | temp=insert(temp,h); |
---|
| 1540 | G=list2ideal(temp); |
---|
[7f3ad4] | 1541 | G=sort(G)[1]; |
---|
[288382] | 1542 | } |
---|
| 1543 | } |
---|
| 1544 | } |
---|
| 1545 | attrib(G,"isSB",1); |
---|
| 1546 | return(G); |
---|
[7f3ad4] | 1547 | } |
---|
| 1548 | example |
---|
[288382] | 1549 | { |
---|
| 1550 | "EXAMPLE:"; echo = 2; |
---|
| 1551 | ring r=3,(x(1..3)),dp; |
---|
[7f3ad4] | 1552 | |
---|
[288382] | 1553 | //generate all 3-vectors over GF(3) |
---|
[6e118cf] | 1554 | list points=pointsGen(3,1); |
---|
[7f3ad4] | 1555 | |
---|
[6e118cf] | 1556 | list points2=convPoints(points); |
---|
[7f3ad4] | 1557 | |
---|
[288382] | 1558 | //grasps the first 11 points |
---|
[6e118cf] | 1559 | list p=graspList(points2,1,11); |
---|
[288382] | 1560 | print(p); |
---|
[7f3ad4] | 1561 | |
---|
[288382] | 1562 | //construct the vanishing ideal |
---|
| 1563 | ideal id=vanishId(p); |
---|
| 1564 | print(id); |
---|
| 1565 | } |
---|
| 1566 | |
---|
[fba86f8] | 1567 | /////////////////////////////////////////////////////////////////////////////// |
---|
[7f3ad4] | 1568 | // construct the list of all vectors of length m with elements in p^e, where p |
---|
[fba86f8] | 1569 | //is characteristics |
---|
[6e118cf] | 1570 | proc pointsGen (int m, int e) |
---|
[288382] | 1571 | { |
---|
| 1572 | if (e>1) |
---|
| 1573 | { |
---|
| 1574 | list result; |
---|
| 1575 | int count=1; |
---|
| 1576 | int i,j; |
---|
| 1577 | list l=ringlist(basering); |
---|
| 1578 | int charac=l[1][1]; |
---|
| 1579 | number a=par(1); |
---|
| 1580 | list tmp; |
---|
| 1581 | for (i=1; i<=charac^(e*m); i++) |
---|
| 1582 | { |
---|
| 1583 | result[i]=tmp; |
---|
| 1584 | } |
---|
| 1585 | if (m==1) |
---|
| 1586 | { |
---|
| 1587 | result[count][m]=0; |
---|
| 1588 | count++; |
---|
| 1589 | for (j=1; j<=charac^(e)-1; j++) |
---|
[7f3ad4] | 1590 | { |
---|
[288382] | 1591 | result[count][m]=a^j; |
---|
| 1592 | count++; |
---|
| 1593 | } |
---|
| 1594 | return(result); |
---|
| 1595 | } |
---|
[7f3ad4] | 1596 | list prev=pointsGen(m-1,e); |
---|
[288382] | 1597 | for (i=1; i<=size(prev); i++) |
---|
| 1598 | { |
---|
| 1599 | result[count]=prev[i]; |
---|
| 1600 | result[count][m]=0; |
---|
| 1601 | count++; |
---|
| 1602 | for (j=1; j<=charac^(e)-1; j++) |
---|
| 1603 | { |
---|
| 1604 | result[count]=prev[i]; |
---|
| 1605 | result[count][m]=a^j; |
---|
| 1606 | count++; |
---|
| 1607 | } |
---|
| 1608 | } |
---|
| 1609 | return(result); |
---|
| 1610 | } |
---|
[7f3ad4] | 1611 | |
---|
[288382] | 1612 | if (e==1) |
---|
| 1613 | { |
---|
| 1614 | list result; |
---|
| 1615 | int count=1; |
---|
| 1616 | int i,j; |
---|
| 1617 | list l=ringlist(basering); |
---|
[7f3ad4] | 1618 | int charac=l[1][1]; |
---|
[288382] | 1619 | list tmp; |
---|
| 1620 | for (i=1; i<=charac^m; i++) |
---|
| 1621 | { |
---|
| 1622 | result[i]=tmp; |
---|
| 1623 | } |
---|
| 1624 | if (m==1) |
---|
| 1625 | { |
---|
| 1626 | for (j=0; j<=charac-1; j++) |
---|
[7f3ad4] | 1627 | { |
---|
[288382] | 1628 | result[count][m]=number(j); |
---|
| 1629 | count++; |
---|
| 1630 | } |
---|
| 1631 | return(result); |
---|
| 1632 | } |
---|
[7f3ad4] | 1633 | list prev=pointsGen(m-1,e); |
---|
[288382] | 1634 | for (i=1; i<=size(prev); i++) |
---|
| 1635 | { |
---|
| 1636 | for (j=0; j<=charac-1; j++) |
---|
| 1637 | { |
---|
| 1638 | result[count]=prev[i]; |
---|
| 1639 | result[count][m]=number(j); |
---|
| 1640 | count++; |
---|
| 1641 | } |
---|
| 1642 | } |
---|
| 1643 | return(result); |
---|
| 1644 | } |
---|
[7f3ad4] | 1645 | |
---|
[288382] | 1646 | } |
---|
| 1647 | |
---|
[6e118cf] | 1648 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1649 | // convert list to a column vector |
---|
[288382] | 1650 | static proc list2vec (list l) |
---|
| 1651 | { |
---|
| 1652 | matrix m[size(l)][1]; |
---|
| 1653 | for (int i=1; i<=size(l); i++) |
---|
| 1654 | { |
---|
| 1655 | m[i,1]=l[i]; |
---|
| 1656 | } |
---|
| 1657 | return(m); |
---|
| 1658 | } |
---|
| 1659 | |
---|
[6e118cf] | 1660 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1661 | // convert all the point in the list with list2vec |
---|
| 1662 | proc convPoints (list points) |
---|
[288382] | 1663 | { |
---|
| 1664 | for (int i=1; i<=size(points); i++) |
---|
| 1665 | { |
---|
| 1666 | points[i]=list2vec(points[i]); |
---|
| 1667 | } |
---|
| 1668 | return(points); |
---|
| 1669 | } |
---|
| 1670 | |
---|
[6e118cf] | 1671 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1672 | // extracts elements from l in the range m..n |
---|
| 1673 | proc graspList (list l, int m, int n) |
---|
[288382] | 1674 | { |
---|
| 1675 | list result; |
---|
| 1676 | int count=1; |
---|
| 1677 | for (int i=m; i<=n; i++) |
---|
| 1678 | { |
---|
| 1679 | result[count]=l[i]; |
---|
| 1680 | count++; |
---|
| 1681 | } |
---|
| 1682 | return(result); |
---|
| 1683 | } |
---|
| 1684 | |
---|
[6e118cf] | 1685 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1686 | // "characteristic" polynomial |
---|
[288382] | 1687 | static proc xi_gen (matrix p, int e, int s) |
---|
| 1688 | { |
---|
| 1689 | poly prod=1; |
---|
[7f3ad4] | 1690 | list rl=ringlist(basering); |
---|
[288382] | 1691 | int charac=rl[1][1]; |
---|
| 1692 | int l; |
---|
| 1693 | for (l=1; l<=s; l++) |
---|
| 1694 | { |
---|
| 1695 | prod=prod*(1-(x(l)-p[l,1])^(charac^e-1)); |
---|
| 1696 | } |
---|
| 1697 | return(prod); |
---|
| 1698 | } |
---|
| 1699 | |
---|
[6e118cf] | 1700 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1701 | // generating polynomials in Fitzgerald-Lax construction |
---|
[288382] | 1702 | static proc gener_funcs (matrix check, list points, int e, ideal id, int s) |
---|
| 1703 | { |
---|
| 1704 | int n=ncols(check); |
---|
| 1705 | if (n!=size(points)) {ERROR("Incompatible sizes of check and points");} |
---|
| 1706 | ideal xi; |
---|
| 1707 | int i,j; |
---|
| 1708 | for (i=1; i<=n; i++) |
---|
| 1709 | { |
---|
| 1710 | xi[i]=xi_gen(points[i],e,s); |
---|
| 1711 | } |
---|
| 1712 | ideal result; |
---|
| 1713 | int m=nrows(check); |
---|
| 1714 | poly sum; |
---|
| 1715 | for (i=1; i<=m; i++) |
---|
| 1716 | { |
---|
| 1717 | sum=0; |
---|
| 1718 | for (j=1; j<=n; j++) |
---|
| 1719 | { |
---|
| 1720 | sum=sum+check[i,j]*xi[j]; |
---|
| 1721 | } |
---|
| 1722 | result[i]=NF(sum,id); |
---|
| 1723 | } |
---|
| 1724 | return(result); |
---|
| 1725 | } |
---|
| 1726 | |
---|
[6e118cf] | 1727 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1728 | |
---|
[288382] | 1729 | proc sysFL (matrix check, matrix y, int t, int e, int s) |
---|
[fba86f8] | 1730 | "USAGE: sysFL (check,y,t,e,s); check,y matrix, t,e,s int |
---|
| 1731 | @format |
---|
[7f3ad4] | 1732 | - check is a check matrix of the code, |
---|
| 1733 | - y is a received word, |
---|
[fba86f8] | 1734 | - t the number of errors to correct, |
---|
[7f3ad4] | 1735 | - e is the extension degree, |
---|
[fba86f8] | 1736 | - s is the dimension of the point for the vanishing ideal |
---|
| 1737 | @end format |
---|
[3599e9] | 1738 | RETURN: the system of Fitzgerald-Lax for the given decoding problem |
---|
[7f3ad4] | 1739 | THEORY: Based on 'check' of the given linear code, the procedure constructs |
---|
[3599e9] | 1740 | the corresponding ideal constructed with a generalization of |
---|
| 1741 | Cooper's philosophy. For basics of the method @ref{Fitzgerald-Lax method}. |
---|
| 1742 | SEE ALSO: sysQE |
---|
[288382] | 1743 | EXAMPLE: example sysFL; shows an example |
---|
| 1744 | " |
---|
| 1745 | { |
---|
[7f3ad4] | 1746 | list rl=ringlist(basering); |
---|
| 1747 | int charac=rl[1][1]; |
---|
[288382] | 1748 | int n=ncols(check); |
---|
[7f3ad4] | 1749 | int m=nrows(check); |
---|
[6e118cf] | 1750 | list points=pointsGen(s,e); |
---|
| 1751 | list points2=convPoints(points); |
---|
[7f3ad4] | 1752 | list p=graspList(points2,1,n); |
---|
| 1753 | ideal id=vanishId(p,e); |
---|
| 1754 | ideal funcs=gener_funcs(check,p,e,id,s); |
---|
| 1755 | |
---|
[288382] | 1756 | ideal result; |
---|
| 1757 | poly temp; |
---|
| 1758 | int i,j,k; |
---|
[7f3ad4] | 1759 | |
---|
[6e118cf] | 1760 | //--------------- add vanishing realtions --------------------- |
---|
[288382] | 1761 | for (i=1; i<=t; i++) |
---|
[7f3ad4] | 1762 | { |
---|
[288382] | 1763 | for (j=1; j<=size(id); j++) |
---|
| 1764 | { |
---|
[7f3ad4] | 1765 | temp=id[j]; |
---|
[288382] | 1766 | for (k=1; k<=s; k++) |
---|
| 1767 | { |
---|
| 1768 | temp=subst(temp,x(k),x_var(i,k,s)); |
---|
[7f3ad4] | 1769 | } |
---|
[288382] | 1770 | result=result,temp; |
---|
| 1771 | } |
---|
[7f3ad4] | 1772 | } |
---|
| 1773 | |
---|
[6e118cf] | 1774 | //--------------- add field equations -------------------- |
---|
[288382] | 1775 | for (i=1; i<=t; i++) |
---|
| 1776 | { |
---|
| 1777 | for (k=1; k<=s; k++) |
---|
| 1778 | { |
---|
| 1779 | result=result,x_var(i,k,s)^(charac^e)-x_var(i,k,s); |
---|
| 1780 | } |
---|
| 1781 | } |
---|
| 1782 | for (i=1; i<=t; i++) |
---|
| 1783 | { |
---|
| 1784 | result=result,e(i)^(charac^e-1)-1; |
---|
| 1785 | } |
---|
[7f3ad4] | 1786 | |
---|
[288382] | 1787 | result=simplify(result,8); |
---|
[7f3ad4] | 1788 | |
---|
[6e118cf] | 1789 | //--------------- add check realtions -------------------- |
---|
[288382] | 1790 | poly sum; |
---|
[7f3ad4] | 1791 | matrix syn[m][1]=syndrome(check,y); |
---|
[288382] | 1792 | for (i=1; i<=size(funcs); i++) |
---|
[7f3ad4] | 1793 | { |
---|
[288382] | 1794 | sum=0; |
---|
| 1795 | for (j=1; j<=t; j++) |
---|
| 1796 | { |
---|
| 1797 | temp=funcs[i]; |
---|
| 1798 | for (k=1; k<=s; k++) |
---|
| 1799 | { |
---|
| 1800 | temp=subst(temp,x(k),x_var(j,k,s)); |
---|
| 1801 | } |
---|
[7f3ad4] | 1802 | sum=sum+temp*e(j); |
---|
[288382] | 1803 | } |
---|
| 1804 | result=result,sum-syn[i,1]; |
---|
| 1805 | } |
---|
[7f3ad4] | 1806 | |
---|
[288382] | 1807 | result=simplify(result,2); |
---|
[7f3ad4] | 1808 | |
---|
[288382] | 1809 | points=points2; |
---|
[7f3ad4] | 1810 | export points; |
---|
[288382] | 1811 | return(result); |
---|
[7f3ad4] | 1812 | } |
---|
[fba86f8] | 1813 | example |
---|
[288382] | 1814 | { |
---|
| 1815 | "EXAMPLE:"; echo = 2; |
---|
[fba86f8] | 1816 | intvec vopt = option(get); |
---|
[7f3ad4] | 1817 | |
---|
[288382] | 1818 | list l=FLpreprocess(3,1,11,2,""); |
---|
| 1819 | def r=l[1]; |
---|
| 1820 | setring r; |
---|
[7f3ad4] | 1821 | int s_work=l[2]; |
---|
| 1822 | |
---|
[288382] | 1823 | //the check matrix of [11,6,5] ternary code |
---|
| 1824 | matrix h[5][11]=1,0,0,0,0,1,1,1,-1,-1,0, |
---|
| 1825 | 0,1,0,0,0,1,1,-1,1,0,-1, |
---|
| 1826 | 0,0,1,0,0,1,-1,1,0,1,-1, |
---|
| 1827 | 0,0,0,1,0,1,-1,0,1,-1,1, |
---|
| 1828 | 0,0,0,0,1,1,0,-1,-1,1,1; |
---|
| 1829 | matrix g=dual_code(h); |
---|
| 1830 | matrix x[1][6]; |
---|
| 1831 | matrix y[1][11]=encode(x,g); |
---|
| 1832 | //disturb with 2 errors |
---|
[00142a] | 1833 | matrix rec[1][11]=errorInsert(y,list(2,4),list(1,-1)); |
---|
[288382] | 1834 | |
---|
[7f3ad4] | 1835 | //the Fitzgerald-Lax system |
---|
[288382] | 1836 | ideal sys=sysFL(h,rec,2,1,s_work); |
---|
| 1837 | print(sys); |
---|
| 1838 | option(redSB); |
---|
| 1839 | ideal red_sys=std(sys); |
---|
[7f3ad4] | 1840 | red_sys; |
---|
[fba86f8] | 1841 | // read the solutions from this redGB |
---|
[288382] | 1842 | // the points are (0,0,1) and (0,1,0) with error values 1 and -1 resp. |
---|
| 1843 | // use list points to find error positions; |
---|
[fba86f8] | 1844 | points; |
---|
[7f3ad4] | 1845 | |
---|
[fba86f8] | 1846 | option(set,vopt); |
---|
[288382] | 1847 | } |
---|
| 1848 | |
---|
[6e118cf] | 1849 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1850 | // preprocessing steps for the Fitzgerald-Lax scheme |
---|
[288382] | 1851 | proc FLpreprocess (int p, int e, int n, int t, string minp) |
---|
| 1852 | { |
---|
| 1853 | ring r1=p,x,dp; |
---|
| 1854 | int s=1; |
---|
| 1855 | while(p^(s*e)<n) |
---|
| 1856 | { |
---|
| 1857 | s++; |
---|
[7f3ad4] | 1858 | } |
---|
[288382] | 1859 | list var_ord; |
---|
| 1860 | int i,j; |
---|
| 1861 | int count=1; |
---|
| 1862 | for (i=s; i>=1; i--) |
---|
| 1863 | { |
---|
| 1864 | var_ord[count]=string("x("+string(i)+")"); |
---|
| 1865 | count++; |
---|
| 1866 | } |
---|
| 1867 | for (i=t; i>=1; i--) |
---|
| 1868 | { |
---|
| 1869 | var_ord[count]=string("e("+string(i)+")"); |
---|
| 1870 | count++; |
---|
| 1871 | for (j=s; j>=1; j--) |
---|
| 1872 | { |
---|
| 1873 | var_ord[count]=string("x1("+string(s*(i-1)+j)+")"); |
---|
| 1874 | count++; |
---|
| 1875 | } |
---|
[7f3ad4] | 1876 | } |
---|
| 1877 | |
---|
[288382] | 1878 | list rl; |
---|
| 1879 | list tmp; |
---|
[7f3ad4] | 1880 | |
---|
[288382] | 1881 | if (e>1) |
---|
| 1882 | { |
---|
| 1883 | rl[1]=tmp; |
---|
| 1884 | rl[1][1]=p; |
---|
| 1885 | rl[1][2]=tmp; |
---|
| 1886 | rl[1][2][1]=string("a"); |
---|
| 1887 | rl[1][3]=tmp; |
---|
[7f3ad4] | 1888 | rl[1][3][1]=tmp; |
---|
[288382] | 1889 | rl[1][3][1][1]=string("lp"); |
---|
| 1890 | rl[1][3][1][2]=1; |
---|
| 1891 | rl[1][4]=ideal(0); |
---|
| 1892 | } else { |
---|
| 1893 | rl[1]=p; |
---|
| 1894 | } |
---|
[7f3ad4] | 1895 | |
---|
[288382] | 1896 | rl[2]=var_ord; |
---|
[7f3ad4] | 1897 | |
---|
[288382] | 1898 | rl[3]=tmp; |
---|
[7f3ad4] | 1899 | rl[3][1]=tmp; |
---|
[288382] | 1900 | rl[3][1][1]=string("lp"); |
---|
| 1901 | intvec v=1; |
---|
| 1902 | for (i=1; i<=size(var_ord)-1; i++) |
---|
| 1903 | { |
---|
| 1904 | v=v,1; |
---|
| 1905 | } |
---|
[7f3ad4] | 1906 | rl[3][1][2]=v; |
---|
[288382] | 1907 | rl[3][2]=tmp; |
---|
[7f3ad4] | 1908 | rl[3][2][1]=string("C"); |
---|
[288382] | 1909 | rl[3][2][2]=intvec(0); |
---|
[7f3ad4] | 1910 | |
---|
| 1911 | rl[4]=ideal(0); |
---|
| 1912 | |
---|
[288382] | 1913 | def r2=ring(rl); |
---|
| 1914 | setring r2; |
---|
[7f3ad4] | 1915 | list l=ringlist(r2); |
---|
[288382] | 1916 | if (e>1) |
---|
| 1917 | { |
---|
[7f3ad4] | 1918 | execute(string("poly f="+minp)); |
---|
| 1919 | ideal id=f; |
---|
[288382] | 1920 | l[1][4]=id; |
---|
| 1921 | } |
---|
[7f3ad4] | 1922 | |
---|
[288382] | 1923 | def r=ring(l); |
---|
[7f3ad4] | 1924 | setring r; |
---|
| 1925 | |
---|
[288382] | 1926 | return(list(r,s)); |
---|
| 1927 | } |
---|
| 1928 | |
---|
[6e118cf] | 1929 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1930 | // imitating two indeces |
---|
[288382] | 1931 | static proc x_var (int i, int j, int s) |
---|
[7f3ad4] | 1932 | { |
---|
[288382] | 1933 | return(x1(s*(i-1)+j)); |
---|
| 1934 | } |
---|
| 1935 | |
---|
[6e118cf] | 1936 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1937 | // random vector of length n with entries from p^e, p the characteristic |
---|
[fba86f8] | 1938 | static proc randomvector(int n, int e) |
---|
[7f3ad4] | 1939 | { |
---|
[288382] | 1940 | int i; |
---|
| 1941 | matrix result[n][1]; |
---|
| 1942 | for (i=1; i<=n; i++) |
---|
| 1943 | { |
---|
[fba86f8] | 1944 | result[i,1]=asElement(random_prime_vector(e)); |
---|
[288382] | 1945 | } |
---|
[7f3ad4] | 1946 | return(result); |
---|
[288382] | 1947 | } |
---|
| 1948 | |
---|
[fba86f8] | 1949 | /////////////////////////////////////////////////////////////////////////////// |
---|
[7f3ad4] | 1950 | // "convert" representation of an element from the field extension from vector |
---|
[fba86f8] | 1951 | //to an elelemnt |
---|
[288382] | 1952 | static proc asElement(list l) |
---|
| 1953 | { |
---|
| 1954 | number s; |
---|
| 1955 | int i; |
---|
| 1956 | number w=1; |
---|
[7f3ad4] | 1957 | if (size(l)>1) {w=par(1);} |
---|
[288382] | 1958 | for (i=0; i<=size(l)-1; i++) |
---|
| 1959 | { |
---|
| 1960 | s=s+w^i*l[i+1]; |
---|
| 1961 | } |
---|
| 1962 | return(s); |
---|
| 1963 | } |
---|
| 1964 | |
---|
[6e118cf] | 1965 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1966 | // random vector of length n with entries from p, p the characteristic |
---|
[fba86f8] | 1967 | static proc random_prime_vector (int n) |
---|
[7f3ad4] | 1968 | { |
---|
[fba86f8] | 1969 | list rl=ringlist(basering); |
---|
| 1970 | int i, charac; |
---|
| 1971 | for (i=2; i<=rl[1][1]; i++) |
---|
[288382] | 1972 | { |
---|
[fba86f8] | 1973 | if (rl[1][1] mod i ==0) |
---|
[7f3ad4] | 1974 | { |
---|
[fba86f8] | 1975 | break; |
---|
| 1976 | } |
---|
[288382] | 1977 | } |
---|
[7f3ad4] | 1978 | charac=i; |
---|
| 1979 | |
---|
[288382] | 1980 | list l; |
---|
[7f3ad4] | 1981 | |
---|
[288382] | 1982 | for (i=1; i<=n; i++) |
---|
| 1983 | { |
---|
| 1984 | l=l+list(random(0,charac-1)); |
---|
| 1985 | } |
---|
| 1986 | return(l); |
---|
| 1987 | } |
---|
| 1988 | |
---|
[6e118cf] | 1989 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1990 | |
---|
[fba86f8] | 1991 | proc decodeRandomFL(int n, int redun, int p, int e, int t, int ncodes, int ntrials, string minpol) |
---|
[7f3ad4] | 1992 | "USAGE: decodeRandomFL(redun,p,e,n,t,ncodes,ntrials,minpol); |
---|
[fba86f8] | 1993 | @format |
---|
[7f3ad4] | 1994 | - n is length of codes generated, redun = redundancy of codes |
---|
| 1995 | generated, |
---|
| 1996 | - p is characteristics, |
---|
| 1997 | - e is the extension degree, |
---|
| 1998 | - t is the number of errors to correct, |
---|
[fba86f8] | 1999 | - ncodes is the number of random codes to be processed, |
---|
| 2000 | - ntrials is the number of received vectors per code to be corrected, |
---|
[7f3ad4] | 2001 | - minpol: due to some pecularities of SINGULAR one needs to provide |
---|
[fba86f8] | 2002 | minimal polynomial for the extension explicitly |
---|
| 2003 | @end format |
---|
[288382] | 2004 | RETURN: nothing |
---|
[fba86f8] | 2005 | EXAMPLE: example decodeRandomFL; shows an example |
---|
| 2006 | " |
---|
[288382] | 2007 | { |
---|
[fba86f8] | 2008 | intvec vopt = option(get); |
---|
| 2009 | |
---|
[7f3ad4] | 2010 | list l=FLpreprocess(p,e,n,t,minpol); |
---|
| 2011 | |
---|
[288382] | 2012 | def r=l[1]; |
---|
[7f3ad4] | 2013 | int s_work=l[2]; |
---|
[288382] | 2014 | export(s_work); |
---|
| 2015 | setring r; |
---|
[7f3ad4] | 2016 | |
---|
[288382] | 2017 | int i,j; |
---|
| 2018 | matrix h, g, word, y, rec; |
---|
| 2019 | int dist, tim, tim2, tim3, timdist, timdec, timdist2, timdec2, timdec3; |
---|
| 2020 | ideal sys, sys2, sys3; |
---|
[7f3ad4] | 2021 | list tmp; |
---|
| 2022 | |
---|
| 2023 | option(redSB); |
---|
| 2024 | matrix z[1][n]; |
---|
[288382] | 2025 | |
---|
| 2026 | for (i=1; i<=ncodes; i++) |
---|
| 2027 | { |
---|
[7f3ad4] | 2028 | h=randomCheck(redun,n,e); |
---|
[288382] | 2029 | g=dual_code(h); |
---|
[7f3ad4] | 2030 | tim2=rtimer; |
---|
[6e118cf] | 2031 | tim3=rtimer; |
---|
| 2032 | |
---|
[7f3ad4] | 2033 | //---------------- generate the template system ----------------------- |
---|
| 2034 | sys=sysFL(h,z,t,e,s_work); |
---|
[288382] | 2035 | timdec3=timdec3+rtimer-tim3; |
---|
[7f3ad4] | 2036 | |
---|
[fba86f8] | 2037 | //------ modifying the template according to the received word --------- |
---|
[288382] | 2038 | for (j=1; j<=ntrials; j++) |
---|
| 2039 | { |
---|
| 2040 | word=randomvector(n-redun,1); |
---|
| 2041 | y=encode(transpose(word),g); |
---|
| 2042 | rec=errorRand(y,t,e); |
---|
[7f3ad4] | 2043 | sys2=LF_add_synd(rec,h,sys); |
---|
[288382] | 2044 | tim=rtimer; |
---|
| 2045 | sys3=std(sys2); |
---|
| 2046 | timdec=timdec+rtimer-tim; |
---|
[7f3ad4] | 2047 | } |
---|
[288382] | 2048 | timdec2=timdec2+rtimer-tim2; |
---|
[7f3ad4] | 2049 | } |
---|
| 2050 | |
---|
[288382] | 2051 | printf("Time for decoding: %p", timdec2); |
---|
| 2052 | printf("Time for GB in decoding: %p", timdec); |
---|
[7f3ad4] | 2053 | printf("Time for generating Fitzgerald-Lax system during decoding: %p", timdec3); |
---|
| 2054 | |
---|
[fba86f8] | 2055 | option(set,vopt); |
---|
[7f3ad4] | 2056 | } |
---|
[fba86f8] | 2057 | example |
---|
[288382] | 2058 | { |
---|
| 2059 | "EXAMPLE:"; echo = 2; |
---|
[7f3ad4] | 2060 | |
---|
| 2061 | // correcting one error for one random binary code of length 25, |
---|
[fba86f8] | 2062 | //redundancy 14; 300 words are processed |
---|
| 2063 | decodeRandomFL(25,14,2,1,1,1,300,""); |
---|
[288382] | 2064 | } |
---|
| 2065 | |
---|
[6e118cf] | 2066 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 2067 | // add syndrome values to the template system in FL |
---|
[288382] | 2068 | static proc LF_add_synd (matrix rec, matrix check, ideal sys) |
---|
| 2069 | { |
---|
| 2070 | int redun=nrows(check); |
---|
| 2071 | ideal result=sys; |
---|
| 2072 | matrix s[redun][1]=syndrome(check,rec); |
---|
| 2073 | for (int i=size(sys)-redun+1; i<=size(sys); i++) |
---|
| 2074 | { |
---|
| 2075 | result[i]=result[i]-s[i-size(sys)+redun,1]; |
---|
| 2076 | } |
---|
| 2077 | return(result); |
---|
| 2078 | } |
---|
| 2079 | |
---|
| 2080 | |
---|
| 2081 | /* |
---|
[00142a] | 2082 | ////////////// SOME RELATIVELY EASY EXAMPLES ////////////// |
---|
| 2083 | /////////////////// THAT RUN AROUND ONE MINUTE //////////////// |
---|
| 2084 | |
---|
| 2085 | "EXAMPLE:"; echo = 2; |
---|
[7f3ad4] | 2086 | int q=128; int n=120; int redun=n-30; |
---|
[00142a] | 2087 | ring r=(q,a),x,dp; |
---|
[fba86f8] | 2088 | decodeRandom(n,redun,1,1,6); |
---|
[00142a] | 2089 | |
---|
[7f3ad4] | 2090 | int q=128; int n=120; int redun=n-20; |
---|
[00142a] | 2091 | ring r=(q,a),x,dp; |
---|
[fba86f8] | 2092 | decodeRandom(n,redun,1,1,9); |
---|
[00142a] | 2093 | |
---|
| 2094 | int q=128; int n=120; int redun=n-10; |
---|
| 2095 | ring r=(q,a),x,dp; |
---|
[fba86f8] | 2096 | decodeRandom(n,redun,1,1,19); |
---|
[00142a] | 2097 | |
---|
| 2098 | int q=256; int n=150; int redun=n-10; |
---|
| 2099 | ring r=(q,a),x,dp; |
---|
[fba86f8] | 2100 | decodeRandom(n,redun,1,1,22); |
---|
[00142a] | 2101 | |
---|
[288382] | 2102 | ////////////// SOME HARD EXAMPLES ////////////////////// |
---|
| 2103 | ////// THAT MAYBE WILL BE DOABLE LATER /////////////// |
---|
| 2104 | |
---|
| 2105 | 1.) These random instances are not doable in <=1000 sec. |
---|
| 2106 | |
---|
| 2107 | "EXAMPLE:"; echo = 2; |
---|
[7f3ad4] | 2108 | int q=128; int n=120; int redun=n-40; |
---|
[288382] | 2109 | ring r=(q,a),x,dp; |
---|
[fba86f8] | 2110 | decodeRandom(n,redun,1,1,6); |
---|
[288382] | 2111 | |
---|
| 2112 | redun=n-30; |
---|
[fba86f8] | 2113 | decodeRandom(n,redun,1,1,8); |
---|
[288382] | 2114 | |
---|
| 2115 | redun=n-20; |
---|
[fba86f8] | 2116 | decodeRandom(n,redun,1,1,12); |
---|
[288382] | 2117 | |
---|
| 2118 | redun=n-10; |
---|
[fba86f8] | 2119 | decodeRandom(n,redun,1,1,24); |
---|
[288382] | 2120 | |
---|
[7f3ad4] | 2121 | int q=256; int n=150; int redun=n-10; |
---|
[288382] | 2122 | ring r=(q,a),x,dp; |
---|
[fba86f8] | 2123 | decodeRandom(n,redun,1,1,26); |
---|
[288382] | 2124 | |
---|
| 2125 | |
---|
| 2126 | 2.) Generic decoding is hard! |
---|
| 2127 | |
---|
| 2128 | int q=32; int n=31; int redun=n-16; int t=3; |
---|
| 2129 | ring r=(q,a),(V(1..n),U(n..1),s(redun..1)),(dp(n),lp(n),dp(redun)); |
---|
| 2130 | matrix check[redun][n]= 1,1,0,1,1,0,0,0,1,0,1,0,0,1,0,0,1,0,0,0,0, |
---|
| 2131 | 0,0,0,0,0,0,0,0,0,0,0,1,1,0,1,1,0,0,0,1, |
---|
| 2132 | 0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0, |
---|
| 2133 | 0,0,1,1,0,1,1,0,0,0,1,0,1,0,0,1,0,0,1,0,0, |
---|
| 2134 | 0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,1,1,0,0, |
---|
| 2135 | 0,1,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0, |
---|
| 2136 | 0,0,0,0,1,1,0,1,1,0,0,0,1,0,1,0,0,1,0,0,1, |
---|
| 2137 | 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,1,1, |
---|
| 2138 | 0,0,0,1,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0, |
---|
| 2139 | 0,0,0,0,0,0,1,1,0,1,1,0,0,0,1,0,1,0,0,1,0, |
---|
| 2140 | 0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0, |
---|
| 2141 | 1,1,0,0,0,1,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0, |
---|
| 2142 | 0,0,0,0,0,0,0,0,1,1,0,1,1,0,0,0,1,0,1,0,0, |
---|
| 2143 | 1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1, |
---|
| 2144 | 1,0,1,1,0,0,0,1,0,1,0,0,1,0,0,1,0,0,0,0,0, |
---|
| 2145 | 0,0,0,0,0,0,0,0,0,0,1,1,0,1,1,0,0,0,1,0,1, |
---|
| 2146 | 0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0, |
---|
| 2147 | 0,1,1,0,1,1,0,0,0,1,0,1,0,0,1,0,0,1,0,0,0, |
---|
| 2148 | 0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,1,1,0,0,0,1, |
---|
| 2149 | 0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0, |
---|
| 2150 | 0,0,0,1,1,0,1,1,0,0,0,1,0,1,0,0,1,0,0,1,0, |
---|
| 2151 | 0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,1,1,0,0, |
---|
| 2152 | 0,1,0,1,0,0,1,0,0,1; |
---|
| 2153 | matrix rec[1][n]; |
---|
| 2154 | |
---|
| 2155 | def A=sysQE(check,rec,t,1,2); |
---|
| 2156 | setring A; |
---|
| 2157 | print(qe); |
---|
[7f3ad4] | 2158 | ideal red_qe=stdfglm(qe); |
---|
[288382] | 2159 | |
---|
[7f3ad4] | 2160 | */ |
---|