/ erasure_code / share.go
share.go
1 package erasure_code 2 3 import "fmt" 4 5 func init() { 6 galoisInit() 7 } 8 9 // Finite fields 10 // ============= 11 12 type ZeroDivisionError struct { 13 } 14 15 func (e *ZeroDivisionError) Error() string { 16 return "division by zero" 17 } 18 19 // should panic with ZeroDivisionError when dividing by zero 20 type Field interface { 21 Add(b Field) Field 22 Sub(b Field) Field 23 Mul(b Field) Field 24 Div(b Field) Field 25 Value() int 26 Factory() FieldFactory 27 } 28 type FieldFactory interface { 29 Construct(v int) Field 30 } 31 32 // per-byte 2^8 Galois field 33 // Note that this imposes a hard limit that the number of extended chunks can 34 // be at most 256 along each dimension 35 type Galois struct { 36 v uint8 37 } 38 39 var gexptable [255]uint8 40 var glogtable [256]uint8 41 42 func galoisTpl(a uint8) uint8 { 43 r := a ^ (a << 1) // a * (x+1) 44 if (a & (1 << 7)) != 0 { 45 // would overflow (have an x^8 term); reduce by the AES polynomial, 46 // x^8 + x^4 + x^3 + x + 1 47 return r ^ 0x1b 48 } else { 49 return r 50 } 51 } 52 53 func galoisInit() { 54 var v uint8 = 1 55 for i := uint8(0); i < 255; i++ { 56 glogtable[v] = i 57 gexptable[i] = v 58 v = galoisTpl(v) 59 } 60 } 61 62 func (a *Galois) Add(_b Field) Field { 63 b := _b.(*Galois) 64 return &Galois{a.v ^ b.v} 65 } 66 func (a *Galois) Sub(_b Field) Field { 67 b := _b.(*Galois) 68 return &Galois{a.v ^ b.v} 69 } 70 func (a *Galois) Mul(_b Field) Field { 71 b := _b.(*Galois) 72 if a.v == 0 || b.v == 0 { 73 return &Galois{0} 74 } 75 return &Galois{gexptable[(int(glogtable[a.v])+ 76 int(glogtable[b.v]))%255]} 77 } 78 func (a *Galois) Div(_b Field) Field { 79 b := _b.(*Galois) 80 if b.v == 0 { 81 panic(ZeroDivisionError{}) 82 } 83 if a.v == 0 { 84 return &Galois{0} 85 } 86 return &Galois{gexptable[(int(glogtable[a.v])+255- 87 int(glogtable[b.v]))%255]} 88 } 89 func (a *Galois) Value() int { 90 return int(a.v) 91 } 92 func (a *Galois) String() string { 93 return fmt.Sprintf("%d",a.v) 94 } 95 96 type galoisFactory struct { 97 } 98 99 func GaloisFactory() FieldFactory { 100 return &galoisFactory{} 101 } 102 func (self *Galois) Factory() FieldFactory { 103 return GaloisFactory() 104 } 105 func (self *galoisFactory) Construct(v int) Field { 106 return &Galois{uint8(v)} 107 } 108 109 // Modular arithmetic class 110 type modulo struct { 111 v uint 112 n uint // the modulus 113 } 114 115 func (a *modulo) Add(_b Field) Field { 116 b := _b.(*modulo) 117 return &modulo{(a.v + b.v) % a.n, a.n} 118 } 119 120 func (a *modulo) Sub(_b Field) Field { 121 b := _b.(*modulo) 122 return &modulo{(a.v + a.n - b.v) % a.n, a.n} 123 } 124 125 func (a *modulo) Mul(_b Field) Field { 126 b := _b.(*modulo) 127 return &modulo{(a.v * b.v) % a.n, a.n} 128 } 129 130 func powmod(b uint, e uint, m uint) uint { 131 var r uint = 1 132 for e > 0 { 133 if (e & 1) == 1 { 134 r = (r * b) % m 135 } 136 b = (b * b) % m 137 e >>= 1 138 } 139 return r 140 } 141 142 func (a *modulo) Div(_b Field) Field { 143 b := _b.(*modulo) 144 return &modulo{(a.v * powmod(b.v, a.n-2, a.n)) % a.n, a.n} 145 } 146 147 func (self *modulo) Value() int { 148 return int(self.v) 149 } 150 151 type moduloFactory struct { 152 n uint 153 } 154 155 func (self *modulo) Factory() FieldFactory { 156 return &moduloFactory{self.n} 157 } 158 159 func (self *moduloFactory) Construct(v int) Field { 160 return &modulo{uint(v), self.n} 161 } 162 163 func MakeModuloFactory(n uint) FieldFactory { 164 return &moduloFactory{n} 165 } 166 167 func zero(f FieldFactory) Field { 168 return f.Construct(0) 169 } 170 func one(f FieldFactory) Field { 171 return f.Construct(1) 172 } 173 174 // Helper functions 175 // ================ 176 177 // Evaluates a polynomial p in little-endian form (e.g. x^2 + 3x + 2 is 178 // represented as [2, 3, 1]) at coordinate x, 179 func EvalPolyAt(poly []Field, x Field) Field { 180 arithmetic := x.Factory() 181 r, xi := zero(arithmetic), one(arithmetic) 182 for _, ci := range poly { 183 r = r.Add(xi.Mul(ci)) 184 xi = xi.Mul(x) 185 } 186 return r 187 } 188 189 // Given p+1 y values and x values with no errors, recovers the original 190 // p+1 degree polynomial. For example, 191 // LagrangeInterp({51.0, 59.0, 66.0}, {1, 3, 4}) = {50.0, 0, 1.0} 192 // (or it would be, if floats were Fields) 193 func LagrangeInterp(pieces []Field, xs []Field) []Field { 194 arithmetic := pieces[0].Factory() 195 zero, one := zero(arithmetic), one(arithmetic) 196 197 // `size` is the number of datapoints; the degree of the result polynomial 198 // is then `size-1` 199 size := len(pieces) 200 201 root := []Field{one} // initially just the polynomial "1" 202 // build up the numerator polynomial, `root`, by taking the product of (x-v) 203 // (implemented as convolving repeatedly with [-v, 1]) 204 for _, v := range xs { 205 // iterate backward since new root[i] depends on old root[i-1] 206 for i := len(root) - 1; i >= 0; i-- { 207 root[i] = root[i].Mul(zero.Sub(v)) 208 if i > 0 { 209 root[i] = root[i].Add(root[i-1]) 210 } 211 } 212 // polynomial is always monic so save an extra multiply by doing this 213 // after 214 root = append(root, one) 215 } 216 217 // generate per-value numerator polynomials by dividing the master 218 // polynomial back by each x coordinate 219 nums := make([][]Field, size) 220 for i, v := range xs { 221 // divide `root` by (x-v) to get a degree size-1 polynomial 222 // (i.e. with `size` coefficients) 223 num := make([]Field, size) 224 // compute the x^0, x^1, ..., x^(p-2) coefficients by long division 225 last := one 226 num[len(num)-1] = last // still always a monic polynomial 227 for j := size - 2; j >= 0; j-- { 228 last = root[j+1].Add(last.Mul(v)) 229 num[j] = last 230 } 231 nums[i] = num 232 } 233 234 // generate denominators by evaluating numerator polys at their x 235 denoms := make([]Field, size) 236 for i, x := range xs { 237 denoms[i] = EvalPolyAt(nums[i], x) 238 } 239 240 // generate output polynomial by taking the sum over i of 241 // (nums[i] * pieces[i] / denoms[i]) 242 sum := make([]Field, size) 243 for i := range sum { 244 sum[i] = zero 245 } 246 for i, y := range pieces { 247 factor := y.Div(denoms[i]) 248 // add nums[i] * factor to sum, as a vector 249 for j := 0; j < size; j++ { 250 sum[j] = sum[j].Add(nums[i][j].Mul(factor)) 251 } 252 } 253 return sum 254 } 255 256 // Given two linear equations, eliminates the first variable and returns 257 // the resulting equation. 258 // 259 // An equation of the form a_1 x_1 + ... + a_n x_n + b = 0 260 // is represented as the array [a_1, ..., a_n, b]. 261 func elim(a []Field, b []Field) []Field { 262 result := make([]Field, len(a)-1) 263 for i := range result { 264 result[i] = a[i+1].Mul(b[0]).Sub(b[i+1].Mul(a[0])) 265 } 266 return result 267 } 268 269 // Given one homogeneous linear equation and the values of all but the first 270 // variable, solve for the value of the first variable. 271 // 272 // For an equation of the form 273 // a_1 x_1 + ... + a_n x_n = 0 274 // pass two arrays, [a_1, ..., a_n] and [x_2, ..., x_n]. 275 func evaluate(coeffs []Field, vals []Field) Field { 276 total := zero(coeffs[0].Factory()) 277 for i, val := range vals { 278 total = total.Sub(coeffs[i+1].Mul(val)) 279 } 280 return total.Div(coeffs[0]) 281 } 282 283 // Given an n*n system of inhomogeneous linear equations, solve for the value of 284 // every variable. 285 // 286 // For equations of the form 287 // a_1,1 x_1 + ... + a_1,n x_n + b_1 = 0 288 // a_2,1 x_1 + ... + a_2,n x_n + b_2 = 0 289 // ... 290 // a_n,1 x_1 + ... + a_n,n x_n + b_n = 0 291 // pass a two-dimensional array 292 // [[a_1,1, ..., a_1,n, b_1], ..., [a_n,1, ..., a_n,n, b_n]]. 293 // 294 // Returns the values of [x_1, ..., x_n]. 295 func SysSolve(eqs [][]Field) []Field { 296 arithmetic := eqs[0][0].Factory() 297 backEqs := make([][]Field, 1, len(eqs)) 298 backEqs[0] = eqs[0] 299 300 for len(eqs) > 1 { 301 neweqs := make([][]Field, len(eqs)-1) 302 for i := 0; i < len(eqs)-1; i++ { 303 neweqs[i] = elim(eqs[i], eqs[i+1]) 304 } 305 eqs = neweqs 306 // find a row with a nonzero first entry 307 i := 0 308 for i+1 < len(eqs) && eqs[i][0].Value() == 0 { 309 i++ 310 } 311 backEqs = append(backEqs, eqs[i]) 312 } 313 314 kvals := make([]Field, len(backEqs)+1) 315 kvals[len(backEqs)] = one(arithmetic) 316 // back-substitute in reverse order 317 // (smallest to largest equation) 318 for i := len(backEqs) - 1; i >= 0; i-- { 319 kvals[i] = evaluate(backEqs[i], kvals[i+1:]) 320 } 321 322 return kvals[:len(kvals)-1] 323 } 324 325 // Divide two polynomials with nonzero leading terms. 326 // T should be a field. 327 func PolyDiv(Q []Field, E []Field) []Field { 328 if len(Q) < len(E) { 329 return []Field{} 330 } 331 div := make([]Field, len(Q)-len(E)+1) 332 for i := len(div) - 1; i >= 0; i-- { 333 factor := Q[len(Q)-1].Div(E[len(E)-1]) 334 div[i] = factor 335 // subtract factor * E * x^i from Q 336 Q = Q[:len(Q)-1] // the highest term should cancel 337 for j := 0; j < len(E)-1; j++ { 338 Q[i+j] = Q[i+j].Sub(factor.Mul(E[j])) 339 } 340 } 341 return div 342 } 343 344 func trySysSolve(eqs [][]Field) (ret []Field, ok bool) { 345 defer func() { 346 err := recover() 347 if err == nil { 348 return 349 } 350 switch err := err.(type) { 351 case ZeroDivisionError: 352 ret = nil 353 ok = false 354 default: 355 panic(err) 356 } 357 }() 358 return SysSolve(eqs), true 359 } 360 361 type NotEnoughData struct { } 362 func (self *NotEnoughData) Error() string { 363 return "Not enough data!" 364 } 365 type TooManyErrors struct { } 366 func (self *TooManyErrors) Error() string { 367 return "Answer doesn't match (too many errors)!" 368 } 369 370 // Given a set of y coordinates and x coordinates, and the degree of the 371 // original polynomial, determines the original polynomial even if some of the y 372 // coordinates are wrong. If m is the minimal number of pieces (ie. degree + 373 // 1), t is the total number of pieces provided, then the algo can handle up to 374 // (t-m)/2 errors. 375 func BerlekampWelchAttempt(pieces []Field, xs []Field, masterDegree int) ([]Field, error) { 376 errorLocatorDegree := (len(pieces) - masterDegree - 1) / 2 377 arithmetic := pieces[0].Factory() 378 zero, one := zero(arithmetic), one(arithmetic) 379 // Set up the equations for y[i]E(x[i]) = Q(x[i]) 380 // degree(E) = errorLocatorDegree 381 // degree(Q) = masterDegree + errorLocatorDegree - 1 382 eqs := make([][]Field, 2*errorLocatorDegree+masterDegree+1) 383 for i := range eqs { 384 eq := []Field{} 385 x := xs[i] 386 piece := pieces[i] 387 neg_x_j := zero.Sub(one) 388 for j := 0; j < errorLocatorDegree+masterDegree+1; j++ { 389 eq = append(eq, neg_x_j) 390 neg_x_j = neg_x_j.Mul(x) 391 } 392 x_j := one 393 for j := 0; j < errorLocatorDegree+1; j++ { 394 eq = append(eq, x_j.Mul(piece)) 395 x_j = x_j.Mul(x) 396 } 397 eqs[i] = eq 398 } 399 // Solve the equations 400 // Assume the top error polynomial term to be one 401 errors := errorLocatorDegree 402 ones := 1 403 var polys []Field 404 for errors >= 0 { 405 if p, ok := trySysSolve(eqs); ok { 406 for i := 0; i < ones; i++ { 407 p = append(p, one) 408 } 409 polys = p 410 break 411 } 412 // caught ZeroDivisionError 413 eqs = eqs[:len(eqs)-1] 414 for i, eq := range eqs { 415 eq[len(eq)-2] = eq[len(eq)-2].Add(eq[len(eq)-1]) 416 eqs[i] = eq[:len(eq)-1] 417 } 418 errors-- 419 ones++ 420 } 421 if errors < 0 { 422 return nil, &NotEnoughData{} 423 } 424 // divide the polynomials... 425 split := errorLocatorDegree + masterDegree + 1 426 div := PolyDiv(polys[:split], polys[split:]) 427 corrects := 0 428 for i := 0; i < len(xs); i++ { 429 if EvalPolyAt(div, xs[i]).Value() == pieces[i].Value() { 430 corrects++ 431 } 432 } 433 if corrects < masterDegree+errors { 434 return nil, &TooManyErrors{} 435 } 436 return div, nil 437 } 438 439 // Extends a list of integers in [0 ... 255] (if using Galois arithmetic) by 440 // adding n redundant error-correction values 441 func Extend(data []int, n int, arithmetic FieldFactory) ([]int, error) { 442 size := len(data) 443 444 dataF := make([]Field, size) 445 for i, d := range data { 446 dataF[i] = arithmetic.Construct(d) 447 } 448 449 xs := make([]Field, size) 450 for i := range xs { 451 xs[i] = arithmetic.Construct(i) 452 } 453 454 poly, err := BerlekampWelchAttempt(dataF, xs, size-1) 455 if err != nil { 456 return nil, err 457 } 458 459 for i := 0; i < n; i++ { 460 data = append(data, 461 int(EvalPolyAt(poly, arithmetic.Construct(size+i)).Value())) 462 } 463 return data, nil 464 } 465 466 // Repairs a list of integers in [0 ... 255]. Some integers can be erroneous, 467 // and you can put -1 in place of an integer if you know that a certain 468 // value is defective or missing. Uses the Berlekamp-Welch algorithm to 469 // do error-correction 470 func Repair(data []int, datasize int, arithmetic FieldFactory) ([]int, error) { 471 vs := make([]Field, 0, len(data)) 472 xs := make([]Field, 0, len(data)) 473 for i, d := range data { 474 if d >= 0 { 475 vs = append(vs, arithmetic.Construct(d)) 476 xs = append(xs, arithmetic.Construct(i)) 477 } 478 } 479 480 poly, err := BerlekampWelchAttempt(vs, xs, datasize-1) 481 if err != nil { 482 return nil, err 483 } 484 485 result := make([]int, len(data)) 486 for i := range result { 487 result[i] = int(EvalPolyAt(poly, arithmetic.Construct(i)).Value()) 488 } 489 return result, nil 490 } 491 492 func transpose(d [][]int) [][]int { 493 width := len(d[0]) 494 result := make([][]int, width) 495 for i := range result { 496 col := make([]int, len(d)) 497 for j := range col { 498 col[j] = d[j][i] 499 } 500 result[i] = col 501 } 502 return result 503 } 504 505 func extractColumn(d [][]int, j int) []int { 506 result := make([]int, len(d)) 507 for i, row := range d { 508 result[i] = row[j] 509 } 510 return result 511 } 512 513 // Extends a list of bytearrays 514 // eg. ExtendChunks([map(ord, 'hello'), map(ord, 'world')], 2) 515 // n is the number of redundant error-correction chunks to add 516 func ExtendChunks(data [][]int, n int, arithmetic FieldFactory) ([][]int, error) { 517 width := len(data[0]) 518 o := make([][]int, width) 519 for i := 0; i < width; i++ { 520 row, err := Extend(extractColumn(data, i), n, arithmetic) 521 if err != nil { 522 return nil, err 523 } 524 o[i] = row 525 } 526 return transpose(o), nil 527 } 528 529 // Repairs a list of bytearrays. Use an empty array in place of a missing array. 530 // Individual arrays can contain some missing or erroneous data. 531 func RepairChunks(data [][]int, datasize int, arithmetic FieldFactory) ([][]int, error) { 532 var width int 533 for _, row := range data { 534 if len(row) > 0 { 535 width = len(row) 536 break 537 } 538 } 539 filledData := make([][]int, len(data)) 540 for i, row := range data { 541 if len(row) == 0 { 542 filledData[i] = make([]int, width) 543 for j := range filledData[i] { 544 filledData[i][j] = -1 545 } 546 } else { 547 filledData[i] = row 548 } 549 } 550 o := make([][]int, width) 551 for i := range o { 552 row, err := Repair(extractColumn(data, i), datasize, arithmetic) 553 if err != nil { 554 return nil, err 555 } 556 o[i] = row 557 } 558 return transpose(o), nil 559 }