`nx_fec_gf256.nx` source

↩ module page · 173 lines · 5871 B
1// nx_fec_gf256.nx -- the FEC CEILING: Reed-Solomon erasure coding over
2// GF(256).  MDS (maximum-distance-separable): an RS(n=K+R, K) code
3// recovers ANY R lost symbols, in ANY positions -- not just the specific
4// patterns XOR parity handles.  This is the DSN/CCSDS/RaptorQ-class
5// rateless recovery ceiling.  Sovereign, integer-only: GF(256) tables are
6// BUILT AT RUNTIME from the primitive polynomial 0x11D (sidesteps the
7// const-array-literal gap, same trick as nx_crc32c), and a CAUCHY parity
8// matrix guarantees every K x K submatrix is invertible (true MDS).
9//
10// Encoding matrix = [ I_K ; C ] where C[j][i] = 1 / ((K+j) ^ i) in GF(256)
11// (disjoint Cauchy point sets -> non-singular).  Decode: take the K rows
12// for any K surviving symbols -> K x K matrix M -> data = M^-1 * survivors
13// (Gauss-Jordan over GF(256)).  Reconstructs every erased symbol.
14//
15// license_tier: ORIGINAL
16
17import "nx_syscalls.nx"
18
19const GF_POLY: i64 = 0x11D     // x^8+x^4+x^3+x^2+1, generator 2
20
21// ---- GF(256) tables (exp sized 512 so multiply needs no mod) ----
22func gf_build(exp: *i64, log: *i64) -> i64 {
  var x: i64 = 1
  var i: i64 = 0
  while i < 255 {
      exp[i] = x
      log[x] = i
      x = x << 1
      if (x & 0x100) != 0 { x = x ^ GF_POLY }
      i = i + 1
  }
  var j: i64 = 255
  while j < 512 { exp[j] = exp[j - 255]; j = j + 1 }
  return 0
35}
36func gf_mul(exp: *i64, log: *i64, a: i64, b: i64) -> i64 {
  if a == 0 { return 0 }
  if b == 0 { return 0 }
  return exp[log[a] + log[b]]
40}
41func gf_inv(exp: *i64, log: *i64, a: i64) -> i64 { return exp[255 - log[a]] }
42
43func rs_cauchy(exp: *i64, log: *i64, j: i64, i: i64, K: i64) -> i64 {
  return gf_inv(exp, log, (K + j) ^ i)
45}
46
47// ---- encode: K data packets (S bytes each) -> R parity packets ----
48func rs_encode(exp: *i64, log: *i64, data: *u8, K: i64, R: i64, S: i64, parity: *u8) -> i64 {
  var jj: i64 = 0
  while jj < R {
      var b: i64 = 0
      while b < S {
          var acc: i64 = 0
          var ii: i64 = 0
          while ii < K {
              let c: i64 = rs_cauchy(exp, log, jj, ii, K)
              let dv: i64 = (data[ii * S + b] as i64) & 0xff
              acc = acc ^ gf_mul(exp, log, c, dv)
              ii = ii + 1
          }
          parity[jj * S + b] = (acc & 0xff) as u8
          b = b + 1
      }
      jj = jj + 1
  }
  return 0
67}
68
69// ---- Gauss-Jordan invert a K x K GF(256) matrix.  0 ok, -1 singular ----
70func gf_mat_invert(exp: *i64, log: *i64, M: *i64, K: i64, Minv: *i64) -> i64 {
  let W: i64 = 2 * K
  let A: *i64 = sys_mmap(K * W * 8 + 16) as *i64
  var r: i64 = 0
  while r < K {
      var c: i64 = 0
      while c < K { A[r * W + c] = M[r * K + c]; c = c + 1 }
      c = 0
      while c < K { if c == r { A[r * W + K + c] = 1 } else { A[r * W + K + c] = 0 } c = c + 1 }
      r = r + 1
  }
  var p: i64 = 0
  while p < K {
      if A[p * W + p] == 0 {
          var sr: i64 = p + 1
          var found: i64 = 0 - 1
          while sr < K { if A[sr * W + p] != 0 { found = sr; sr = K } else { sr = sr + 1 } }
          if found == (0 - 1) { return 0 - 1 }
          var c2: i64 = 0
          while c2 < W {
              let tmp: i64 = A[p * W + c2]
              A[p * W + c2] = A[found * W + c2]
              A[found * W + c2] = tmp
              c2 = c2 + 1
          }
      }
      let pivinv: i64 = gf_inv(exp, log, A[p * W + p])
      var c3: i64 = 0
      while c3 < W { A[p * W + c3] = gf_mul(exp, log, A[p * W + c3], pivinv); c3 = c3 + 1 }
      var rr: i64 = 0
      while rr < K {
          if rr != p {
              let factor: i64 = A[rr * W + p]
              if factor != 0 {
                  var c4: i64 = 0
                  while c4 < W {
                      A[rr * W + c4] = A[rr * W + c4] ^ gf_mul(exp, log, factor, A[p * W + c4])
                      c4 = c4 + 1
                  }
              }
          }
          rr = rr + 1
      }
      p = p + 1
  }
  r = 0
  while r < K {
      var c: i64 = 0
      while c < K { Minv[r * K + c] = A[r * W + K + c]; c = c + 1 }
      r = r + 1
  }
  return 0
122}
123
124// ---- decode from erasures.  recv = n=K+R symbols (data 0..K-1, parity
125// K..K+R-1) each S bytes; present[] flags survivors.  Reconstructs all K
126// data packets into out.  Returns 0, or -1 if < K survived. ----
127func rs_decode_erasures(exp: *i64, log: *i64, recv: *u8, present: *i64,
                      K: i64, R: i64, S: i64, out: *u8) -> i64 {
  let n: i64 = K + R
  let surv: *i64 = sys_mmap(K * 8 + 16) as *i64
  var cnt: i64 = 0
  var s: i64 = 0
  while s < n {
      if present[s] == 1 { if cnt < K { surv[cnt] = s; cnt = cnt + 1 } }
      s = s + 1
  }
  if cnt < K { return 0 - 1 }
  let M: *i64 = sys_mmap(K * K * 8 + 16) as *i64
  var r: i64 = 0
  while r < K {
      let sym: i64 = surv[r]
      var i: i64 = 0
      while i < K {
          if sym < K {
              if i == sym { M[r * K + i] = 1 } else { M[r * K + i] = 0 }
          } else {
              M[r * K + i] = rs_cauchy(exp, log, sym - K, i, K)
          }
          i = i + 1
      }
      r = r + 1
  }
  let Minv: *i64 = sys_mmap(K * K * 8 + 16) as *i64
  if gf_mat_invert(exp, log, M, K, Minv) != 0 { return 0 - 1 }
  var di: i64 = 0
  while di < K {
      var b: i64 = 0
      while b < S {
          var acc: i64 = 0
          var rr: i64 = 0
          while rr < K {
              let mv: i64 = Minv[di * K + rr]
              let val: i64 = (recv[surv[rr] * S + b] as i64) & 0xff
              acc = acc ^ gf_mul(exp, log, mv, val)
              rr = rr + 1
          }
          out[di * S + b] = (acc & 0xff) as u8
          b = b + 1
      }
      di = di + 1
  }
  return 0
173}
nx_fec_gf256.nx source

`nx_fec_gf256.nx` source
