///|
/// Exact affine erasure analysis. Only bits shared by every completion are marked
/// recovered; a particular valid completion is separate and is NOT evidence of uniqueness.
pub struct RecoveryReport {
  values : Array[Int32?]
  recovered : Int32
  unresolved : Array[Int32]
  ambiguity_dimension : Int32
  unique : Bool
  particular_word : Array[Int32]
} derive(Debug)

///|
/// Gaussian fallback for BEC stopping sets. This uses dense GF(2), unlike linear-time
/// peeling, and is limited to 2048 unknowns/checks and one million augmented cells.
pub fn ParityCheck::recover_erasures(
  self : ParityCheck,
  observations : Array[Int32?],
) -> RecoveryReport raise CodeError {
  if observations.length() != self.n {
    raise LengthMismatch(self.n, observations.length())
  }
  let unknown_positions : Array[Int32] = []
  let unknown_index = Array::make(self.n, -1)
  let particular_word = Array::make(self.n, 0)
  for v, observation in observations {
    match observation {
      None => {
        unknown_index[v] = unknown_positions.length()
        unknown_positions.push(v)
      }
      Some(bit) => {
        if bit != 0 && bit != 1 {
          raise InvalidBit(v, bit)
        }
        particular_word[v] = bit
      }
    }
  }
  let unknown_count = unknown_positions.length()
  if unknown_count > 2048 ||
    self.checks() > 2048 ||
    (unknown_count + 1).to_int64() * self.checks().to_int64() > 1048576L {
    raise ResourceLimit("erasure elimination dimensions")
  }
  let rows = Array::makei(self.checks(), fn(c) {
    let row = Array::make(unknown_count + 1, 0)
    let mut rhs = 0
    for v in self.rows[c] {
      let index = unknown_index[v]
      if index >= 0 {
        row[index] = 1
      } else {
        rhs = rhs ^ particular_word[v]
      }
    }
    row[unknown_count] = rhs
    row
  })
  let pivots : Array[Int32] = []
  let mut rank = 0
  let mut work = 0L
  for col in 0.. 50000000L {
          raise ResourceLimit("erasure elimination work")
        }
        for k in col..<(unknown_count + 1) {
          rows[r][k] = rows[r][k] ^ rows[rank][k]
        }
      }
    }
    pivots.push(col)
    rank += 1
    if rank == rows.length() {
      break
    }
  }
  for r in rank.. ParityCheck raise CodeError {
  if self.n > 2048 ||
    self.checks() > 2048 ||
    self.n.to_int64() * self.checks().to_int64() > 1048576L {
    raise ResourceLimit("puncture elimination dimensions")
  }
  let deleted = Array::make(self.n, false)
  for v in removed {
    if v < 0 || v >= self.n {
      raise InvalidIndex(-1, v)
    }
    if deleted[v] {
      raise DuplicateEntry(-1, v)
    }
    deleted[v] = true
  }
  if removed.length() == self.n {
    raise InvalidDimensions("puncturing must retain at least one coordinate")
  }
  let dense = Array::makei(self.checks(), fn(c) {
    let row = Array::make(self.n, 0)
    for v in self.rows[c] {
      row[v] = 1
    }
    row
  })
  let mut rank = 0
  let mut work = 0L
  // Eliminate ONLY removed columns. Pivot rows constrain discarded coordinates
  // and vanish under projection; remaining rows are the constraints on retained bits.
  for col in 0.. 50000000L {
          raise ResourceLimit("puncture elimination work")
        }
        // Retained columns before this pivot also matter, so xor the full row.
        for k in 0.. Json::null()
      Some(bit) => Json(bit)
    }
  })
  {
    "values": Json(values),
    "recovered": self.recovered,
    "unresolved": Json(self.unresolved),
    "ambiguity_dimension": self.ambiguity_dimension,
    "unique": self.unique,
    "particular_word": Json(self.particular_word),
  }
}