1############################################################################# 2## 3#A codegen.gi GUAVA library Reinald Baart 4#A &Jasper Cramwinckel 5#A &Erik Roijackers 6## &David Joyner 7## 8## This file contains info/functions for generating codes 9## 10## BCHCode modified 9-2004 11## added 11-2004: 12## comment in GoppaCode 13## GeneralizedReedSolomonCode and record fields 14## points, degree, ring, 15## (added SpecialDecoder field later: SetSpecialDecoder(C, GRSDecoder);) 16## EvaluationCode and record fields 17## points, basis, ring, 18## OnePointAGCode and record fields 19## points, curve, multiplicity, ring 20## weighted option for GeneralizedReedSolomonCode and record fields 21## points, degree, ring, weights 22## minor bug in EvaluationCode fixed 11-26-2004 23## GeneratorMatCodeNC added 12-18-2004 (after release of guava 2.0) 24## added 5-10-2005: SetSpecialDecoder(C, CyclicDecoder); to 25## the cyclic codes which are not BCH codes 26## in codegen.gi, line 2066, replace C.GeneratorMat 27## by C.EvaluationMat 28## added 5-13-2005: QQRCode 29## added 11-27-2005: CheckMatCodeMutable with *mutable* 30## generator and check matrices 31## added 1-10-2006: QQRCodeNC Greg Coy and David Joyner 32## added 1-2006: FerreroDesignCode, written with Peter Mayr 33## added 12-2007 (CJ): ExtendedReedSolomonCode, QuasiCyclicCode, 34## CyclicMDSCode, FourNegacirculantSelfDualCode functions. 35## 36 37DeclareRepresentation( "IsCodeDefaultRep", 38 IsAttributeStoringRep and IsComponentObjectRep and IsCode, 39 ["name", "history", "lowerBoundMinimumDistance", 40 "upperBoundMinimumDistance", "boundsCoveringRadius"]); 41 42DeclareHandlingByNiceBasis( "IsLinearCodeRep", "for linear codes" ); 43 44InstallTrueMethod( IsLinearCode, IsLinearCodeRep ); 45 46#T The following is of course a hack, as is much of the 47#T ``pretend as if you were a vector space'' stuff. 48#T (`IsLinearCode' changes harmless codes into vector spaces; 49#T `AsLinearCode' would be clean, but now it is too late ...) 50InstallTrueMethod( IsFreeLeftModule, IsLinearCodeRep ); 51 52### functions to work with NiceBasis functionality for Linear Codes. 53InstallHandlingByNiceBasis( "IsLinearCodeRep", rec( 54 detect:= function( R, gens, V, zero ) 55 return IsCodewordCollection( V ); 56 end, 57 58 NiceFreeLeftModuleInfo:= ReturnFalse, 59 60 NiceVector:= function( C, w ) 61 return VectorCodeword( w ); 62 end, 63 64 UglyVector:= function( C, v ) 65 return Codeword( v, Length( v ), LeftActingDomain( C ) ); 66 end ) ); 67 68 69############################################################################# 70## 71#F ElementsCode( <L> [, <name> ], <F> ) . . . . . . code from list of words 72## 73 74InstallMethod(ElementsCode, "list,name,Field", true, 75 [IsList,IsString,IsField], 0, 76function(L, name, F) 77 local test, C; 78 L := Codeword(L, F); 79 test := WordLength(L[1]); 80 if ForAny(L, i->WordLength(i) <> test) then 81 Error("All elements must have the same length"); 82 fi; 83 test := FamilyObj(L[1]); 84 if ForAny(L, i->FamilyObj(i) <> test) then 85 Error ("All elements must have the same family"); 86 fi; 87 L := Set(L); 88 C := Objectify(NewType(CollectionsFamily(test), IsCodeDefaultRep), rec()); 89 SetLeftActingDomain(C, F); 90 SetAsSSortedList(C, L); 91 C!.name := name; 92 return C; 93end); 94 95InstallOtherMethod(ElementsCode, "list,Field", true, [IsList, IsField], 0, 96function(L,F) 97 return ElementsCode(L, "user defined unrestricted code", F); 98end); 99 100InstallOtherMethod(ElementsCode, "list,name,fieldsize", true, 101 [IsList, IsString, IsInt], 0, 102function(L, name, q) 103 return ElementsCode(L, name, GF(q)); 104end); 105 106InstallOtherMethod(ElementsCode, "list,size of field", true, [IsList,IsInt], 0, 107function(L, q) 108 return ElementsCode(L, "user defined unrestricted code", GF(q)); 109end); 110 111##Create a linear code from a list of Codeword generators 112LinearCodeByGenerators := function(F, gens) 113 114 local V; 115 V:= Objectify( NewType( FamilyObj( gens ), 116 IsLeftModule and 117 IsLinearCodeRep and IsCodeDefaultRep ), 118 rec() ); 119 SetLeftActingDomain( V, F ); 120 SetGeneratorsOfLeftModule( V, AsList( One(F)*gens ) ); 121 SetIsLinearCode(V, true); 122 SetWordLength(V, Length(gens[1])); 123 return V; 124 125end; 126 127 128############################################################################# 129## 130#F RandomCode( <n>, <M> [, <F>] ) . . . . . . . . random unrestricted code 131## 132InstallMethod(RandomCode, "wordlength,codesize,field", true, 133 [IsInt, IsInt, IsField], 0, 134function(n, M, F) 135 local L; 136 if Size(F)^n < M then 137 Error(Size(F),"^",n," < ",M); 138 fi; 139 L := []; 140 while Length(L) < M do 141 AddSet(L, List([1..n], i -> Random(F))); 142 od; 143 return ElementsCode(L, "random unrestricted code", F); 144end); 145 146InstallOtherMethod(RandomCode, "wordlength,codesize", true, [IsInt, IsInt], 0, 147function(n, M) 148 return RandomCode(n, M, GF(2)); 149end); 150 151InstallOtherMethod(RandomCode, "wordlength,codesize,fieldsize", true, 152 [IsInt, IsInt, IsInt], 0, 153function(n, M, q) 154 return RandomCode(n, M, GF(q)); 155end); 156 157 158############################################################################# 159## 160#F HadamardCode( <H | n> [, <t>] ) . Hadamard code of <t>'th kind, order <n> 161## 162 163InstallMethod(HadamardCode, "matrix, kind (int)", true, [IsMatrix, IsInt], 0, 164function(H, t) 165 local n, vec, C; 166 n := Length(H); 167 if H * TransposedMat(H) <> n * IdentityMat(n,n) then 168 Error("The matrix is not a Hadamard matrix"); 169 fi; 170 H := (H-1)/(-2); 171 if t = 1 then 172 H := TransposedMat(TransposedMat(H){[2..n]}); 173 C := ElementsCode(H, Concatenation("Hadamard code of order ", 174 String(n)), GF(2) ); 175 C!.lowerBoundMinimumDistance := n/2; 176 C!.upperBoundMinimumDistance := n/2; 177 vec := NullVector(n); 178 # this was ... := NullVector(n+1); 179 # but this seems to be wrong -- Eric Minkes 180 vec[1] := 1; 181 vec[n/2+1] := Size(C) - 1; 182 SetInnerDistribution(C, vec); 183 elif t = 2 then 184 H := ShallowCopy( TransposedMat(TransposedMat(H){[2..n]}) ); 185 Append(H, 1 - H); 186 C := ElementsCode(H, Concatenation("Hadamard code of order ", 187 String(n)), GF(2) ); 188 C!.lowerBoundMinimumDistance := n/2 - 1; 189 C!.upperBoundMinimumDistance := n/2 - 1; 190 vec := NullVector(n); 191 vec[1] := 1; 192 vec[n] := 1; 193 # this was ... [n+1]... 194 # but this seems to be wrong -- Eric Minkes 195 vec[n/2] := Size(C)/2 - 1; 196 vec[n/2+1] := Size(C)/2 - 1; 197 SetInnerDistribution(C, vec); 198 else 199 Append(H, 1 - H); 200 C := ElementsCode( H, Concatenation("Hadamard code of order ", 201 String(n)), GF(2) ); 202 C!.lowerBoundMinimumDistance := n/2; 203 C!.upperBoundMinimumDistance := n/2; 204 vec := NullVector(n); 205 vec[1] := 1; 206 vec[n+1] := 1; 207 vec[n/2+1] := Size(C) - 2; 208 SetInnerDistribution(C, vec); 209 fi; 210 return(C); 211end); 212 213InstallOtherMethod(HadamardCode, "order, kind (int)", true, [IsInt, IsInt], 0, 214function(n, t) 215 return HadamardCode(HadamardMat(n), t); 216end); 217 218InstallOtherMethod(HadamardCode, "matrix", true, [IsMatrix], 0, 219function(H) 220 return HadamardCode(H, 3); 221end); 222 223InstallOtherMethod(HadamardCode, "order", true, [IsInt], 0, 224function(n) 225 return HadamardCode(HadamardMat(n), 3); 226end); 227 228InstallOtherMethod(HadamardCode, "matrix, kind (string)", true, 229 [IsMatrix, IsString], 0, 230function(H, t) 231 if t in ["a", "A", "1"] then 232 return HadamardCode(H, 1); 233 elif t in ["b", "B", "2"] then 234 return HadamardCode(H, 2); 235 else 236 return HadamardCode(H, 3); 237 fi; 238end); 239 240InstallOtherMethod(HadamardCode, "order, kind (string)", true, 241 [IsInt, IsString], 0, 242function(n, t) 243 if t in ["a", "A", "1"] then 244 return HadamardCode(HadamardMat(n), 1); 245 elif t in ["b", "B", "2"] then 246 return HadamardCode(HadamardMat(n), 2); 247 else 248 return HadamardCode(HadamardMat(n), 3); 249 fi; 250end); 251 252############################################################################# 253## 254#F ConferenceCode( <n | M> ) . . . . . . . . . . code from conference matrix 255## 256 257InstallMethod(ConferenceCode, "matrix", true, [IsMatrix], 0, 258function(S) 259 local n, I, J, F, els, w, wd, C; 260 n := Length(S); 261 if S*TransposedMat(S) <> (n-1)*IdentityMat(n) or 262 TransposedMat(S) <> S then 263 Error("argument must be a symmetric conference matrix"); 264 fi; 265 # Normalize S by multiplying rows and columns: 266 for I in [2..n] do 267 if S[I][1] <> 1 then 268 for J in [1..n] do 269 S[I][J] := S[I][J] * -1; 270 od; 271 fi; 272 od; 273 for J in [2..n] do 274 if S[1][J] <> 1 then 275 for I in [1..n] do 276 S[I][J] := S[I][J] * -1; 277 od; 278 fi; 279 od; 280 # Strip first row and first column: 281 S := List([2..n], i-> S[i]{[2..n]}); 282 n := n - 1; 283 284 F := GF(2); 285 els := [NullWord(n, F)]; 286 I := IdentityMat(n); 287 J := NullMat(n,n) + 1; 288 Append(els, Codeword(1/2 * (S+I+J), F)); 289 Append(els, Codeword(1/2 * (-S+I+J), F)); 290 Add( els, Codeword( List([1..n], x -> One(F) ), n, One(F) ) ); 291 C := ElementsCode( els, "conference code", F); 292 w := Weight(AsSSortedList(C)[2]); 293 wd := List([1..n+1], x -> 0); 294 wd[1] := 1; wd[n+1] := 1; 295 wd[w+1] := Size(C) - 2; 296 SetWeightDistribution(C, wd); 297 C!.lowerBoundMinimumDistance := (n-1)/2; 298 C!.upperBoundMinimumDistance := (n-1)/2; 299 return C; 300end); 301 302InstallOtherMethod(ConferenceCode, "integer", true, [IsInt], 0, 303function(n) 304 local LegendreSym, zero, QRes, E, I, S, F, els, J, w, wd, C; 305 306 LegendreSym := function(i) 307 if i = zero then 308 return 0; 309 elif i in QRes then 310 return 1; 311 else 312 return -1; 313 fi; 314 end; 315 316 if (not IsPrimePowerInt(n)) or (n mod 4 <> 1) then 317 Error ("n must be a primepower and n mod 4 = 1"); 318 fi; 319 E := AsSSortedList(GF(n)); 320 zero := E[1]; 321 QRes := []; 322 for I in E do 323 AddSet(QRes, I^2); 324 od; 325 S := List(E, i->List(E, j->LegendreSym(j-i))); 326 327 328 F := GF(2); 329 els := [NullWord(n, F)]; 330 I := IdentityMat(n); 331 J := NullMat(n,n) + 1; 332 Append(els, Codeword(1/2 * (S+I+J), F)); 333 Append(els, Codeword(1/2 * (-S+I+J), F)); 334 Add( els, Codeword( List([1..n], x -> One(F) ), n, One(F) ) ); 335 C := ElementsCode( els, "conference code", F); 336 w := Weight(AsSSortedList(C)[2]); 337 wd := List([1..n+1], x -> 0); 338 wd[1] := 1; wd[n+1] := 1; 339 wd[w+1] := Size(C) - 2; 340 SetWeightDistribution(C, wd); 341 C!.lowerBoundMinimumDistance := (n-1)/2; 342 C!.upperBoundMinimumDistance := (n-1)/2; 343 return C; 344end); 345 346 347############################################################################# 348## 349#F MOLSCode( [ <n>, ] <q> ) . . . . . . . . . . . . . . . . code from MOLS 350## 351## MOLSCode([n, ] q) returns a (n, q^2, n-1) code over GF(q) 352## by creating n-2 mutually orthogonal latin squares of size q. 353## If n is omitted, a wordlength of 4 will be set. 354## If there are no n-2 MOLS known, the code will return an error 355## 356 357InstallMethod(MOLSCode, "wordlength,size of field", true, [IsInt, IsInt], 0, 358function(n,q) 359 local M, els, i, j, k, wd, C; 360 M := MOLS(q, n-2); 361 if M = false then 362 Error("No ", n-2, " MOLS of order ", q, " are known"); 363 else 364 els := []; 365 for i in [1..q] do 366 for j in [1..q] do 367 els[(i-1)*q+j] := []; 368 els[(i-1)*q+j][1] := i-1; 369 els[(i-1)*q+j][2]:= j-1; 370 for k in [3..n] do 371 els[(i-1)*q+j][k]:= M[k-2][i][j]; 372 od; 373 od; 374 od; 375 C := ElementsCode( els, Concatenation("code generated by ", 376 String(n-2), " MOLS of order ", String(q)), GF(q) ); 377 C!.lowerBoundMinimumDistance := n - 1; 378 C!.upperBoundMinimumDistance := n - 1; 379 wd := List( [1..n+1], x -> 0 ); 380 wd[1] := 1; 381 wd[n] := (q-1) * n; 382 wd[n+1] := (q-1) * (q + 1 - n); 383 SetWeightDistribution(C, wd); 384 return C; 385 fi; 386end); 387 388InstallOtherMethod(MOLSCode, "wordlength,field", true, [IsInt, IsField], 0, 389function(n, F) 390 return MOLSCode(n, Size(F)); 391end); 392 393InstallOtherMethod(MOLSCode, "fieldsize", true, [IsInt], 0, 394function(q) 395 return MOLSCode(4, q); 396end); 397 398InstallOtherMethod(MOLSCode, "field", true, [IsField], 0, 399function(F) 400 return MOLSCode(4, Size(F)); 401end); 402 403 404## Was commented out in gap 3.4.4 Guava version. 405 406############################################################################# 407## 408#F QuadraticCosetCode( <Q> ) . . . . . . . . . . coset of RM(1,m) in R(2,m) 409## 410## QuadraticCosetCode(Q) returns a coset of the ReedMullerCode of 411## order 1 (R(1,m)) in R(2,m) where m is the size of square matrix Q. 412## Q is the upper triangular matrix that defines the quadratic part of 413## the boolean functions that are in the coset. 414## 415#QuadraticCosetCode := function(arg) 416# local Q, m, V, R, RM1, k, f, C, h, wd; 417# if Length(arg) = 1 and IsMat(arg[1]) then 418# Q := arg[1]; 419# else 420# Error("usage: QuadraticCosetCode( <mat> )"); 421# fi; 422# m := Length(Q); 423# V := Tuples(Elements(GF(2)), m); 424# R := V*Q*TransposedMat(V); 425# f := List([1..2^m], i->R[i][i]); 426# RM1 := Concatenation(NullMat(1,2^m,GF(2))+GF(2).one, 427# TransposedMat(Tuples(Elements(GF(2)), m))); 428# k := Length(RM1); 429# C := rec( 430# isDomain := true, 431# isCode := true, 432# operations := CodeOps, 433# baseField := GF(2), 434# wordLength := 2^m, 435# elements := Codeword(List(Tuples(Elements(GF(2)), k) * RM1, t-> t+f)), 436# lowerBoundMinimumDistance := 2^(m-1), 437# upperBoundMinimumDistance := 2^(m-1) 438# ); 439# h := RankMat(Q + TransposedMat(Q))/2; 440# wd := NullMat(1, 2^m+1)[1]; 441# wd[2^(m-1) - 2^(m-h-1) + 1] := 2^(2*h); 442# wd[2^(m-1) + 1] := 2^(m+1) - 2^(2*h + 1); 443# wd[2^(m-1) + 2^(m-h-1) + 1] := 2^(2*h); 444# C.weightDistribution := wd; 445# C.name := "quadratic coset code"; 446# return C; 447#end; 448 449 450############################################################################# 451## 452#F GeneratorMatCode( <G> [, <name> ], <F> ) . . code from generator matrix 453## 454 455InstallMethod(GeneratorMatCode, "generator matrix, name, field", true, 456 [IsMatrix, IsString, IsField], 0, 457function(G, name, F) 458 local C; 459 if (Length(G) = 0) or (Length(G[1]) = 0) then 460 Error("Use NullCode to generate a code with dimension 0"); 461 fi; 462 G := G*One(F); 463 G := BaseMat(G); 464 C := LinearCodeByGenerators(F, Codeword(G,F)); 465 SetGeneratorMat(C, G); 466 SetWordLength(C, Length(G[1])); 467 C!.name := name; 468 return C; 469end); 470 471InstallOtherMethod(GeneratorMatCode, "generator matrix, field", true, 472 [IsMatrix, IsField], 0, 473function(G, F) 474 return GeneratorMatCode(G, "code defined by generator matrix", F); 475end); 476 477InstallOtherMethod(GeneratorMatCode, "generator matrix, name, size of field", 478 true, [IsMatrix, IsString, IsInt], 0, 479function(G, name, q) 480 return GeneratorMatCode(G, name, GF(q)); 481end); 482 483InstallOtherMethod(GeneratorMatCode, "generator matrix, size of field", true, 484 [IsMatrix, IsInt], 0, 485function(G, q) 486 return GeneratorMatCode(G, "code defined by generator matrix", GF(q)); 487end); 488 489############################################################################# 490## 491#F GeneratorMatCodeNC( <G> , <F> ) . . code from generator matrix 492## 493## same as GeneratorMatCode but does not compute upper + lower bounds 494## for the minimum distance or covering radius 495## 496 497InstallMethod(GeneratorMatCodeNC, "generator matrix, field", true, 498 [IsMatrix, IsField], 0, 499function(G, F) 500 return GeneratorMatCode(G, "code defined by generator matrix, NC", F); 501end); 502 503 504 505############################################################################# 506## 507#F CheckMatCodeMutable( <H> [, <name> ], <F> ) . . code from check matrix 508## The generator matrix is mutable 509## 510 511InstallMethod(CheckMatCodeMutable, "check matrix, name, field", true, 512 [IsMatrix, IsString, IsField], 0, 513function(H, name, F) 514 local G, H2, C, dimC; 515 if (Length(H) = 0) or (Length(H[1]) = 0) then 516 Error("use WholeSpaceCode to generate a code with redundancy 0"); 517 fi; 518 H := VectorCodeword(Codeword(H, F)); 519 H2 := BaseMat(H); 520 if Length(H2) < Length(H) then 521 H := H2; 522 fi; 523 524 # Get generator matrix from H 525 if IsInStandardForm(H, false) then 526 dimC := Length(H[1]) - Length(H); 527 if dimC = 0 then 528 G := []; 529 else 530 G := TransposedMat(Concatenation(IdentityMat( dimC, F ), 531 List(-H, x->x{[1..dimC]}))); 532 fi; 533 else 534 G := NullspaceMat(TransposedMat(H)); 535 fi; 536 if Length(G) = 0 then # Call NullCode. Length(H=I) = n-k, in this case. 537 # Length(G) = k = 0 => n = Length(H). 538 C := NullCode(Length(H), F); 539 C!.name := name; 540 return C; 541 fi; 542 543 C := LinearCodeByGenerators(F,Codeword(G, F)); 544 SetGeneratorMat(C, ShallowCopy(G)); 545 SetCheckMat(C, ShallowCopy(H)); 546 SetWordLength(C, Length(G[1])); 547 C!.name := name; 548 return C; 549end); 550 551InstallOtherMethod(CheckMatCodeMutable, "check matrix, field", 552true, [IsMatrix, IsField], 0, 553function(H, F) 554 return CheckMatCode(H, "code defined by check matrix", F); 555end); 556 557InstallOtherMethod(CheckMatCodeMutable, "check matrix, name, size of field", true, [IsMatrix, IsString, IsInt], 0, 558function(H, name, q) 559 return CheckMatCode(H, name, GF(q)); 560end); 561 562InstallOtherMethod(CheckMatCodeMutable, "check matrix, size of field", 563true, [IsMatrix, IsInt], 0, 564function(H, q) 565 return CheckMatCodeMutable(H, "code defined by check matrix", GF(q)); 566end); 567 568 569 570############################################################################# 571## 572#F CheckMatCode( <H> [, <name> ], <F> ) . . . . . . code from check matrix 573## 574 575InstallMethod(CheckMatCode, "check matrix, name, field", true, 576 [IsMatrix, IsString, IsField], 0, 577function(H, name, F) 578 local G, H2, C, dimC; 579 if (Length(H) = 0) or (Length(H[1]) = 0) then 580 Error("use WholeSpaceCode to generate a code with redundancy 0"); 581 fi; 582 H := VectorCodeword(Codeword(H, F)); 583 H2 := BaseMat(H); 584 if Length(H2) < Length(H) then 585 H := H2; 586 fi; 587 588 # Get generator matrix from H 589 if IsInStandardForm(H, false) then 590 dimC := Length(H[1]) - Length(H); 591 if dimC = 0 then 592 G := []; 593 else 594 G := TransposedMat(Concatenation(IdentityMat( dimC, F ), 595 List(-H, x->x{[1..dimC]}))); 596 fi; 597 else 598 G := NullspaceMat(TransposedMat(H)); 599 fi; 600 if Length(G) = 0 then # Call NullCode. Length(H=I) = n-k, in this case. 601 # Length(G) = k = 0 => n = Length(H). 602 C := NullCode(Length(H), F); 603 C!.name := name; 604 return C; 605 fi; 606 607 C := LinearCodeByGenerators(F,Codeword(G, F)); 608 SetGeneratorMat(C, G); 609 SetCheckMat(C, H); 610 SetWordLength(C, Length(G[1])); 611 C!.name := name; 612 return C; 613end); 614 615InstallOtherMethod(CheckMatCode, "check matrix, field", true, 616 [IsMatrix, IsField], 0, 617function(H, F) 618 return CheckMatCode(H, "code defined by check matrix", F); 619end); 620 621InstallOtherMethod(CheckMatCode, "check matrix, name, size of field", true, 622 [IsMatrix, IsString, IsInt], 0, 623function(H, name, q) 624 return CheckMatCode(H, name, GF(q)); 625end); 626 627InstallOtherMethod(CheckMatCode, "check matrix, size of field", true, 628 [IsMatrix, IsInt], 0, 629function(H, q) 630 return CheckMatCode(H, "code defined by check matrix", GF(q)); 631end); 632 633 634############################################################################# 635## 636#F RandomLinearCode( <n>, <k> [, <F>] ) . . . . . . . . random linear code 637## 638 639InstallMethod(RandomLinearCode, "n,k,field", true, [IsInt, IsInt, IsField], 0, 640function(n, k, F) 641 if k = 0 then 642 return NullCode( n, F ); 643 else 644 return GeneratorMatCode(PermutedCols(List(IdentityMat(k,F), i -> 645 Concatenation(i,List([k+1..n],j->Random(F)))), 646 Random(SymmetricGroup(n))), "random linear code", F); 647 fi; 648end); 649 650InstallOtherMethod(RandomLinearCode, "n,k,size of field", true, 651 [IsInt, IsInt, IsInt], 0, 652function(n, k, q) 653 return RandomLinearCode(n, k, GF(q)); 654end); 655 656InstallOtherMethod(RandomLinearCode, "n,k", true, [IsInt, IsInt], 0, 657function(n, k) 658 return RandomLinearCode(n, k, GF(2)); 659end); 660 661 662############################################################################# 663## 664#F HammingCode( <r> [, <F>] ) . . . . . . . . . . . . . . . . Hamming code 665## 666 667InstallMethod(HammingCode, "r,F", true, [IsInt, IsField], 0, 668function(r, F) 669 local q, H, H2, C, i, j, n, TupAllow, wd; 670 671 TupAllow := function(W) 672 local l; 673 l := 1; 674 while (W[l] = Zero(F)) and l < Length(W) do 675 l := l + 1; 676 od; 677 return (W[l] = One(F)); 678 end; 679 680 q := Size(F); 681 if not IsPrimePowerInt(q) then 682 Error("q must be prime power"); 683 fi; 684 H := Tuples(AsSSortedList(F), r); 685 H2 := []; 686 j := 1; 687 for i in [1..Length(H)] do 688 if TupAllow(H[i]) then 689 H2[j] := H[i]; 690 j := j + 1; 691 fi; 692 od; 693 n := (q^r-1)/(q-1); 694 C := CheckMatCode(TransposedMat(H2), Concatenation("Hamming (", String(r), 695 ",", String(q), ") code"), F); 696 C!.lowerBoundMinimumDistance := 3; 697 C!.upperBoundMinimumDistance := 3; 698 C!.boundsCoveringRadius := [ 1 ]; 699 SetIsPerfectCode(C, true); 700 SetIsSelfDualCode(C, false); 701 SetSpecialDecoder(C, HammingDecoder); 702 if q = 2 then 703 SetIsNormalCode(C, true); 704 wd := [1, 0]; 705 for i in [2..n] do 706 Add(wd, 1/i * (Binomial(n, i-1) - wd[i] - (n-i+2)*wd[i-1])); 707 od; 708 SetWeightDistribution(C, wd); 709 fi; 710 return C; 711end); 712 713InstallOtherMethod(HammingCode, "r,size of field", true, [IsInt, IsInt], 0, 714function(r, q) 715 return HammingCode(r, GF(q)); 716end); 717 718InstallOtherMethod(HammingCode, "r", true, [IsInt], 0, 719function(r) 720 return HammingCode(r, GF(2)); 721end); 722 723 724############################################################################# 725## 726#F SimplexCode( <r>, <F> ) . The SimplexCode is the Dual of the HammingCode 727## 728 729InstallMethod(SimplexCode, "redundancy, field", true, [IsInt, IsField], 0, 730function(r, F) 731 local C; 732 C := DualCode( HammingCode(r,F) ); 733 C!.name := "simplex code"; 734 if F = GF(2) then 735 SetIsNormalCode(C, true); 736 C!.boundsCoveringRadius := [ 2^(r-1) - 1 ]; 737 fi; 738 return C; 739end); 740 741InstallOtherMethod(SimplexCode, "redundancy, fieldsize", true,[IsInt,IsInt],0, 742function(r, q) 743 return SimplexCode(r, GF(q)); 744end); 745 746InstallOtherMethod(SimplexCode, "redundancy", true, [IsInt], 0, 747function(r) 748 return SimplexCode(r, GF(2)); 749end); 750 751 752############################################################################# 753## 754#F ReedMullerCode( <r>, <k> ) . . . . . . . . . . . . . . Reed-Muller code 755## 756## ReedMullerCode(r, k) creates a binary Reed-Muller code of dimension k, 757## order r; 0 <= r <= k 758## 759 760InstallMethod(ReedMullerCode, "dimension, order", true, [IsInt,IsInt], 0, 761function(r,k) 762 local mat,c,src,dest,index,num,dim,C,wd, h,t,A, bcr; 763 764 if r > k then 765 return ReedMullerCode(k, r); #for compatibility with older versions 766 #of guava, It used to be 767 #ReedMullerCode(k, r); 768 fi; 769 mat := [ [] ]; 770 num := 2^k; 771 dim := Sum(List([0..r], x->Binomial(k,x))); 772 for t in [1..num] do 773 mat[1][t] := Z(2)^0; 774 od; 775 if r > 0 then 776 Append(mat, TransposedMat(Tuples ([0*Z(2), Z(2)^0], k))); 777 for t in [2..r] do 778 for index in Combinations([1..k], t) do 779 dest := List([1..2^k], i->Product(index, j->mat[j+1][i])); 780 Append(mat, [dest]); 781 od; 782 od; 783 fi; 784 C := GeneratorMatCode( mat, Concatenation("Reed-Muller (", String(r), ",", 785 String(k), ") code"), GF(2) ); 786 C!.lowerBoundMinimumDistance := 2^(k-r); 787 C!.upperBoundMinimumDistance := 2^(k-r); 788 SetIsPerfectCode(C, false); 789 SetIsSelfDualCode(C, (2*r = k-1)); 790 if r = 0 then 791 wd := List([1..num + 1], x -> 0); 792 wd[1] := 1; 793 wd[num+1] := 1; 794 SetWeightDistribution(C, wd); 795 elif r = 1 then 796 wd := List([1..num + 1], x -> 0); 797 wd[1] := 1; 798 wd[num + 1] := 1; 799 wd[num / 2 + 1] := Size(C) - 2; 800 SetWeightDistribution(C, wd); 801 elif r = 2 then 802 wd := List([1..num + 1], x -> 0); 803 wd[1] := 1; 804 wd[num + 1] := 1; 805 for h in [1..QuoInt(k,2)] do 806 A := 2^(h*(h+1)); 807 for t in [0..2*h-1] do 808 A := A*(2^(k-t)-1); 809 od; 810 for t in [1..h] do 811 A := A/(2^(2*t)-1); 812 od; 813 wd[2^(k-1)+2^(k-1-h)+1] := A; 814 wd[2^(k-1)-2^(k-1-h)+1] := A; 815 od; 816 wd[2^(k-1)+1] := Size(C)-Sum(wd); 817 SetWeightDistribution(C, wd); 818 fi; 819 820 bcr := BoundsCoveringRadius( C ); 821 822 if 0 <= r and r <= k-3 then 823 if IsEvenInt( r ) then 824 C!.boundsCoveringRadius := 825 [ Maximum( 2^(k-r-3)*(r+4), bcr[1] ) 826 .. Maximum( bcr ) ]; 827 else 828 C!.boundsCoveringRadius := 829 [ Maximum( 2^(k-r-3)*(r+5), bcr[ 1 ] ) 830 .. Maximum( bcr ) ]; 831 fi; 832 fi; 833 834 if r = k then 835 SetCoveringRadius(C, 0); 836 C!.boundsCoveringRadius := [ 0 ]; 837 elif r = k - 1 then 838 SetCoveringRadius(C, 1); 839 C!.boundsCoveringRadius := [ 1 ]; 840 elif r = k - 2 then 841 SetCoveringRadius(C, 2); 842 C!.boundsCoveringRadius := [ 2 ]; 843 elif r = k - 3 then 844 if IsEvenInt( k ) then 845 SetCoveringRadius(C, k + 2 ); 846 C!.boundsCoveringRadius := [ k + 2 ]; 847 else 848 SetCoveringRadius(C, k + 1 ); 849 C!.boundsCoveringRadius := [ k + 1 ]; 850 fi; 851 elif r = 0 then 852 SetCoveringRadius(C, 2^(k-1) ); 853 C!.boundsCoveringRadius := [ 2^(k-1) ]; 854 elif r = 1 then 855 if IsEvenInt( k ) then 856 SetCoveringRadius(C, 2^(k-1) - 2^(k/2-1) ); 857 C!.boundsCoveringRadius := [ 2^(k-1) - 2^(k/2-1) ]; 858 elif k = 5 then 859 SetCoveringRadius(C, 12 ); 860 C!.boundsCoveringRadius := [ 12 ]; 861 elif k = 7 then 862 SetCoveringRadius(C, 56 ); 863 C!.boundsCoveringRadius := [ 56 ]; 864 elif k >= 15 then 865 C!.boundsCoveringRadius := [ 2^(k-1) - 2^((k+1)/2)*27/64 866 .. 2^(k-1) - 2^( QuoInt(k,2)-1 ) ]; 867 else 868 C!.boundsCoveringRadius := [ 2^(k-1) - 2^((k+1)/2)/2 869 .. 2^(k-1) - 2^( QuoInt(k,2)-1 ) ]; 870 fi; 871 elif r = 2 then 872 if k = 6 then 873 SetCoveringRadius(C, 18 ); 874 C!.boundsCoveringRadius := [ 18 ]; 875 elif k = 7 then 876 C!.boundsCoveringRadius := [ 40 .. 44 ]; 877 elif k = 8 then 878 C!.boundsCoveringRadius := [ 84 879 .. Maximum( bcr ) ]; 880 fi; 881 elif r = 3 then 882 if k = 7 then 883 C!.boundsCoveringRadius := [ 20 .. 23 ]; 884 fi; 885 elif r = 4 then 886 if k = 8 then 887 C!.boundsCoveringRadius := [ 22 888 .. Maximum( bcr ) ]; 889 fi; 890 fi; 891 892 if r = 1 and 893 ( IsEvenInt( k ) or k = 3 or k = 5 or k = 7 ) then 894 SetIsNormalCode(C, true); 895 fi; 896 897 return C; 898end); 899 900 901############################################################################# 902## 903#F LexiCode( <M | n>, <d>, <F> ) . . . . . Greedy code with standard basis 904## 905 906InstallMethod(LexiCode, "basis,distance,field", true, 907 [IsMatrix, IsInt, IsField], 0, 908function(M, d, F) 909 local base, n, k, one, zero, elms, Sz, vec, word, i, dist, carry, pos, C; 910 base := VectorCodeword(Codeword(M, F)); 911 n := Length(base[1]); 912 k := Length(base); 913 one := One(F); 914 zero := Zero(F); 915 elms := []; 916 Sz := 0; 917 vec := NullVector(k,F); 918 repeat 919 word := vec * base; 920 i := 1; 921 dist := d; 922 while (dist>=d) and (i <= Sz) do 923 dist := DistanceVecFFE(word, elms[i]); 924 i := i + 1; 925 od; 926 if dist >= d then 927 Add(elms,ShallowCopy(word)); 928 Sz := Sz + 1; 929 fi; 930 # generate the (lexicographical) next word in F^k 931 carry := true; 932 pos := k; 933 while carry and (pos > 0) do 934 if vec[pos] = zero then 935 carry := false; 936 vec[pos] := one; 937 else 938 vec[pos] := PrimitiveRoot(F)^(LogFFE(vec[pos],PrimitiveRoot(F))+1); 939 if vec[pos] = one then 940 vec[pos] := zero; 941 else 942 carry := false; 943 fi; 944 fi; 945 pos := pos - 1; 946 od; 947 until carry; 948 if Size(F) = 2 then # or even (2^(2^LogInt(LogInt(q,2),2)) = q) ? 949 C := GeneratorMatCode(elms, "lexicode", F); 950 else 951 C := ElementsCode(elms, "lexicode", F); 952 fi; 953 C!.lowerBoundMinimumDistance := d; 954 return C; 955end); 956 957InstallOtherMethod(LexiCode, "basis,distance,size of field", true, 958 [IsMatrix, IsInt, IsInt], 0, 959function(M, d, q) 960 return LexiCode(M, d, GF(q)); 961end); 962 963InstallOtherMethod(LexiCode, "wordlength,distance,field", true, 964 [IsInt, IsInt, IsField], 0, 965function(n, d, F) 966 return LexiCode(IdentityMat(n,F), d, F); 967end); 968 969InstallOtherMethod(LexiCode, "wordlength,distance,size of field", true, 970 [IsInt, IsInt, IsInt], 0, 971function(n, d, q) 972 return LexiCode(IdentityMat(n, GF(q)), d, GF(q)); 973end); 974 975 976############################################################################# 977## 978#F GreedyCode( <M>, <d> [, <F>] ) . . . . Greedy code from list of elements 979## 980 981InstallMethod(GreedyCode, "matrix,design distance,Field", true, 982 [IsMatrix, IsInt, IsField], 0, 983function(M,d,F) 984 local space, n, word, elms, Sz, i, dist, C; 985 space := VectorCodeword(Codeword(M,F)); 986 n := Length(space[1]); 987 elms := [space[1]]; 988 Sz := 1; 989 for word in space do 990 i := 1; 991 repeat 992 dist := DistanceVecFFE(word, elms[i]); 993 i := i + 1; 994 until dist < d or i > Sz; 995 if dist >= d then 996 Add(elms, word); 997 Sz := Sz + 1; 998 fi; 999 od; 1000 C := ElementsCode(elms, "Greedy code, user defined basis", F); 1001 C!.lowerBoundMinimumDistance := d; 1002 return C; 1003end); 1004 1005InstallOtherMethod(GreedyCode, "matrix,design distance,size of field", true, 1006 [IsMatrix, IsInt, IsInt], 0, 1007function(M,d,q) 1008 return GreedyCode(M,d,GF(q)); 1009end); 1010 1011InstallOtherMethod(GreedyCode, "matrix,design distance", true, 1012 [IsMatrix,IsInt], 0, 1013function(M,d) 1014 return GreedyCode(M,d,DefaultField(Flat(M))); 1015end); 1016 1017 1018############################################################################# 1019## 1020#F AlternantCode( <r>, <Y> [, <alpha>], <F> ) . . . . . . . Alternant code 1021## 1022 1023InstallMethod(AlternantCode, "redundancy, Y, alpha, field", true, 1024 [IsInt, IsList, IsList, IsField], 0, 1025function(r, Y, els, F) 1026 local C, n, q, i, temp; 1027 n := Length(Y); 1028 els := Set(VectorCodeword(Codeword(els, F))); 1029 Y := VectorCodeword(Codeword(Y, F) ); 1030 1031 if ForAny(Y, i-> i = Zero(F)) then 1032 Error("Y contains zero"); 1033 elif Length(els) <> Length(Y) then 1034 Error("<Y> and <alpha> have inequal length or <alpha> is not distinct"); 1035 fi; 1036 q := Characteristic(F); 1037 temp := NullMat(n, n, F); 1038 for i in [1..n] do 1039 temp[i][i] := Y[i]; 1040 od; 1041 Y := temp; 1042 C := CheckMatCode( BaseMat(VerticalConversionFieldMat( List([0..r-1], 1043 i -> List([1..n], j-> els[j]^i)) * Y)), "alternant code", F ); 1044 C!.lowerBoundMinimumDistance := r + 1; 1045 return C; 1046end); 1047 1048InstallOtherMethod(AlternantCode, "redundancy, Y, alpha, fieldsize", true, 1049 [IsInt, IsList, IsList, IsInt], 0, 1050function(r, Y, a, q) 1051 return AlternantCode(r, Y, a, GF(q)); 1052end); 1053 1054InstallOtherMethod(AlternantCode, "redundancy, Y, field", true, 1055 [IsInt, IsList, IsField], 0, 1056function(r, Y, F) 1057 return AlternantCode(r, Y, AsSSortedList(F){[2..Length(Y)+1]}, F); 1058end); 1059 1060InstallOtherMethod(AlternantCode, "redundancy, Y, fieldsize", true, 1061 [IsInt, IsList, IsInt], 0, 1062function(r, Y, q) 1063 return AlternantCode(r, Y, AsSSortedList(GF(q)){[2..Length(Y)+1]}, GF(q)); 1064end); 1065 1066 1067############################################################################# 1068## 1069#F GoppaCode( <G>, <L | n> ) . . . . . . . . . . . . . . classical Goppa code 1070## 1071 1072InstallGlobalFunction(GoppaCode, 1073function(arg) 1074 local C, GP, F, L, n, q, m, r, zero, temp; 1075 1076 GP := PolyCodeword(Codeword(arg[1])); 1077 F := CoefficientsRing(DefaultRing(GP)); 1078 q := Characteristic(F); 1079 m := Dimension(F); 1080 F := GF(q); 1081 zero := Zero(F); 1082 r := DegreeOfLaurentPolynomial(GP); 1083 1084 # find m 1085 if IsInt(arg[2]) then 1086 n := arg[2]; 1087 m := Maximum(m, LogInt(n,q)); 1088 repeat 1089 L := Filtered(AsSSortedList(GF(q^m)),i -> Value(GP,i) <> zero); 1090 m := m + 1; 1091 until Length(L) >= n; 1092 m := m - 1; 1093 L := L{[1..n]}; 1094 else 1095 L := arg[2]; 1096 n := Length(L); 1097 m := Maximum(m, Dimension(DefaultField(L))); 1098 fi; 1099 C := CheckMatCode( BaseMat(VerticalConversionFieldMat( List([0..r-1], 1100 i-> List(L, j-> (j)^i / Value(GP, j) )) )), "classical Goppa code", F); 1101 1102 # Make the code 1103 temp := Factors(GP); 1104 if (q = 2) and (Length(temp) = Length(Set(temp))) then 1105 # second condition checks if the roots of G are distinct 1106 C!.lowerBoundMinimumDistance := Minimum(n, 2*r + 1); 1107 else 1108 C!.lowerBoundMinimumDistance := Minimum(n, r + 1); 1109 fi; 1110 return C; 1111end); 1112 1113 1114############################################################################# 1115## 1116#F CordaroWagnerCode( <n> ) . . . . . . . . . . . . . . Cordaro-Wagner code 1117## 1118 1119InstallMethod(CordaroWagnerCode, "length", true, [IsInt], 0, 1120function(n) 1121 local r, C, zero, one, F, d, wd; 1122 if n < 2 then 1123 Error("n must be 2 or more"); 1124 fi; 1125 r := Int((n+1)/3); 1126 d := (2 * r - Int( (n mod 3) / 2) ); 1127 F := GF(2); 1128 zero := Zero(F); 1129 one := One(F); 1130 C := GeneratorMatCode( [Concatenation(List([1..r],i -> zero), 1131 List([r+1..n],i -> one)), 1132 Concatenation(List([r+1..n],i -> one), List([1..r], 1133 i -> zero))], "Cordaro-Wagner code", F ); 1134 C!.lowerBoundMinimumDistance := d; 1135 C!.upperBoundMinimumDistance := d; 1136 wd := List([1..n+1], i-> 0); 1137 wd[1] := 1; 1138 wd[2*r+1] := 1; 1139 wd[n-r+1] := wd[n-r+1] + 2; 1140 SetWeightDistribution(C, wd); 1141 return C; 1142end); 1143 1144 1145############################################################################# 1146## 1147#F GeneralizedSrivastavaCode( <a>, <w>, <z> [, <t>] [, <F>] ) . . . . . . 1148## 1149 1150InstallMethod(GeneralizedSrivastavaCode, "a,w,z,t,F", true, 1151 [IsList, IsList, IsList, IsInt, IsField], 0, 1152function(a,w,z,t,F) 1153 local C, n, s, i, H; 1154 a := VectorCodeword(Codeword(a, F)); 1155 w := VectorCodeword(Codeword(w, F)); 1156 z := VectorCodeword(Codeword(z, F)); 1157 n := Length(a); 1158 s := Length(w); 1159 if Length(Set(Concatenation(a,w))) <> n + s then 1160 Error("<alpha> and w are not distinct"); 1161 fi; 1162 if ForAny(z,i -> i = Zero(F)) then 1163 Error("<z> must be nonzero"); 1164 fi; 1165 1166 H := []; 1167 for i in List([1..s], index -> List([1..t], vert -> List([1..n], 1168 hor -> z[hor]/(a[hor] - w[index])^vert))) do 1169 Append(H, i); 1170 od; 1171 C := CheckMatCode( BaseMat(VerticalConversionFieldMat(H)), 1172 "generalized Srivastava code", GF(Characteristic(F)) ); 1173 C!.lowerBoundMinimumDistance := s + 1; 1174 return C; 1175end); 1176 1177InstallOtherMethod(GeneralizedSrivastavaCode, "a,w,z,t,q", true, 1178 [IsList, IsList, IsList, IsInt, IsInt], 0, 1179function(a, w, z, t, q) 1180 return GeneralizedSrivastavaCode(a, w, z, t, GF(q)); 1181end); 1182 1183InstallOtherMethod(GeneralizedSrivastavaCode, "a,w,z,t", true, 1184 [IsList, IsList, IsList, IsInt], 0, 1185function(a, w, z, t) 1186 return GeneralizedSrivastavaCode(a, w, z, t, 1187 DefaultField(Concatenation(a,w,z))); 1188end); 1189 1190InstallOtherMethod(GeneralizedSrivastavaCode, "a,w,z,F", true, 1191 [IsList, IsList, IsList, IsField], 0, 1192function(a, w, z, F) 1193 return GeneralizedSrivastavaCode(a, w, z, 1, F); 1194end); 1195 1196InstallOtherMethod(GeneralizedSrivastavaCode, "a, w, z", true, 1197 [IsList, IsList, IsList], 0, 1198function(a, w, z) 1199 return GeneralizedSrivastavaCode(a, w, z, 1, 1200 DefaultField(Concatenation(a, w, z))); 1201end); 1202 1203 1204############################################################################# 1205## 1206#F SrivastavaCode( <a>, <w> [, <mu>] [, <F>] ) . . . . . . . Srivastava code 1207## 1208 1209InstallMethod(SrivastavaCode, "a,w,mu,F", true, 1210 [IsList,IsList,IsInt, IsField], 0, 1211function(a, w, mu, F) 1212 local C, n, s, i, zero, TheMat; 1213 a := VectorCodeword(Codeword(a, F)); 1214 w := VectorCodeword(Codeword(w, F)); 1215 n := Length(a); 1216 s := Length(w); 1217 if Length(Set(Concatenation(a,w))) <> n + s then 1218 Error("the elements of <alpha> and w are not distinct"); 1219 fi; 1220 zero := Zero(F); 1221 for i in [1.. n] do 1222 if a[i]^mu = zero then 1223 Error("z[",i,"] = ",a[i],"^",mu," = ",zero); 1224 fi; 1225 od; 1226 TheMat := List([1..s], j -> List([1..n], i -> a[i]^mu/(a[i] - w[j]) )); 1227 C := CheckMatCode( BaseMat(VerticalConversionFieldMat(TheMat)), 1228 "Srivastava code", GF(Characteristic(F)) ); 1229 C!.lowerBoundMinimumDistance := s + 1; 1230 return C; 1231end); 1232 1233InstallOtherMethod(SrivastavaCode, "a,w,mu,q", true, 1234 [IsList,IsList,IsInt,IsInt], 0, 1235function(a, w, mu, q) 1236 return SrivastavaCode(a, w, mu, GF(q)); 1237end); 1238 1239InstallOtherMethod(SrivastavaCode, "a,w,mu", true, 1240 [IsList,IsList,IsInt], 0, 1241function(a, w, mu) 1242 return SrivastavaCode(a, w, mu, DefaultField(Concatenation(a,w))); 1243end); 1244 1245InstallOtherMethod(SrivastavaCode, "a,w,F", true, 1246 [IsList,IsList,IsField], 0, 1247function(a, w, F) 1248 return SrivastavaCode(a, w, 1, F); 1249end); 1250 1251InstallOtherMethod(SrivastavaCode, "a,w", true, [IsList,IsList], 0, 1252function(a, w) 1253 return SrivastavaCode(a, w, 1, DefaultField(Concatenation(a,w))); 1254end); 1255 1256 1257############################################################################# 1258## 1259#F ExtendedBinaryGolayCode( ) . . . . . . . . . extended binary Golay code 1260## 1261InstallMethod(ExtendedBinaryGolayCode, "only method", true, [], 0, 1262function() 1263 local C; 1264 C := ExtendedCode(BinaryGolayCode()); 1265 C!.name := "extended binary Golay code"; 1266 Unbind( C!.history ); 1267 SetIsCyclicCode(C, false); 1268 SetIsPerfectCode(C, false); 1269 SetIsSelfDualCode(C, true); 1270 C!.boundsCoveringRadius := [ 4 ]; 1271 SetIsNormalCode(C, true); 1272 SetWeightDistribution(C, 1273 [1,0,0,0,0,0,0,0,759,0,0,0,2576,0,0,0,759,0,0,0,0,0,0,0,1]); 1274 #SetAutomorphismGroup(C, M24); 1275 C!.lowerBoundMinimumDistance := 8; 1276 C!.upperBoundMinimumDistance := 8; 1277 return C; 1278end); 1279 1280 1281############################################################################# 1282## 1283#F ExtendedTernaryGolayCode( ) . . . . . . . . . extended ternary Golay code 1284## 1285InstallMethod(ExtendedTernaryGolayCode, "only method", true, [], 0, 1286function() 1287 local C; 1288 C := ExtendedCode(TernaryGolayCode()); 1289 SetIsCyclicCode(C, false); 1290 SetIsPerfectCode(C, false); 1291 SetIsSelfDualCode(C, true); 1292 C!.boundsCoveringRadius := [ 3 ]; 1293 SetIsNormalCode(C, true); 1294 C!.name := "extended ternary Golay code"; 1295 Unbind( C!.history ); 1296 SetWeightDistribution(C, [1,0,0,0,0,0,264,0,0,440,0,0,24]); 1297 #SetAutomorphismGroup(C, M12); 1298 C!.lowerBoundMinimumDistance := 6; 1299 C!.upperBoundMinimumDistance := 6; 1300 return C; 1301end); 1302 1303 1304############################################################################# 1305## 1306#F BestKnownLinearCode( <n>, <k> [, <F>] ) . returns best known linear code 1307#F BestKnownLinearCode( <rec> ) 1308## 1309## L describs how to create a code. L is a list with two elements: 1310## L[1] is a function and L[2] is a list of arguments for L[1]. 1311## One or more of the argumenst of L[2] may again be such descriptions and 1312## L[2] can be an empty list. 1313## The field .construction contains such a list or false if the code is not 1314## yet in the apropiatelibrary file (/tbl/codeq.g) 1315## 1316 1317InstallMethod(BestKnownLinearCode, "method for bounds/construction record", 1318 true, [IsRecord], 0, 1319function(bds) 1320 local MakeCode, C; 1321 1322 # L describs how to create a code. L is a list with two elements: 1323 # L[1] is a function and L[2] is a list of arguments for L[1]. 1324 # One or more of the argumenst of L[2] may again be such descriptions and 1325 # L[2] can be an empty list. 1326 MakeCode := function(L) 1327 #beware: this is the most beautiful function in GUAVA (according to J) 1328 if IsList(L) and IsBound(L[1]) and IsFunction(L[1]) then 1329 return CallFuncList( L[1], List( L[2], i -> MakeCode(i) ) ); 1330 else 1331 return L; 1332 fi; 1333 end; 1334 1335 if not IsBound(bds.construction) then 1336 bds := BoundsMinimumDistance(bds.n, bds.k, bds.q); 1337 fi; 1338 if bds.construction = false then 1339 Print("Code not yet in library\n"); 1340 return fail; 1341 else 1342 C := MakeCode(bds.construction); 1343 if LowerBoundMinimumDistance(C) > bds.lowerBound then 1344 Print("New table entry found!\n"); 1345 fi; 1346 C!.lowerBoundMinimumDistance := Maximum(bds.lowerBound, 1347 LowerBoundMinimumDistance(C)); 1348 C!.upperBoundMinimumDistance := Minimum(bds.upperBound, 1349 UpperBoundMinimumDistance(C)); 1350 return C; 1351 fi; 1352end); 1353 1354InstallOtherMethod(BestKnownLinearCode, "n, k, q", true, 1355 [IsInt, IsInt, IsInt], 0, 1356function(n, k, q) 1357 local r; 1358 r := BoundsMinimumDistance(n, k, q); 1359 return BestKnownLinearCode(r); 1360end); 1361 1362InstallOtherMethod(BestKnownLinearCode, "n, k, F", true, 1363 [IsInt, IsInt, IsField], 0, 1364function(n, k, F) 1365 local r; 1366 r := BoundsMinimumDistance(n, k, Size(F)); 1367 return BestKnownLinearCode(r); 1368end); 1369 1370InstallOtherMethod(BestKnownLinearCode, "n, k", true, 1371 [IsInt, IsInt], 0, 1372function(n, k) 1373 local r; 1374 r := BoundsMinimumDistance(n, k); 1375 return BestKnownLinearCode(r); 1376end); 1377 1378 1379 1380## Helper functions for cyclic code creation 1381GeneratorMatrixFromPoly := function(p,n) 1382 local coeffs, j, res, r, zero; 1383 coeffs := CoefficientsOfLaurentPolynomial(p); 1384 coeffs := ShiftedCoeffs(coeffs[1], coeffs[2]); 1385 r := DegreeOfLaurentPolynomial(p); 1386 zero := Zero(Field(coeffs)); 1387 res := []; 1388 res[1] := []; 1389 # first row 1390 Append(res[1], coeffs); 1391 Append(res[1], List([r+2..n], i->zero)); 1392 # 2..last-1 1393 if n-r > 2 then 1394 for j in [2..(n-r-1)] do 1395 res[j] := []; 1396 Append(res[j], List([1..j-1], i->zero)); 1397 Append(res[j], coeffs); 1398 Append(res[j], List([r+1+j..n], i->zero)); 1399 od; 1400 fi; 1401 # last row 1402 if n-r > 1 then 1403 res[n-r] := []; 1404 Append(res[n-r], List([1..n-r-1], i->zero)); 1405 Append(res[n-r], coeffs); 1406 fi; 1407 return res; 1408end; 1409 1410 1411CyclicCodeByGenerator := function(F, n, G) 1412 local C, GM; 1413 ## for now, using linear code representation. 1414 ## Get generator matrix and call linear code creation function 1415 ## Note input is a codeword, for consistency. 1416 ## Further note GenMatFromPoly doesn't handle NullPoly case well, 1417 ## so calling NullMat instead. Once using poly rep, this should be 1418 ## unnecessary. 1419 ## And GMFP doesn't handle p = x^n-1 case. 1420 1421 if (G = NullWord(n,F)) or 1422 (VectorCodeword(G) = []) or 1423 (G = Codeword(Indeterminate(F)^n-One(F), F)) then 1424 GM := NullMat(1,n,F); 1425 else 1426 GM := GeneratorMatrixFromPoly(PolyCodeword(G), n); 1427 fi; 1428 C := LinearCodeByGenerators(F, Codeword(One(F)*GM, F)); 1429 #C := LinearCodeByGenerators(F, One(F)*GM); 1430 SetGeneratorMat(C, GM); 1431 SetIsCyclicCode(C, true); 1432 SetSpecialDecoder(C, CyclicDecoder); 1433 return C; 1434end; 1435 1436 1437############################################################################# 1438## 1439#F GeneratorPolCode( <G>, <n> [, <name> ], <F> ) . code from generator poly 1440## 1441 1442InstallMethod(GeneratorPolCode, "Poly, wordlength,name,field", true, 1443 [IsUnivariatePolynomial, IsInt, IsString, IsField], 0, 1444function(G, n, name, F) 1445 local R; 1446 #G := PolyCodeword( Codeword(G, F) ); 1447 if not IsZero(G) then 1448 G := Gcd(G,(Indeterminate(F)^n-One(F))); 1449 fi; 1450 R := CyclicCodeByGenerator(F, n, Codeword(G, F)); 1451 SetGeneratorPol(R, G); 1452 R!.name := name; 1453 return R; 1454end); 1455 1456InstallOtherMethod(GeneratorPolCode, "Poly,wordlength,field", true, 1457 [IsUnivariatePolynomial, IsInt, IsField], 0, 1458function(G, n, F) 1459 return GeneratorPolCode(G,n,"code defined by generator polynomial", F); 1460end); 1461 1462InstallOtherMethod(GeneratorPolCode, "Poly,wordlength,name,fieldsize", true, 1463 [IsUnivariatePolynomial, IsInt, IsString, IsInt], 0, 1464function(G, n, name, q) 1465 return GeneratorPolCode(G, n, name, GF(q)); 1466end); 1467 1468InstallOtherMethod(GeneratorPolCode, "Poly, wordlength, fieldsize", true, 1469 [IsUnivariatePolynomial, IsInt, IsInt], 0, 1470function(G, n, q) 1471 return GeneratorPolCode(G, n, "code defined by generator polynomial", 1472 GF(q)); 1473end); 1474 1475 1476############################################################################# 1477## 1478#F CheckPolCode( <H>, <n> [, <name> ], <F> ) . . code from check polynomial 1479## 1480 1481InstallMethod(CheckPolCode, "check poly, wordlength, name, field", true, 1482 [IsUnivariatePolynomial, IsInt, IsString, IsField], 0, 1483function(H, n, name, F) 1484 local R,G; 1485 H := PolyCodeword( Codeword(H, F) ); 1486 H := Gcd(H, (Indeterminate(F)^n-One(F))); 1487 # Get generator pol 1488 G := EuclideanQuotient((Indeterminate(F)^n-One(F)), H); 1489 1490 R := CyclicCodeByGenerator(F, n, Codeword(G,F)); 1491 SetCheckPol(R, H); 1492 SetGeneratorPol(R, G); 1493 R!.name := name; 1494 return R; 1495end); 1496 1497InstallOtherMethod(CheckPolCode, "check poly, wordlength, field", true, 1498 [IsUnivariatePolynomial, IsInt, IsField], 0, 1499function(H, n, F) 1500 return CheckPolCode(H, n, "code defined by check polynomial", F); 1501end); 1502 1503InstallOtherMethod(CheckPolCode, "check poly, wordlength, name, fieldsize", 1504 true, [IsUnivariatePolynomial, IsInt, IsString, IsInt], 0, 1505function(H, n, name, q) 1506 return CheckPolCode(H, n, name, GF(q)); 1507end); 1508 1509InstallOtherMethod(CheckPolCode, "check poly, wordlength, fieldsize", true, 1510 [IsUnivariatePolynomial, IsInt, IsInt], 0, 1511function(H, n, q) 1512 return CheckPolCode(H, n, "code defined by check polynomial", GF(q)); 1513end); 1514 1515 1516############################################################################# 1517## 1518#F RepetitionCode( <n> [, <F>] ) . . . . . . . repetition code of length <n> 1519## 1520 1521InstallMethod( RepetitionCode, "wordlength, Field", true, [IsInt, IsField], 0, 1522function(n, F) 1523 local C, q, wd, p; 1524 q :=Size(F); 1525 p := LaurentPolynomialByCoefficients( ElementsFamily(FamilyObj(F)), 1526 List([1..n], t->One(F)), 0); 1527 C := GeneratorPolCode(p, n, "repetition code", F); 1528 C!.lowerBoundMinimumDistance := n; 1529 C!.upperBoundMinimumDistance := n; 1530 if n = 2 and q = 2 then 1531 SetIsSelfDualCode(C, true); 1532 else 1533 SetIsSelfDualCode(C, false); 1534 fi; 1535 wd := NullVector(n+1); 1536 wd[1] := 1; 1537 wd[n+1] := q-1; 1538 SetWeightDistribution(C, wd); 1539 if n < 260 then 1540 SetAutomorphismGroup(C, SymmetricGroup(n)); 1541 fi; 1542 if F = GF(2) then 1543 SetIsNormalCode(C, true); 1544 fi; 1545 if (n mod 2 = 0) or (F <> GF(2)) then 1546 SetIsPerfectCode(C, false); 1547 C!.boundsCoveringRadius := [ Minimum(n-1,QuoInt((q-1)*n,q)) ]; 1548 else 1549 C!.boundsCoveringRadius := [ QuoInt(n,2) ]; 1550 SetIsPerfectCode(C, true); 1551 fi; 1552 return C; 1553end); 1554 1555InstallOtherMethod(RepetitionCode, "wordlength, fieldsize", true, 1556 [IsInt, IsInt], 0, 1557function(n, q) 1558 return RepetitionCode(n, GF(q)); 1559end); 1560 1561InstallOtherMethod(RepetitionCode, "wordlength", true, [IsInt], 0, 1562function(n) 1563 return RepetitionCode(n, GF(2)); 1564end); 1565 1566 1567############################################################################# 1568## 1569#F WholeSpaceCode( <n> [, <F>] ) . . . . . . . . . . returns <F>^<n> as code 1570## 1571 1572InstallMethod(WholeSpaceCode, "wordlength, Field", true, [IsInt, IsField], 0, 1573function(n, F) 1574 local C, index, q; 1575 C := GeneratorPolCode( Indeterminate(F)^0, n, "whole space code", F); 1576 C!.lowerBoundMinimumDistance := 1; 1577 C!.upperBoundMinimumDistance := 1; 1578 SetAutomorphismGroup(C, SymmetricGroup(n)); 1579 C!.boundsCoveringRadius := [ 0 ]; 1580 if F = GF(2) then 1581 SetIsNormalCode(C, true); 1582 fi; 1583 SetIsPerfectCode(C, true); 1584 SetIsSelfDualCode(C, false); 1585 q := Size(F) - 1; 1586 SetWeightDistribution(C, List([0..n], i-> q^i*Binomial(n, i)) ); 1587 return C; 1588end); 1589 1590InstallOtherMethod(WholeSpaceCode, "wordlength, fieldsize", true, 1591 [IsInt, IsInt], 0, 1592function(n, q) 1593 return WholeSpaceCode(n, GF(q)); 1594end); 1595 1596InstallOtherMethod(WholeSpaceCode, "wordlength", true, [IsInt], 0, 1597function(n) 1598 return WholeSpaceCode(n, GF(2)); 1599end); 1600 1601 1602############################################################################# 1603## 1604#F CyclicCodes( <n> ) . . returns a list of all cyclic codes of length <n> 1605## 1606 1607InstallMethod(CyclicCodes, "wordlength, Field", true, 1608 [IsInt, IsField], 0, 1609function(n, F) 1610 local f, Pl; 1611 f := Factors(Indeterminate(F) ^ n - One(F)); 1612 Pl := List(Combinations(f), c->Product(c)*Indeterminate(F)^0); 1613 return List(Pl, p->GeneratorPolCode(p, n, "enumerated code", F)); 1614end); 1615 1616InstallOtherMethod(CyclicCodes, "wordlength, k, Field", true, 1617 [IsInt, IsInt, IsField], 0, 1618function(n, k, F) 1619 local r, f, Pl, codes; 1620 r := n - k; 1621 f := Factors(Indeterminate(F)^n-One(F)); 1622 Pl := []; 1623 1624 codes := function(f, g) 1625 local i, tempf; 1626 if Degree(g) < r then 1627 i := 1; 1628 while i <= Length(f) and Degree(g)+Degree(f[i][1]) <= r do 1629 if f[i][2] = 1 then 1630 tempf := f{[ i+1 .. Length(f) ]}; 1631 else 1632 tempf := ShallowCopy( f ); 1633 tempf[i][2] := f[i][2] - 1; 1634 fi; 1635 codes( tempf, g * f[i][1] ); 1636 i := i + 1; 1637 od; 1638 elif Degree(g) = r then 1639 Add( Pl, g ); 1640 fi; 1641 end; 1642 1643 codes( Collected( f ), Indeterminate(F)^0 ); 1644 return List(Pl, p->GeneratorPolCode(p,n,"enumerated code",F)); 1645end); 1646 1647InstallOtherMethod(CyclicCodes, "wordlength, fieldsize", true, 1648 [IsInt, IsInt], 0, 1649function(n, q) 1650 return CyclicCodes(n, GF(q)); 1651end); 1652 1653InstallOtherMethod(CyclicCodes, "wordlength, k, fieldsize", true, 1654 [IsInt, IsInt, IsInt], 0, 1655function(n, k, q) 1656 return CyclicCodes(n, k, GF(q)); 1657end); 1658 1659 1660############################################################################# 1661## 1662#F NrCyclicCodes( <n>, <F>) . . . number of cyclic codes of length <n> 1663## 1664 1665InstallMethod(NrCyclicCodes, "wordlength, Field", true, [IsInt, IsField], 0, 1666function(n, F) 1667 return NrCombinations(Factors(Indeterminate(F)^n-One(F))); 1668end); 1669 1670InstallOtherMethod(NrCyclicCodes, "wordlength, fieldsize", true, 1671 [IsInt, IsInt], 0, 1672function(n, q) 1673 return NrCyclicCodes(n, GF(q)); 1674end); 1675 1676 1677############################################################################# 1678## 1679#F BCHCode( <n> [, <b>], <delta> [, <F>] ) . . . . . . . . . . . . BCH code 1680## 1681## BCHCode (n [, b ], delta [, F]) returns the BCH code over F with 1682## wordlength n, designedDistance delta, constructed from powers 1683## x^b, x^(b+1), ..., x^(b+delta-2), where x is a primitive n'th power root 1684## of unity; b = 1 by default; the function returns a narrow sense BCH code 1685## Gcd(n,q) = 1 and 2<=delta<=n-b+1 1686 1687InstallMethod(BCHCode, "wordlength, start, designed distance, fieldsize", true, 1688 [IsInt, IsInt, IsInt, IsInt], 0, 1689function(n, start, del, q) 1690 local stop, m, b, C, test, Cyclo, PowerSet, t, 1691 zero, desdist, G, superfl, i, BCHTable; 1692 1693 BCHTable := [ [31,11,11], [63,36,11], [63,30,13], [127,92,11], 1694 [127,85,13], [255,223,9], [255,215,11], [255,207,13], 1695 [255,187,19], [255,171,23], [255,155,27], [255,99,47], 1696 [255,79,55], [255,29,95], [255,21,111] ]; 1697########### increase size of table? - wdj 1698 stop := start + del - 2; 1699 if Gcd(n,q) <> 1 then 1700 Error ("n and q must be relative primes"); 1701 fi; 1702 zero := Zero(GF(q)); 1703 m := OrderMod(q,n); 1704 b := PrimitiveUnityRoot(q,n); 1705 PowerSet := [start..stop]; 1706 G := Indeterminate(GF(q))^0; 1707 while Length(PowerSet) > 0 do 1708 test := PowerSet[1] mod n; ############### changed 8-1-2004 1709 G := G * MinimalPolynomial(GF(q), b^test); 1710 t := (q*test) mod n; 1711 while t <> (test mod n) do ############### changed 7-31-2004 1712 RemoveSet(PowerSet, t); 1713 t := (q*t) mod n; 1714 od; 1715 RemoveSet(PowerSet, test); 1716 od; 1717###################################### 1718###### This loop computes the product of the 1719###### min polys of the elements when it should only 1720###### compute the lcm of them 1721###### Why not use cyclotomic codests and roots of code? 1722###################################### 1723 C := GeneratorPolCode(G, n, GF(q)); 1724########### this removes extra powers by taking the gcd with x^n-1 1725 SetSpecialDecoder(C, BCHDecoder); 1726########### this overwrites SetSpecialDecoder(C, CyclicDecoder); 1727 # Calculate Bose distance: 1728 Cyclo := CyclotomicCosets(q,n); 1729######## why not do this in the above loop? 1730 PowerSet := []; 1731 for t in [start..stop] do 1732 for test in [1..Length(Cyclo)] do 1733 if (t mod n) in Cyclo[test] then ##### added 8-1-2004 1734 AddSet(PowerSet, Cyclo[test]); 1735 fi; 1736 od; 1737 od; 1738 PowerSet := Flat(PowerSet); 1739 while stop + 1 in PowerSet do 1740 stop := stop + 1; 1741 od; 1742 while start - 1 in PowerSet do 1743 start := start - 1; 1744 od; 1745 desdist := stop - start + 2; 1746 if desdist > n then 1747 Error("invalid designed distance"); 1748 fi; 1749 # In some cases the true minimumdistance is known: 1750 SetDesignedDistance(C, desdist); 1751 C!.lowerBoundMinimumDistance := desdist; 1752 if (q=2) and (n mod desdist = 0) and (start = 1) then 1753 C!.upperBoundMinimumDistance := desdist; 1754 elif q=2 and desdist mod 2 = 0 and (n=2^m - 1) and (start=1) and 1755 (Sum(List([0..QuoInt(desdist-1, 2) + 1], i -> Binomial(n, i))) > 1756 (n + 1) ^ QuoInt(desdist-1, 2)) then 1757 C!.upperBoundMinimumDistance := desdist; 1758 elif (n = q^m - 1) and (desdist = q^OrderMod(q,desdist) - 1) 1759 and (start=1) then 1760 C!.upperBoundMinimumDistance := desdist; 1761 fi; 1762 # Look up this code in the table 1763########## table only for q=2^r, add this to if statement 1764########## to speed this up ?? - wdj 1765 if start=1 then 1766 for i in BCHTable do 1767 if i[1] = n and i[2] = Dimension(C) then 1768 C!.lowerBoundMinimumDistance := i[3]; 1769 C!.upperBoundMinimumDistance := i[3]; 1770 fi; 1771 od; 1772 fi; 1773 # Calculate minimum of q*desdist - 1 for primitive n.s. BCH code 1774 if q^m - 1 = n and start = 1 then 1775 PowerSet := [start..stop]; 1776 superfl := true; 1777 i := PowerSet[Length(PowerSet)] * q mod n; 1778 while superfl do 1779 while i <> PowerSet[Length(PowerSet)] and not i in PowerSet do 1780 i := i * q mod n; 1781 od; 1782 if i = PowerSet[Length(PowerSet)] then 1783 superfl := false; 1784 else 1785 PowerSet := PowerSet{[1..Length(PowerSet)-1]}; 1786 i := PowerSet[Length(PowerSet)] * q mod n; 1787 fi; 1788 od; 1789 C!.upperBoundMinimumDistance := Minimum(UpperBoundMinimumDistance(C), 1790 q * (Length(PowerSet) + 1) - 1); 1791 fi; 1792 C!.name := Concatenation("BCH code, delta=", 1793 String(desdist), ", b=", String(start)); 1794 return C; 1795end); 1796 1797InstallOtherMethod(BCHCode, "wordlength, start, designed distance, Field", 1798 true, [IsInt, IsInt, IsInt, IsField], 0, 1799function(n, b, del, F) 1800 return BCHCode(n, b, del, Size(F)); 1801end); 1802 1803InstallOtherMethod(BCHCode, "wordlength, designed distance", true, 1804 [IsInt, IsInt], 0, 1805function(n, del) 1806 return BCHCode(n, 1, del, 2); 1807end); 1808 1809InstallOtherMethod(BCHCode, "wordlength, start, designed distance", true, 1810 [IsInt, IsInt, IsInt], 0, 1811function(n, b, del) 1812 return BCHCode(n, b, del, 2); 1813end); 1814 1815InstallOtherMethod(BCHCode, "wordlength, designed distance, field", true, 1816 [IsInt, IsInt, IsField], 0, 1817function(n, del, F) 1818 return BCHCode(n, 1, del, Size(F)); 1819end); 1820 1821 1822############################################################################# 1823## 1824#F ReedSolomonCode( <n>, <d> ) . . . . . . . . . . . . . . Reed-Solomon code 1825## 1826## ReedSolomonCode (n, d) returns a primitive narrow sense BCH code with 1827## wordlength n, over alphabet q = n+1, designed distance d 1828 1829InstallMethod(ReedSolomonCode, "wordlength, designed distance", true, 1830 [IsInt, IsInt], 0, 1831function(n, d) 1832 local C,b,q,wd,w; 1833 q := n+1; 1834 if not IsPrimePowerInt(q) then 1835 Error("q = n+1 must be a prime power"); 1836 fi; 1837 b := Z(q); 1838 1839 C := GeneratorPolCode(Product([1..d-1], i-> (Indeterminate(GF(q))-b^i)), n, 1840 "Reed-Solomon code", GF(q) ); 1841 SetRootsOfCode(C, List([1..d-1], i->b^i)); 1842 C!.lowerBoundMinimumDistance := d; 1843 C!.upperBoundMinimumDistance := d; 1844 SetDesignedDistance(C, d); 1845 SetSpecialDecoder(C, BCHDecoder); 1846 IsMDSCode(C); # Calculate weightDistribution field 1847 return C; 1848end); 1849 1850InstallMethod(ExtendedReedSolomonCode, "wordlength, designed distance", true, 1851 [IsInt, IsInt], 0, 1852function(n, d) 1853 local i, j, s, C, G, Ce; 1854 C := ReedSolomonCode(n-1, d-1); 1855 G := MutableCopyMat( GeneratorMat(C) ); 1856 TriangulizeMat(G); 1857 for i in [1..Size(G)] do; 1858 s := 0; 1859 for j in [1..Size(G[i])] do; 1860 s := s + G[i][j]; 1861 od; 1862 Append(G[i], [-s]); 1863 od; 1864 Ce := GeneratorMatCodeNC(G, LeftActingDomain(C)); 1865 Ce!.name := "extended Reed Solomon code"; 1866 Ce!.lowerBoundMinimumDistance := d; 1867 Ce!.upperBoundMinimumDistance := d; 1868 IsMDSCode(Ce); 1869 return Ce; 1870end); 1871 1872## RootsCode implementation expunged and rewritten for Guava 3.11 1873## J. E. Fields 1/15/2012 (with help from WDJ) 1874############################################################################# 1875## 1876#F RootsCode( <n>, <list>, <field> ) code constructed from roots of polynomial 1877## 1878## RootsCode (n, rootlist, F) or RootsCode (n, <powerlist>, fieldsize) or 1879## RootsCode(n, rootlist) returns the 1880## code with generator polynomial equal to the least common multiple of 1881## the minimal polynomials of the n'th roots of unity in the list. 1882## The code has wordlength n. 1883## 1884 1885InstallMethod(RootsCode, "method for n, rootlist, field", true, [IsInt, IsList, IsField], 0, 1886function(n, L, F) 1887 local g, C, num, power, p, q, z, i, j, rootslist, powerlist, max, cc, CC, CCz; 1888 L := Set(L); 1889 q := Size(Field(L)); 1890 z := Z(q); 1891 g:=One(F); 1892 if List(L, i->i^n) <> NullVector(Length(L), F) + z^0 then 1893 Error("powers must all be n'th roots of unity"); 1894 fi; 1895 p := Characteristic(Field(L)); 1896 CC:=CyclotomicCosets(p,q-1); 1897 CCz:=List([1..Length(CC)],i->List(CC[i],j->z^j)); 1898 ## this is the set of cyclotomic cosets, represented 1899 ## as powers of a primitive element z 1900 powerlist := []; 1901 for i in [1..Length(L)] do 1902 for cc in CCz do 1903 if L[i] in cc then 1904 Append(powerlist,[cc]); 1905 fi; 1906 od; 1907 od; 1908 ########### add a cyclotomic coset into powerlist 1909 ########### if there is an element of L in that coset 1910 g:=Product([1..Length(powerlist)],i->MinimalPolynomial(F,powerlist[i][1])); 1911 C:=GeneratorPolCode(g,n,"code defined by roots",F); 1912 1913 rootslist := []; 1914 for cc in powerlist do 1915 Append(rootslist, cc); 1916 od; 1917 rootslist := Set(rootslist); 1918 SetRootsOfCode(C, rootslist); 1919 1920 # Find the largest number of successive powers for BCH bound 1921 max := 1; 1922 i := 1; 1923 num := Length(powerlist); 1924 for z in [2..num] do 1925 if powerlist[z] <> powerlist[i] + z-i then 1926 max := Maximum(max, z - i); 1927 i := z; 1928 fi; 1929 od; 1930 C!.lowerBoundMinimumDistance := Maximum(max, num+1 - i) + 1; 1931return C; 1932end); 1933 1934 1935InstallOtherMethod(RootsCode, "method for n, rootlist", true, [IsInt, IsList], 0, 1936function(n, L) 1937local p,q,z,L1,L2,F; 1938 L1 := Set(L); 1939 p := Characteristic(Field(L1)); 1940 z := Z(Size(Field(L1))); 1941 F := GF(p); 1942 #L2 := List([1..Length(L1)],i-> LogFFE(L[i],z)); 1943 #Print(L1, p, z, F, L2); 1944 return RootsCode(n, L1, F); 1945end); 1946 1947InstallOtherMethod(RootsCode, "method for n, powerlist, fieldsize", true, [IsInt, IsList, IsInt], 0, 1948function(n, P, q) 1949# , "method for n, powerlist, fieldsize", true, [IsInt, IsList, IsInt], 0, 1950local z, L; 1951 z := PrimitiveUnityRoot(q,n); 1952 L := Set(List(P, i->z^i)); 1953 return RootsCode(n, L, GF(q)); 1954end); 1955 1956############################################################################# 1957## 1958#F QRCode( <n> [, <F>] ) . . . . . . . . . . . . . . quadratic residue code 1959## 1960 1961InstallMethod(QRCode, "modulus, fieldsize", true, [IsInt, IsInt], 0, 1962function(n, q) 1963 local m, b, Q, N, t, g, lower, upper, C, F, coeffs; 1964 if Jacobi(q,n) <> 1 then 1965 Error("q must be a quadratic residue modulo n"); 1966 elif not IsPrimeInt(n) then 1967 Error("n must be a prime"); 1968 elif not IsPrimeInt(q) then 1969 Error("q must be a prime"); 1970 fi; 1971 m := OrderMod(q,n); 1972 F := GF(q^m); 1973 b := PrimitiveUnityRoot(q,n); 1974 Q := []; 1975 N := [1..n]; 1976 for t in [1..n-1] do 1977 AddSet(Q, t^2 mod n); 1978 od; 1979 for t in Q do 1980 RemoveSet(N, t); 1981 od; 1982 g := Product(Q, 1983 i -> LaurentPolynomialByCoefficients(ElementsFamily(FamilyObj(F)), 1984 [-b^i, b^0], 0) ); 1985 coeffs := CoefficientsOfLaurentPolynomial(g); 1986 coeffs := ShiftedCoeffs(coeffs[1], coeffs[2]); 1987 C := GeneratorPolCode( LaurentPolynomialByCoefficients( 1988 ElementsFamily(FamilyObj(GF(q))), coeffs, 0), 1989 n, "quadratic residue code", GF(q) ); 1990 if RootInt(n)^2 = n then 1991 lower := RootInt(n); 1992 else 1993 lower := RootInt(n)+1; 1994 fi; 1995 if n mod 4 = 3 then 1996 while lower^2-lower+1 < n do 1997 lower := lower + 1; 1998 od; 1999 fi; 2000 if (n mod 8 = 7) and (q = 2) then 2001 while lower mod 4 <> 3 do 2002 lower := lower + 1; 2003 od; 2004 fi; 2005 upper := Weight(Codeword(coeffs)); 2006 C!.lowerBoundMinimumDistance := lower; 2007 C!.upperBoundMinimumDistance := upper; 2008 return C; 2009end); 2010 2011InstallOtherMethod(QRCode, "modulus, field", true, [IsInt, IsField], 0, 2012function(n, F) 2013 return QRCode(n, Size(F)); 2014end); 2015 2016InstallOtherMethod(QRCode, "modulus", true, [IsInt], 0, 2017function(n) 2018 return QRCode(n, 2); 2019end); 2020 2021############################################################################# 2022## 2023#F QQRCode( <n> [, <F>] ) . . . . . . . . binary quasi-quadratic residue code 2024## 2025## Code of Bazzi-Mittel (see Bazzi, L. and Mitter, S.K. "Some constructions of 2026## codes from group actions" preprint March 2003 (submitted to IEEE IT) 2027## 2028 2029InstallMethod(QQRCode, "Modulus", true, [IsInt], 0, 2030function(p) 2031local G0,G,Q,N,QN,i,C,F,QuadraticResidueSupports,NonQuadraticResidueSupports; 2032 2033# start local functions: 2034QuadraticResidueSupports:=function(p) 2035 local L,i; 2036 L:=List([1..(p-1)],i->1+Legendre(i,p))/2; 2037 return L; 2038 end; 2039 NonQuadraticResidueSupports:=function(p) 2040 local L,i; 2041 L:=List([1..(p-1)],i->1-Legendre(i,p))/2; 2042 return L; 2043 end; 2044# end local functions 2045 2046 F:=GF(2); 2047 Q:=[Zero(F)]; 2048 Q:=One(F)*Concatenation(Q,QuadraticResidueSupports(p)); 2049 N:=[Zero(F)]; 2050 N:=One(F)*Concatenation(N,NonQuadraticResidueSupports(p)); 2051 QN:=Concatenation(Q,N); 2052 G:=BlockMatrix([[1,1,CirculantMatrix(p,Q)],[1,2,CirculantMatrix(p,N)]],1,2); 2053 if p mod 4 = 1 then 2054 G0:=List([1..(p-1)],i->G[i]); 2055 else 2056 G0:=G; 2057 fi; 2058 C:=GeneratorMatCode(G0,F); 2059 C!.DoublyCirculant:=G0; 2060 C!.GeneratorMat:=G0; 2061 return C; 2062end); 2063 2064############################################################################# 2065## 2066#F QQRCodeNC( <n> [, <F>] ) . . . . . . . . binary quasi-quadratic residue code 2067## 2068## Code of Bazzi-Mittel (see Bazzi, L. and Mitter, S.K. "Some constructions of 2069## codes from group actions" preprint March 2003 (submitted to IEEE IT) 2070## 2071 2072InstallMethod(QQRCodeNC, "Modulus", true, [IsInt], 0, 2073function(p) 2074local G,G0,Q,N,QN,i,C,F,QuadraticResidueSupports,NonQuadraticResidueSupports; 2075 2076# start local functions: 2077QuadraticResidueSupports:=function(p) 2078 local L,i; 2079 L:=List([1..(p-1)],i->1+Legendre(i,p))/2; 2080 return L; 2081 end; 2082 NonQuadraticResidueSupports:=function(p) 2083 local L,i; 2084 L:=List([1..(p-1)],i->1-Legendre(i,p))/2; 2085 return L; 2086 end; 2087# end local functions 2088 2089 F:=GF(2); 2090 Q:=[Zero(F)]; 2091 #Q:=[]; 2092 Q:=One(F)*Concatenation(Q,QuadraticResidueSupports(p)); 2093 N:=[Zero(F)]; 2094 #N:=[]; 2095 N:=One(F)*Concatenation(N,NonQuadraticResidueSupports(p)); 2096 #QN:=Concatenation(Q,N); 2097 G:=BlockMatrix([[1,1,CirculantMatrix(p,N)],[1,2,CirculantMatrix(p,Q)]],1,2); 2098 if p mod 4 = 1 then 2099 G0:=List([1..(p-1)],i->G[i]); 2100 else 2101 G0:=G; 2102 fi; 2103 C:=GeneratorMatCodeNC(G0,F); 2104 C!.DoublyCirculant:=G0; 2105 C!.GeneratorMat:=G0; 2106 return C; 2107end); 2108 2109############################################################################# 2110## 2111#F NullCode( <n> [, <F>] ) . . . . . . . . . . . code consisting only of <0> 2112## 2113 2114InstallMethod(NullCode, "wordlength, field", true, [IsInt, IsField], 0, 2115function(n, F) 2116 local C; 2117 C := ElementsCode([NullWord(n,F)], "nullcode", F); 2118 C!.lowerBoundMinimumDistance := n; 2119 C!.upperBoundMinimumDistance := n; 2120 SetMinimumDistance(C, n); 2121 SetWeightDistribution(C, Concatenation([1], NullVector(n))); 2122 SetAutomorphismGroup(C, SymmetricGroup(n)); 2123 C!.boundsCoveringRadius := [ n ]; 2124 SetCoveringRadius(C, n); 2125 IsCyclicCode(C); # will set all basic linear and cyclic code info 2126 return C; 2127end); 2128 2129InstallOtherMethod(NullCode, "wordlength, fieldsize", true, [IsInt, IsInt], 0, 2130function(n, q) 2131 return NullCode(n, GF(q)); 2132end); 2133 2134InstallOtherMethod(NullCode, "wordlength", true, [IsInt], 0, 2135function(n) 2136 return NullCode(n, GF(2)); 2137end); 2138 2139 2140############################################################################# 2141## 2142#F FireCode( <G>, <b> ) . . . . . . . . . . . . . . . . . . . . . Fire code 2143## 2144## FireCode (G, b) constructs the Fire code that is capable of correcting any 2145## single error burst of length b or less. 2146## G is a primitive polynomial of degree m 2147## 2148 2149InstallMethod(FireCode, "poly, burstlength", true, 2150 [IsUnivariatePolynomial, IsInt], 0, 2151function(G, b) 2152 local m, C, n; 2153 G := PolyCodeword(Codeword(G)); 2154 if CoefficientsRing(DefaultRing(G)) <> GF(2) then 2155 Error("polynomial must be over GF(2)"); 2156 fi; 2157 if Length(Factors(G)) <> 1 then 2158 Error("polynomial G must be primitive"); 2159 fi; 2160 m := DegreeOfLaurentPolynomial(G); 2161 n := Lcm(2^m-1,2*b-1); 2162 C := GeneratorPolCode( G*(Indeterminate(GF(2))^(2*b-1) + One(GF(2))), n, 2163 Concatenation(String(b), " burst error correcting fire code"), 2164 GF(2) ); 2165 return C; 2166end); 2167 2168############################################################################# 2169## 2170#F BinaryGolayCode( ) . . . . . . . . . . . . . . . . . . binary Golay code 2171## 2172InstallMethod(BinaryGolayCode, "only method", true, [], 0, 2173function() 2174 local p,C; 2175 p := LaurentPolynomialByCoefficients( 2176 ElementsFamily(FamilyObj(GF(2))), 2177 Z(2)^0*[1,0,1,0,1,1,1,0,0,0,1,1], 0); 2178 C := CyclicCodeByGenerator(GF(2), 23, Codeword(p)); 2179 SetGeneratorPol(C, p); 2180 SetDimension(C, 12); 2181 SetRedundancy(C, 11); 2182 SetSize(C, 2^12); 2183 C!.name := "binary Golay code"; 2184 C!.lowerBoundMinimumDistance := 7; 2185 C!.upperBoundMinimumDistance := 7; 2186 SetWeightDistribution(C, 2187 [1,0,0,0,0,0,0,253,506,0,0,1288,1288,0,0,506,253,0,0,0,0,0,0,1]); 2188 C!.boundsCoveringRadius := [ 3 ]; 2189 SetIsNormalCode(C, true); 2190 SetIsPerfectCode(C, true); 2191 return C; 2192end); 2193 2194 2195############################################################################# 2196## 2197#F TernaryGolayCode( ) . . . . . . . . . . . . . . . . . ternary Golay code 2198## 2199InstallMethod(TernaryGolayCode, "only method", true, [], 0, 2200function() 2201 local p, C; 2202 p := LaurentPolynomialByCoefficients(ElementsFamily(FamilyObj(GF(3))), 2203 Z(3)^0*[2,0,1,2,1,1], 0); 2204 C := CyclicCodeByGenerator(GF(3), 11, Codeword(p)); 2205 C!.name := "ternary Golay code"; 2206 SetGeneratorPol(C, p); 2207 SetRedundancy(C, 5); 2208 SetDimension(C, 6); 2209 SetSize(C, 3^6); 2210 C!.lowerBoundMinimumDistance := 5; 2211 C!.upperBoundMinimumDistance := 5; 2212 SetWeightDistribution(C, [1,0,0,0,0,132,132,0,330,110,0,24]); 2213 SetIsNormalCode(C, true); 2214 SetIsPerfectCode(C, true); 2215 return C; 2216end); 2217 2218 2219############################################################################# 2220## 2221#F EvaluationCode( <P>, <L>, <R> ) 2222## 2223## P is a list of n points in F^r 2224## L is a list of rational functions in r variables 2225## EvaluationCode returns the image of the evaluation map f->[f(P1),...,f(Pn)], 2226## as f ranges over the vector space of functions spanned by L. 2227## The output is the code whose generator matrix has rows (f(P1)...f(Pn)) where 2228## f is in L and P={P1,..,Pn} 2229## 2230InstallMethod(EvaluationCode,"points, basis functions, multivariate poly ring", true, 2231[IsList, IsList, IsRing], 0, 2232function(P,L,R) 2233 local i, G, n, C, j, k, varsn,varsd, vars, F, ValueExtended; 2234####################### 2235 ValueExtended:=function(f,vars,pt) 2236 local df, nf; 2237 df:=DenominatorOfRationalFunction(f*vars[1]^0); 2238 nf:=NumeratorOfRationalFunction(f*vars[1]^0); 2239 #if (varsn=[] and varsd=[]) then return f; fi; 2240 return Value(nf,vars,pt)*Value(df,vars,pt)^(-1); 2241 end; 2242####################### 2243 vars:=IndeterminatesOfPolynomialRing(R); 2244 F:=CoefficientsRing(R); 2245 n:=Length(P); 2246 k:=Length(L); 2247 G:=ShallowCopy(NullMat(k,n,F)); 2248 for i in [1..k] do 2249 for j in [1..n] do 2250 G[i][j]:=ValueExtended(L[i],vars,P[j]); 2251 od; 2252 od; 2253 C:=GeneratorMatCode(G," evaluation code",F); 2254 C!.EvaluationMat:=ShallowCopy(G); 2255 C!.basis:=L; 2256 C!.points:=P; 2257 C!.ring:=R; 2258 return C; 2259end); 2260 2261############################################################################# 2262## 2263#F GeneralizedReedSolomonCode( <P>, <k>, <R> ) 2264## 2265## P is a list of n points in F 2266## k is an integer 2267## GRSCode returns the image of the evaluation map f->[f(P1),...,f(Pn)], 2268## as f ranges over the vector space of polynomials in 1 variable 2269## of degree < k in the ring R. 2270## The output is the code whose generator matrix has rows (f(P1)...f(Pn)) where 2271## f = x^j, j<k, and P={P1,..,Pn} 2272## 2273InstallMethod(GeneralizedReedSolomonCode,"points, basis functions, univariate poly ring", true, 2274[IsList, IsInt, IsRing], 0, 2275function(P,k,R) 2276local p, L, vars, i, F, f, x, C, R0, P0, G, j, n; 2277 n:=Length(P); 2278 F:=CoefficientsRing(R); 2279 G:=NullMat(k,n,F); 2280 vars:=IndeterminatesOfPolynomialRing(R); 2281 x:=vars[1]; 2282 L:=List([0..(k-1)],i->(x^i)); 2283 for i in [1..k] do 2284 for j in [1..n] do 2285 G[i][j]:=Value(L[i],vars,[P[j]]); 2286 od; 2287 od; 2288 C:=GeneratorMatCode(G," generalized Reed-Solomon code",F); 2289 C!.GeneratorMat:=ShallowCopy(G); 2290 C!.degree:=k; 2291 C!.points:=P; 2292 C!.ring:=R; 2293 SetSpecialDecoder(C, GeneralizedReedSolomonDecoder); 2294 return C; 2295end); 2296 2297############################################################################# 2298## 2299#F GeneralizedReedSolomonCode( <P>, <k>, <R> , <wts> ) 2300## 2301## P is a list of n points in F 2302## k is an integer 2303## GRSCode returns the image of the evaluation map f->[f(P1),...,f(Pn)], 2304## as f ranges over the vector space of polynomials in 1 variable 2305## of degree < k in the ring R. 2306## The output is the code whose generator matrix has rows (w1*f(P1),...,wn*f(Pn)) where 2307## f = x^j, j<k, P={P1,..,Pn}, wts=[w1,...,wn] \in F^n 2308## 2309InstallOtherMethod(GeneralizedReedSolomonCode,"points, basis functions, univariate poly ring", true, 2310[IsList, IsInt, IsRing, IsList], 0, 2311function(P,k,R,wts) 2312local p, L, vars, i, F, f, x, C, R0, P0, G, j, n; 2313 n:=Length(P); 2314 F:=CoefficientsRing(R); 2315 G:=NullMat(k,n,F); 2316 vars:=IndeterminatesOfPolynomialRing(R); 2317 x:=vars[1]; 2318 L:=List([0..(k-1)],i->(x^i)); 2319 for i in [1..k] do 2320 for j in [1..n] do 2321 G[i][j]:=wts[j]*Value(L[i],vars,[P[j]]); 2322 od; 2323 od; 2324 C:=GeneratorMatCode(G," weighted generalized Reed-Solomon code",F); 2325 C!.GeneratorMat:=ShallowCopy(G); 2326 C!.degree:=k; 2327 C!.points:=P; 2328 C!.weights:=wts; 2329 C!.ring:=R; 2330#SetSpecialDecoder(C, GeneralizedReedSolomonDecoder); 2331 return C; 2332end); 2333 2334############################################################################# 2335## 2336#F OnePointAGCode( <crv>, <pts>, <m>, <R> ) 2337## 2338## R = F[x,y] is a polynomial ring over a finite field F 2339## crv is a polynomial in R representing a plane curve 2340## pts is a list of points on the curve 2341## Computes the AG codes associated to the RR space 2342## L(m*infinity) using Proposition VI.4.1 in Stichtenoth 2343## 2344InstallMethod(OnePointAGCode, 2345"polynomial defining planar curve, points, multiplicity, univariate poly ring", 2346true, [IsPolynomial, IsList, IsInt, IsRing], 0, 2347function(crv,pts,m,R) 2348 local F,f,indets,pt,i,j,G,C,degx,degy,basisLD,xx,yy,allpts; 2349 indets := IndeterminatesOfPolynomialRing(R); 2350 xx:=indets[1]; yy:=indets[2]; 2351 F:=CoefficientsRing(R); 2352 allpts:=AffinePointsOnCurve(crv,R,F); 2353 if not(IsSubset(allpts,pts)) then 2354 Error("The points given must be on the curve\n"); 2355 fi; 2356 degx:=DegreeIndeterminate(crv,xx); 2357 degy:=DegreeIndeterminate(crv,yy); 2358 basisLD:=[]; 2359 for i in [0..(degx-1)] do 2360 for j in [0..m] do 2361 if degx*j+degy*i<m+1 then 2362 basisLD:=Concatenation([xx^i*yy^j],basisLD); 2363 fi; 2364 od; 2365 od; 2366 G:=List(pts,pt->List(basisLD,f->Value(f,indets,pt))); 2367 C:=GeneratorMatCode(TransposedMat(G)," one-point AG code",F); 2368 C!.GeneratorMat:=ShallowCopy(TransposedMat(G)); 2369 C!.multiplicity:=m; 2370 C!.points:=pts; 2371 C!.curve:=crv; 2372 C!.ring:=R; 2373 return C; 2374end); 2375 2376############################################################################# 2377## 2378#F FerreroDesignCode( <P>, <m> ) ... code from a Ferrero design 2379## 2380## 2381#InstallMethod(FerreroDesignCode, 2382#"binary linear code constructed using a Ferrero design", 2383#true, [IsList, IsInt], 0, 2384#function( P,m) 2385# **Requires the GAP package SONATA** 2386# Constructs binary linear code arising from the incdence 2387# matrix of a design associated to a "Ferrero pair" arising 2388# from a fixed-point-free (fpf) automorphism groups and Frobenius group. 2389# The designs that we are looking at (from a Frobenius kernel 2390# of order v and a Frobenius complement of order k) have 2391# v*(v-1)/k distinct blocks and they are all of size k. 2392# Moreover each of the v points occurs in exactly v-1distinct 2393# blocks. Hence the rows and the columns of the incidence 2394# matrix M of the design are always of constant weight. 2395# Take a Frobenius (G,+) group with kernel K and complement H. 2396# Consider the design D with point set K and block set 2397# { a^H + b | a, b in K, a <> 0 }. 2398# Here a^H denotes the orbit of a under conjugation by elements 2399# of H. Every planar near-ring design of type "*" can be obtained 2400# in this way from groups. A group K together with a group of 2401# automorphism H of K such that the semidirect product KH is a 2402# Frobenius group with complement H is called a Ferrero pair (K, H) 2403# in SONATA. 2404# 2405# INPUT: P is a list of prime powers describing an abelian group G 2406# m > 0 is an integer such that G admits a cyclic fpf 2407# automorphism group of size m 2408# This means that for all q = p^k in P, OrderMod( p, m ) must divide q 2409# (see the SONATA documentation for FpfAutomorphismGroupsCyclic). 2410# OUTPUT: The binary linear code whose generator matrix is the 2411# incidence matrix of a design associated to a "Ferrero pair" arising 2412# from the fixed-point-free (fpf) automorphism group of G. 2413# The pair (H,K) is called a Ferraro pair and the semidirect product KH is a 2414# Frobenius group with complement H. 2415# AUTHORS: Peter Mayr and David Joyner 2416#local C, f, H, K, M, D; 2417# LoadPackage("sonata"); 2418# f := FpfAutomorphismGroupsCyclic( P, m ); 2419# K := f[2]; 2420# H := Group( f[1][1] ); 2421# D := DesignFromFerreroPair( K, H, "*" ); 2422# M := IncidenceMat( D ); 2423# C:=GeneratorMatCode(M*Z(2), GF(2)); 2424#return C; 2425#end); 2426 2427#Having trouble getting GUAVA to load without errors if 2428#SONATA is not installed. Uncomment this and reload if you 2429#have SONATA. 2430#FerreroDesignCode:=function( P,m) 2431#local C, f, H, K, M, D; 2432# LoadPackage("sonata"); 2433# f := FpfAutomorphismGroupsCyclic( P, m ); 2434# K := f[2]; 2435# H := Group( f[1][1] ); 2436# D := DesignFromFerreroPair( K, H, "*" ); 2437# M := IncidenceMat( D ); 2438# C:=GeneratorMatCode(M*Z(2), GF(2)); 2439#return C; 2440#end; 2441 2442############################################################################# 2443## 2444#F QuasiCyclicCode( <G>, <s>, <F> ) . . . . . . . . . . . quasi cyclic code 2445## 2446## QuasiCyclicCode ( <G>, <s>, <F> ) generates a rate 1/m quasi-cyclic 2447## codes. Note that <G> is a list of univariate polynomial and m is the 2448## cardinality of this list. The integer s is the size of the circulant 2449## and it is not necessarily equal to the code dimension, i.e. k <= s. 2450## The associated field is <F>. 2451## 2452InstallMethod(QuasiCyclicCode, "linear quasi-cyclic code", true, 2453 [IsList, IsInt, IsField], 0, 2454function( L1, s, F ) 2455 # 2456 # A rate 1/m quasi-cyclic code contains m circulant matrices, each of 2457 # the same size, and in this case they are all s x s circulant matrices. 2458 # Each circulant can be specified by a univariate polynomial. 2459 # 2460 local i, m, v, M, C; 2461 2462 # Determine the cardinality of the list L1 2463 m:=Size(L1); 2464 if (m < 2) then 2465 Error("The cardinality of <G> must be at least 2\n"); 2466 fi; 2467 2468 # Make sure that all the list elements are univariate polynomials 2469 for i in [1..m] do; 2470 if (IsUnivariatePolynomial(L1[i]) = false) then 2471 Error("All list elements must be univariate polynomials\n"); 2472 fi; 2473 if (Degree(L1[i]) >= s) then 2474 Error("The degree of the polynomial must be less than s\n"); 2475 fi; 2476 od; 2477 2478 # Convert each univariate polynomial into a circulant matrix and 2479 # concatenate them to generate a generator matrix 2480 M:=[]; 2481 for i in [1..m] do; 2482 v:=ShallowCopy( CoefficientsOfUnivariatePolynomial(L1[i]) ); 2483 Append( v, List([1..(s - (Degree(L1[i])+1))], i->Zero(F)) ); 2484 M:=Concatenation( M, CirculantMatrix(s, v) ); 2485 od; 2486 C := GeneratorMatCode( TransposedMat(M), F ); 2487 C!.name := "quasi-cyclic code"; 2488 return C; 2489end); 2490 2491InstallOtherMethod(QuasiCyclicCode, "binary linear quasi-cyclic code", 2492 true, [IsList, IsInt], 0, 2493function( L1, s ) 2494 2495 local i, j, m, a, t, v, L2, LUT; 2496 2497 LUT:=[ "000", "001", "010", "011", "100", "101", "110", "111" ]; 2498 2499 # Determine the cardinality of the list L1 2500 m := Size(L1); 2501 if (m < 2) then 2502 Error("The cardinality of <G> must be at least 2\n"); 2503 fi; 2504 2505 L2 := []; 2506 for i in [1..m] do; 2507 if (IsInt(L1[i]) = false) then 2508 Error("All list elements must be in octal\n"); 2509 fi; 2510 a := String(L1[i]); 2511 v := []; 2512 for j in [1..Length(a)] do; 2513 t := INT_CHAR(a[j]) - 48; # Conversion of ASCII character to integer 2514 if (t > 7) then 2515 Error("All list elements must be in octal\n"); 2516 fi; 2517 Append(v, LUT[t+1]); 2518 od; 2519 Append(L2, [ReciprocalPolynomial(PolyCodeword(Codeword(v)))]); 2520 od; 2521 2522 return QuasiCyclicCode( L2, s, GF(2) ); 2523end); 2524 2525##################################################################### 2526## 2527#F CyclicMDSCode( <q>, <m>, <k> ) . . . . . . . . . cyclic MDS code 2528## 2529## Construct a [q^m + 1, k, q^m - k + 2] cyclic MDS code over GF(q^m) 2530## 2531InstallMethod(CyclicMDSCode, "method for linear code", true, 2532 [IsInt, IsInt, IsInt], 0, 2533function(q, m, k) 2534 local i, j, l, x, a, r, g, G, F, CS, C, dmin; 2535 2536 if (k < 1) or (k > q^m + 1) then 2537 Error("Incorrect parameter, 1 <= k <= q^m+1.\n"); 2538 fi; 2539 if IsEvenInt(k) and IsOddInt(q^m) then 2540 Error("Cannot construct such code, k must be odd for odd field size.\n"); 2541 fi; 2542 2543 F := GF(q^m); 2544 x := Indeterminate(F, "x"); 2545 a := PrimitiveUnityRoot(F, q^m+1); # Primitive (q^m + 1)-st root of unity 2546 CS := CyclotomicCosets(q^m, q^m + 1); 2547 dmin := q^m - k + 2; 2548 2549 # R. Roth's book (Prob. 8.15, pp. 262) 2550 # If q^m is odd, there exists [q^m + 1, k, q^m - k + 2] cyclic MDS codes for 2551 # odd values of k in the range 1 <= k <= q^m. This is because the cyclotomic 2552 # cosets of q^m mod q^m + 1 are 2553 # { 0, [1,q^m], [2,q^m-2], ..., [(q^m-1)/2, (q^m+3)/2], (q^m+1)/2 }. 2554 # There are two single elements in the above cosets, { 0 } and { (q^m+1)/2 }. 2555 # If k is odd, dmin = q^m - k + 2 is even and delta = dmin-1 is odd (BCH bound). 2556 # This can be easily obtained by including either { 0 } or { (q^m+1)/2 }. 2557 # On the other hand, if k is even, dmin is odd and delta is even. With reference 2558 # to the above cyclotomic cosets, we cannot have exactly delta consecutive integers. 2559 # (Does this definitely mean this kind of code cannot be constructed??) 2560 # 2561 # If q^m is even, there exists [q^m + 1, k, q^m - k + 2] cyclic MDS codes for 2562 # any value of k in the range 1 <= k <= q^m+1. This is because the cyclotomic 2563 # cosets of q^m mod q^m + 1 are 2564 # { 0, [1,q^m], [2,q^m-2], ..., [q^m/2, 1 + q^m/2] }. 2565 # Consequently, we can easily construct a set of odd or even consecutive integers 2566 # using the cyclotomic cosets above. 2567 2568 if IsEvenInt(k) then 2569 g := (x + a^0); 2570 r := [ a^0 ]; 2571 for i in [1..((dmin-2)/2)] do; 2572 g := g * (x + a^(CS[i+1][1])) * (x + a^(CS[i+1][2])); 2573 r := Concatenation(r, [ a^(CS[i+1][1]), a^(CS[i+1][2]) ]); 2574 od; 2575 else 2576 if IsOddInt(q^m) then 2577 g := (x + a^(CS[Size(CS)][1])); 2578 r := [ a^(CS[Size(CS)][1]) ]; 2579 j := Size(CS)-1; 2580 l := (dmin-2)/2; 2581 else 2582 g := 1; 2583 r := []; 2584 j := Size(CS); 2585 l := (dmin-1)/2; 2586 fi; 2587 for i in [0..l-1] do; 2588 g := g * (x + a^(CS[j-i][1])) * (x + a^(CS[j-i][2])); 2589 r := Concatenation(r, [ a^(CS[j-i][1]), a^(CS[j-i][2]) ]); 2590 od; 2591 fi; 2592 2593 G := GeneratorMatrixFromPoly(g, q^m + 1); 2594 C := GeneratorMatCodeNC(G, F); 2595 2596 C!.name := "MDS code"; 2597 2598 # We know these bounds as it is an MDS code 2599 C!.lowerBoundMinimumDistance := dmin; 2600 C!.upperBoundMinimumDistance := dmin; 2601 2602 # Tell it that it is a cyclic code 2603 SetIsCyclicCode(C, true); 2604 SetGeneratorPol(C, g); 2605 2606 # Also tell it that it is an MDS code 2607 IsMDSCode(C); 2608 2609 return C; 2610end); 2611 2612####################################################################### 2613## 2614#F FourNegacirculantSelfDualCode( <ax>, <bx>, <k> ) . . self-dual code 2615## 2616## Construct a [2*k, k, d] self-dual code over F using Harada's 2617## construction. See: 2618## 2619## 1. M. Harada and T. Nishimura, "An extremal singly even self 2620## dual code of length 88", Advances in Mathematics of 2621## Communications, vol 1, no. 2, pp. 261--267, 2007 2622## 2623## 2. M. Harada, W. Holzmann, H. Kharaghani and M. Khorvash, 2624## "Extremal ternary self-dual codes constructed from 2625## negacirculant matrices", Graph and Combinatorics, vol 23, 2626## pp. 401--417, 2007 2627## 2628## 3. M. Harada, "An extremal doubly even self-dual code of 2629## length 112", preprint 2630## 2631## The generator matrix of the code has the following form: 2632## 2633## - - 2634## | : A : B | 2635## G = | I :------:-----| 2636## | : -B^T : A^T | 2637## - - 2638## 2639## Note that the matrices A, B, A^T and B^T are k/2 * k/2 2640## negacirculant matrices. 2641## 2642__G_FourNegacirculantSelfDualCode := function(ax, bx, k) 2643 local i, v, m, x, FA, FB, A, AT, B, BT, G; 2644 2645 if IsOddInt(k) then 2646 Error("k must be an even integer\n"); 2647 fi; 2648 2649 m := k/2; 2650 2651 # Determine field size 2652 FA := Field(VectorCodeword(Codeword(ax))); 2653 FB := Field(VectorCodeword(Codeword(bx))); 2654 if FA <> FB then 2655 Error("Polynomials a(x) and b(x) must have elements from the same field\n"); 2656 fi; 2657 2658 x := Indeterminate(FA); 2659 2660 v := MutableCopyMat( CoefficientsOfUnivariatePolynomial(ax) ); 2661 Append( v, List([1..(m - (Degree(ax)+1))], i->Zero(FA)) ); 2662 A := NegacirculantMatrix(m, v*One(FA)); 2663 AT:= TransposedMat(A); 2664 2665 v := MutableCopyMat( CoefficientsOfUnivariatePolynomial(bx) ); 2666 Append( v, List([1..(m - (Degree(bx)+1))], i->Zero(FA)) ); 2667 B := NegacirculantMatrix(m, v*One(FA)); 2668 BT:= TransposedMat(-B); 2669 2670 G := IdentityMat(k, One(FA)); 2671 2672 # [ A | B ] 2673 for i in [1..m] do; 2674 Append(G[i], A[i]); 2675 Append(G[i], B[i]); 2676 od; 2677 2678 # [ B^T | A^T ] 2679 for i in [1..m] do; 2680 Append(G[m+i], BT[i]); 2681 Append(G[m+i], AT[i]); 2682 od; 2683 2684 return G; 2685end; 2686 2687InstallMethod(FourNegacirculantSelfDualCode, "method for binary linear code", true, 2688 [IsUnivariatePolynomial, IsUnivariatePolynomial, IsInt], 0, 2689function(ax, bx, k) 2690 local G, C, F; 2691 2692 # Obtain the generator matrix 2693 G := __G_FourNegacirculantSelfDualCode(ax, bx, k); 2694 F := Field(VectorCodeword(Codeword(ax))); 2695 2696 C := GeneratorMatCode(G, "four-negacirculant self-dual code", F); 2697 C!.GeneratorMat := ShallowCopy(G); 2698 2699 if (IsSelfDualCode(C) = false) then 2700 Error("Polynomials a(x) and b(x) do not produce a self-dual code\n"); 2701 fi; 2702 2703 return C; 2704end); 2705 2706# Faster version - no minimum distance and covering radius estimation 2707InstallMethod(FourNegacirculantSelfDualCodeNC, "method for binary linear code", true, 2708 [IsUnivariatePolynomial, IsUnivariatePolynomial, IsInt], 0, 2709function(ax, bx, k) 2710 local G, C, F; 2711 2712 # Obtain the generator matrix 2713 G := __G_FourNegacirculantSelfDualCode(ax, bx, k); 2714 F := Field(VectorCodeword(Codeword(ax))); 2715 2716 C := GeneratorMatCodeNC(G, F); 2717 C!.name := "four-negacirculant self-dual code"; 2718 C!.GeneratorMat := ShallowCopy(G); 2719 C!.lowerBoundMinimumDistance := 1; 2720 C!.upperBoundMinimumDistance := k+1; 2721 C!.boundsCoveringRadius := [ 0, WordLength(C) ]; 2722 2723 if (IsSelfDualCode(C) = false) then 2724 Error("Polynomials a(x) and b(x) do not produce a self-dual code\n"); 2725 fi; 2726 2727 return C; 2728end); 2729 2730