Reed–Solomon codes: the polynomial that survives a scratch

1. Introduction

Have you ever noticed how a badly scratched CD still plays, or how a damaged QR code still scans?^[𝓳]i used to think the laser was just "really good at guessing", turns out it's a polynomial doing the job 😂 That's not luck or a better laser — it's a piece of 1960s algebra quietly rebuilding the destroyed bits. The trick is called a Reed-Solomon codes (named after two of their inventors), and the whole idea is almost simple^[𝓳]not that simple tbh if we learn from the first principles 😂: turn your data into a polynomial, then send more points of that polynomial (redundant) than we actually need.

The Error-correcting codes are a major object of study in mathematics, computer science, and electrical engineering; they belong to a field known as "Information Theory", one of the core areas of computer science.

In the Hamming codes blog, the enemy was a single flipped bit, and the question was where. Here the enemy is bigger - whole chunks of data go missing or arrive wrong, and the answer comes from polynomials over a finite field, the same Galois-field machinery I gestured at last time. These are my notes on how it all fits together, built up from first principles.

2. The problem

Picture me sending my friend four packets, each carrying one number. In a perfect world all four arrive. In the real world, the channel is noisy and things go wrong in two distinct flavours:

• an erasure — a packet goes missing and you know which one (the slot is empty);
• an error — a packet arrives corrupted and you don't know which one is lying.

Erasures are easy to tackle with: position known, value unknown. Errors are nastier: position and value unknown. We'll crack erasures first, then come back for errors.

3. First principles I — arithmetic on a clock

Before any of this works, we need a place to do the arithmetic where numbers can't run off to infinity — id est, a finite field. That brings us to the integers modulo a prime p. Two numbers are the same if they leave the same remainder when divided by p:

Think of it as a clock with p evenly-spaced numbers. An ordinary clock only has 1 → 2 → 3 → … → 12. What happens when it's 11 o'clock and two hours pass? Normal arithmetic insists:

But a clock doesn't have a 13. So what do we do? We wrap around: 13 → 1. Our number system pulls the same trick — keep only finitely many numbers, 0, 1, 2, 3, 4, and wrap. So what should 4 + 3 be here? It can't be 7, because 7 doesn't exist on this dial — that's exactly where mod steps in.

Addition is "walk clockwise"; the result always lands back on the dial. For example, starting at 4 and stepping 3 forward in mod 5 lands you on 2:

Subtraction is just adding the additive inverse — the number that brings you back to zero. Since 3 + 2 ≡ 0 (mod 5), the inverse of 3 is 2, so 1 − 3 ≡ 1 + 2 ≡ 3. Multiplication is repeated addition (2 × 3 = 6 ≡ 1), and division is multiplication by the multiplicative inverse:

To divide 4 by 2 mod 5: the inverse of 2 is 3 (because 2 × 3 ≡ 1), so 4 ÷ 2 ≡ 4 × 3 ≡ 12 ≡ 2. Here's the catch that makes p prime matter — if the modulus isn't prime, some numbers have no inverse at all (2 has none mod 4), and division breaks.

4. First principles II — polynomials as functions

The parabola you learned about in school does not apply here; a polynomial is simply a function where you put in an x to get out a y. Since there are only finitely many inputs in the finite field, this "function" is just a few points on the graph. Here's $f\left(x\right)=x^2+2$, first over the real line, then over GF(5):

There are two key facts concerning polynomials that form the foundation of what will follow:

1. A degree-d polynomial can have no more than d roots.
2. A unique degree d polynomial can be defined using d+1 points. Line can be defined by two points, parabola by three, etc.

Second point obviously follows from the first: two different degree-d polynomials P and Q both thread going through the same d+1 points. Let subtract one from the other:

Since P and Q pass through the same d+1 points, they have the same value at each of those points. Therefore, subtraction of equal quantities will always result in zero, thus R equals 0 at all d+1 points. In other words, R has d+1 roots.

However, R(x) remaind degree d or lower polynomial, while we already established that such polynomial cannot have more than d roots. This situation is impossible.

Only one exception is possible in this case — when polynomial is identically zero and equals zero for every x. Thus,R(x), which is 0 for every value of x. So R s the zero polynomial, which means that P(x) - Q(x) = 0, i.e., P(x) = Q(x).

Hence, the two polynomials are actually the same. It is proved that there exists one and only one unique degree-d polynomial that passes through arbitrary d + 1 points.

5. Lagrange interpolation

Unique it may be, but at the same time, the actual construction of the polynomial using the points is necessary. Let us start with the only element required for the entire process – Lagrange basis polynomial. It is a polynomial constructed in such a way that it functions as an on/off switch; it takes value of 1 at its own point and 0 everywhere else.

Building the switch takes two moves. Take three points: (1, 2), (3, 2), (4, -1), and let's make the switch for the first one. Move one — force the zeros: a polynomial is 0 exactly at its roots, so put roots at the other x-values, 3 and 4, giving $(x-3)(x-4)$. Move two — force the 1: at our own point x = 1 this evaluates to (1−3)(1−4) = 6 which is not 1, so divide the whole thing by 6:

Check it: $L_1\left(1\right)=1$, $L_1\left(3\right)=0$, $L_1\left(4\right)=0$. A perfect switch: on at its own point, off at the rest. Each point gets its own $L_i$ the same way.

Now Lagrange interpolation is almost an anticlimax: multiply each switch by the y-value you want at that point, and add them up...

The big product inside is just $L_i\left(x\right)$ written out. There is only one switch in action at any sample point, hence the sum goes through all of them. And this is where comes the magic: this technique uses solely +, −, × and ÷ operations and thus is valid over a finite field as well; in such a case, one deals with Lagrange interpolation modulo p, same algebra but on a clock.

Enough symbols — here is the actual thing to play with. Move the dots around and see how the unique curve passing through all of them changes dynamically in front of your eyes in real time. Then hit show the pieces to see how the polynomial decomposes into the linear combination of basis polynomials $y_i·L_i\left(x\right)$ — each colored curve is 1 at its own dot and 0 at all others and lifts the point. Add or remove points and the degree follows.

Tip: click on any point-chip on the right (or a dot) to select its basis function and have $L_i\left(x\right)$ be displayed with values substituted – just as you would manually do.

6. Why GF(2⁸)? — the field made of bytes

Up until now, everything had been working off a clock with a prime number of hours. However, data does not come in packages labeled "number mod 5" — data comes in packages called bytes. A Byte = 8 bits = 2⁸ possibilities = 256 possibilities = 0 to 255. There are no 257 possibilities. The wish list writes itself: We need a universe with 256 symbols (because bytes), and every symbol except 0 must have an inverse (because division through interpolation). Ordinary modulo cannot give us both. Watch it fail twice.

Attempt one — mod 256. Not prime (256 = 2⁸), and it shows. Try to divide by 2, i.e. hunt for its inverse:

The left side is always even and the right side is odd — no solution, ever. Division by 2 doesn't exist mod 256, and Lagrange breaks. Attempt two — mod 257. 257 is prime, so the arithmetic works beautifully… but computers don't store numbers in base 257. A byte naturally represents 0 … 255 and nothing else. Both doors are shut.

The mathematician's escape is a genuine plot twist: stop treating a byte as a number. Treat it as a polynomial — each bit is a coefficient, and the coefficients live in GF(2) = {0, 1}:

Addition costs nothing. In GF(2), 1 + 1 = 0 — so adding two polynomials just cancels matching terms. Add $(x^3+x+1)+(x^2+x+1)$: the x's kill each other, the 1's kill each other, and you're left with $x^3+x^2$. Now do the same in byte form: 1011 ⊕ 0111 = 1100 — the same answer. Adding polynomials over GF(2) is XOR, the single cheapest instruction a CPU owns. (And since 1 + 1 = 0 means −1 = 1, subtraction is the same operation. One instruction, two jobs.)

Multiplication is school algebra with one twist. Multiply normally, then reduce coefficients mod 2:

— the middle term vanished because 2 ≡ 0 (mod 2). Keep multiplying and everything stays comfortable… until $x^7·x=x^8$. A byte's 8 bits carry x⁰ through x⁷ — there is no bit for x⁸. It's the same wall as before: 200 + 100 = 300, and a byte can't store 300.

The first idea anyone tries is the lazy one: declare x⁸ = 1 and wrap around. After a lot of testing, mathematicians found it breaks the algebra: the result isn't a field, and some elements still end up with no inverse. Back to square one — unless we steal the trick that made the clock work. Mod a prime gave integers a field; what's the polynomial equivalent of a prime?

An irreducible polynomial. Primes refuse to break: 13 = a × b has no answer with both factors above 1, while 15 = 3 × 5 splits. Polynomials can split too — $x^2+3x+2=(x+1)(x+2)$ — but some refuse. Take x² + x + 1 over GF(2): plug in x = 0 and you get 1; plug in x = 1 and you get 1 + 1 + 1 = 1. Never zero, so it has no factors — irreducible. Primes are the atoms of the integers; irreducible polynomials are the atoms of polynomial arithmetic.

So instead of mod 256, work modulo a degree-8 irreducible polynomial. The one the world settled on:

Dividing by a degree-8 polynomial always leaves a remainder of degree at most 7 — which is exactly a byte. And the reduction step is friendlier than it looks. Working mod M(x) means M(x) ≡ 0; rearrange for x⁸, remember minus is plus in GF(2), and out drops a plain substitution rule:

That's the whole reduction rule. Next time a multiplication produces x⁸, just don't panic. Swap it for x⁴ + x³ + x + 1 and carry on. Test it on real bytes: 10000000 (that's x⁷) times 00000010 (that's x) gives x⁸, which the rule turns into x⁴ + x³ + x + 1 — the byte 00011011, better known as 0x1B, or 27.

7. Reed-Solomon, finally

Now stitch it together. To send n numbers and survive up to k erasures:

1. Pick a prime p larger than any value you're sending, and do all arithmetic mod p.
2. Treat your n numbers as the values of a polynomial: $f\left(0\right)$, $f\left(1\right)$, …, $f(n-1)$. Interpolation gives a unique polynomial f of degree n−1.
3. Evaluate f at k extra points and send all n+k values as labelled packets.

On the other end, the receiver might be missing up to k packets — but it still holds at least n. And n points uniquely determine a degree-(n−1) polynomial, so it re-interpolates the exact same f, then read $f\left(0\right)$ … $f(n-1)$ to recover the original numbers. Amazing, isn't it?

8. A worked example

Let's send the four numbers 2, 4, 3, 1 with k = 2, working mod 5. Lagrange interpolation using the points (0, 2), (1, 4), (2, 3), (3, 1) gives us a nice cubic equation:

Now just evaluate it at x = 0 … 5. The first four reproduce the data; the last two are our redundant packets:

Hence, we transmit the six numbers(packets) 2, 4, 3, 1, 0, 2. Say the channel eats packets m₁ and m₆:

Four points and a degree three polynomial — do interpolation mod 5 and you get the exact same f(x) = 2x³ + 2. Plug in x = 0 and x = 5 and the lost packets — 2 and 2 — walk right back in. The data was never really gone; it was implied by the curve the whole time.

9. The whole pipeline, playable

And here's the whole sequence all together, this time over GF(7). Just edit the message by clicking on cells (this bump their value), see how it gets interpolated into a polynomial and evaluated to packets, and click on packets to lose them. It will recover until you remove k packets; after that it won't be able to.

10. From erasures to errors

Erasures were the easy half: you know which packets vanished. The brutal case is a packet that arrives corrupted and pretends to be fine. Now you don't know the position or the value and a single liar, fed into interpolation, corrupt the whole polynomial.

The solution is the Berlekamp-Welch algorithm, and the cost is one additional factor of redundancy: to correct up to k errors we must send n + 2k points (versus n + k for erasures). The algorithm is brilliant. Consider an error-locator polynomial E(x) of degree k that is zero exactly at the corrupted positions. We don't know it yet, but we know it exists. Wherever a received value $r_i$ is honest, it agrees with the true polynomial, so:

And wherever a value is corrupted, $E\left(x_i\right)=0$ kills both sides anyway. So the equation $r_i·E\left(x_i\right)=Q\left(x_i\right)$ holds at every point, where $Q\left(x\right)=P\left(x\right)·E\left(x\right)$. We don't know P or E, but $r_i$ and $x_i$ are known — so each point gives one linear equation in the unknown coefficients of Q and E.

Count them: E has degree k, Q = P·E has degree n+k−1, and the n+2k received points give exactly enough equations to solve the system. Recover Q and E, then divide:

The corrupted packets simply never get a vote. E zeroes them out and the true message polynomial drops out of the division. See below table:

11. How hardware really does it: syndromes, Chien, Forney

Berlekamp-Welch has the nice narrative, yet it doesn't tend to be the decoder used in CD player, a QR scanner, or Voyager's transmitter. Those run an older, leaner pipeline that splits the job into three honest questions — is there a problem?, where?, and how bad is the problem? and answers each with a cheap, parallel pass instead of one big linear solve.

In practice, the RS code differs a little from what we have in our diagram for evaluation: The codewords consist of the multiples of a certain generator polynomial in such a way that each good codeword equals zero in 2k consecutive powers of a primitive element α — one field element whose powers α⁰, α¹, α², … run through every nonzero value before repeating. This is still the same code - just formulated the right way for the mechanics below.

The receiver holds $r\left(x\right)=c\left(x\right)+e\left(x\right)$, the true codeword plus whatever the channel added. Evaluate it at those special points and the codeword part vanishes ($c\left({\alpha }^j\right)=0$), leaving a pure readout of the error — the syndromes:

Each error throws in two unknowns: a locator $X_l={\alpha }^{i_l}$ saying which position it sits at, and a magnitude $Y_l$ saying how far it shoved the value. If all 2k syndromes come out zero, nothing's wrong — pass the word through untouched. Otherwise they're a compressed fingerprint of every error at once, before we know a single position.

Untangling that fingerprint is the job of the Berlekamp–Massey algorithm — a tidy O(k²) recurrence that turns the syndromes into the error-locator polynomial:

Λ is built so its roots land at exactly the inverse locations $1/X_l$. So we now hold a polynomial whose zeros pin the errors — we just have to find those zeros.

That's Chien search, and what a delightfully direct method it is: the field is finite, so try every element.^[𝓳]"just try everything" is brute force in a CS class. finite search space + parallel hardware = guilt-free 😎 Start x through α⁰, α¹, …, α^q−2, evaluate Λ at each, and every time it equals 0, there you've found a root — simply invert it and you have an error position. No fancy root-finding necessary, just brute-force testing; and with each being just multiply-adds, it's simple to test every possible location simultaneously in hardware.

This tell you where but not how big. Forney's algorithm pull out each magnitude directly. Construct the error-evaluator $\Omega \left(x\right)=S\left(x\right)\Lambda \left(x\right)\left(\mathrm{mod}{x}^{2k}\right)$ from the syndrome polynomial $S\left(x\right)$, take the formal derivative ${\Lambda }^{\prime }\left(x\right)$, and then:

The derivative is the trick: at a root of Λ it cancels every factor except the one belonging to this particular error term, resulting in one contribution only — so instead of solving a second linear system, just one calculation is needed per term. Simply subtract the reconstituted $Y_l$ from each and the received word reconstructs to the proper codeword. (Note: The coefficients in this case are those of the narrow-sense case; the particular coefficient will change depending on where you begin the roots.)

12. Where it actually lives

It’s not just a blackboard abstraction; chances are, it's at work somewhere near you right now:

• CDs & DVDs. The cross-interleaved Reed-Solomon distribution spreads your bits around the disc, and any scratches will destroy a collection of symbols instead of an entire packet - exactly what RS is immune to.
• QR codes. Four different error correction levels that can salvage nearly 30% of the code. That's why the company's logo in the center doesn't interfere.
• Deep space. Voyager 1 & 2 used Reed-Solomon to transmit images of outer planets through billions of kilometers of noise.^[𝓳]somewhere past the heliopause, a 1970s polynomial is still doing error correction. job security 😌
• Storage & networks. RAID-6, DSL, data-matrix barcodes — anywhere bits degrade or disappear, a RS code will fix them.

13. Closing thought

What i love about this one is its pedigree. Évariste Galois^[𝓳]died in a duel at 20, writing his theory the night before. the guy had a deadline and met it 😭 worked out finite fields as abstract structure, it answered a question that had nothing to do with telephones or spaceflight (neither of which had been invented). Cauchy reviewed the work and saw it mattered; almost no one saw how it would matter. A century later, this is what's storing your songs, your QR codes, and two probes heading out of our solar system now.

"Mathematicians often discover truths long before the world knows what to do with them. Progress begins when someone else recognizes those same truths hiding inside a different problem."

14. References

Where I picked this up — notes, papers, and the explainer that helps me to understand deep enough to reconstruct:

1. Rao, S., & others. CS70 Note 9 — Polynomials & Error-Correcting Codes. UC Berkeley.
2. Mathematics Stack Exchange. Is there some intuition for Lagrange interpolation formula? (Arturo Magidin’s answer on the on/off-switch construction).
3. Reed, I. S., & Solomon, G. (1960). Polynomial Codes Over Certain Finite Fields. Journal of the SIAM, 8(2), 300-304.
4. Wikipedia. Reed-Solomon error correction.
5. Wikipedia. Berlekamp-Welch algorithm.
6. Wikipedia. Chien search.
7. Wikipedia. Forney algorithm.
8. See also my earlier post, Hamming codes, for the finite-field / Galois background.

note: a few blocks here rephrased with an AI assistant for finding the right word, purely for clearer wording and explanation. the maths, the worked examples, and the playgrounds are the real deal.

Two version of C++ implementation of Reed-Solomon, converted to the browser environment for encoding your own byte sequences, damaging them, and observing how syndromes behave without having to use any compilers at all. basic is the evaluation encoder; prod is the systematic encoder with a generator polynomial and a syndrome check. open standalone ↗

#include <iostream>
          #include <vector>
          #include <cstdint>
          using namespace std;

          // GF(2^8): addition is XOR.
          uint8_t gf_add(uint8_t a, uint8_t b) { return a ^ b; }

          // multiply by x; reduce x^8 using x^8 = x^4 + x^3 + x + 1  (0x1B)
          uint8_t xtime(uint8_t value) {
              bool overflow = value & 0x80;
              value <<= 1;
              if (overflow) value ^= 0x1B;
              return value;
          }

          // Russian-peasant multiply: XOR shifted copies of x.
          uint8_t gf_mul(uint8_t x, uint8_t y) {
              uint8_t result = 0;
              while (y != 0) {
                  if (y & 1) result ^= x;   // include this shifted copy
                  x = xtime(x);             // x, 2x, 4x, ...
                  y >>= 1;
              }
              return result;
          }

          // evaluate P(x) at a point (ascending coefficients).
          uint8_t evaluate_polynomial(const vector<uint8_t>& poly, uint8_t x) {
              uint8_t result = 0, power = 1;
              for (uint8_t coefficient : poly) {
                  result = gf_add(result, gf_mul(coefficient, power));
                  power  = gf_mul(power, x);
              }
              return result;
          }

          // codeword[i] = P(i) for i = 0 .. message+parity-1 : every symbol is a point.
          vector<uint8_t> evaluation_encode(const vector<uint8_t>& message, int paritySymbols) {
              int total = message.size() + paritySymbols;
              vector<uint8_t> codeword;
              for (int i = 0; i < total; i++)
                  codeword.push_back(evaluate_polynomial(message, i));
              return codeword;
          }

          int main() {
              string text = "Hello";
              vector<uint8_t> message(text.begin(), text.end());
              auto codeword = evaluation_encode(message, 10);   // 5 + 10 = 15 points
          }
          

#include <iostream>
          #include <vector>
          #include <cstdint>
          #include <stdexcept>
          using namespace std;

          // ── GF(2^8) arithmetic ──
          uint8_t gf_add(uint8_t a, uint8_t b) { return a ^ b; }

          uint8_t xtime(uint8_t value) {
              bool overflow = value & 0x80;
              value <<= 1;
              if (overflow) value ^= 0x1B;          // x^8 = x^4 + x^3 + x + 1
              return value;
          }

          uint8_t gf_mul(uint8_t x, uint8_t y) {
              uint8_t result = 0;
              while (y != 0) {
                  if (y & 1) result ^= x;
                  x = xtime(x);
                  y >>= 1;
              }
              return result;
          }

          uint8_t gf_inverse(uint8_t value) {
              if (value == 0) throw runtime_error("zero has no multiplicative inverse");
              for (int candidate = 1; candidate < 256; candidate++)
                  if (gf_mul(value, candidate) == 1) return candidate;
              throw runtime_error("inverse not found");
          }

          uint8_t gf_div(uint8_t a, uint8_t b) { return gf_mul(a, gf_inverse(b)); }

          // ── polynomials, ascending order: {3,5,7} = 3 + 5x + 7x^2 ──
          vector<uint8_t> poly_mul(const vector<uint8_t>& a, const vector<uint8_t>& b) {
              vector<uint8_t> result(a.size() + b.size() - 1, 0);
              for (size_t i = 0; i < a.size(); i++)
                  for (size_t j = 0; j < b.size(); j++)            // x^i * x^j = x^(i+j)
                      result[i+j] = gf_add(result[i+j], gf_mul(a[i], b[j]));
              return result;
          }

          vector<uint8_t> poly_mod(vector<uint8_t> dividend, const vector<uint8_t>& divisor) {
              while (dividend.size() >= divisor.size()) {
                  uint8_t lead = dividend.back();
                  if (lead == 0) { dividend.pop_back(); continue; }
                  uint8_t factor = gf_div(lead, divisor.back());   // cancel the top term
                  for (size_t i = 0; i < divisor.size(); i++) {
                      int di = dividend.size() - divisor.size() + i;
                      dividend[di] = gf_add(dividend[di], gf_mul(divisor[i], factor));
                  }
                  while (!dividend.empty() && dividend.back() == 0) dividend.pop_back();
              }
              return dividend;                                      // remainder
          }

          // G(x) = (x + a)(x + a^2) ... (x + a^parity),  a = 2
          vector<uint8_t> generator_polynomial(int paritySymbols) {
              vector<uint8_t> generator = {1};
              uint8_t alpha = 2;
              for (int i = 0; i < paritySymbols; i++) {
                  generator = poly_mul(generator, {alpha, 1});
                  alpha = gf_mul(alpha, 2);
              }
              return generator;
          }

          vector<uint8_t> systematic_encode(const vector<uint8_t>& message, int paritySymbols) {
              auto generator = generator_polynomial(paritySymbols);

              vector<uint8_t> dividend(paritySymbols, 0);
              dividend.insert(dividend.end(), message.begin(), message.end());   // M(x) * x^parity
              auto parity = poly_mod(dividend, generator);
              parity.resize(paritySymbols, 0);

              // Parity goes in the LOW-order coefficients, message above it, so the
              // codeword stays divisible by G(x). (My first draft did message ‖ parity,
              // which left a clean word with NON-zero syndromes — the bug below.)
              vector<uint8_t> codeword = parity;
              codeword.insert(codeword.end(), message.begin(), message.end());
              return codeword;
          }

          uint8_t evaluate_polynomial(const vector<uint8_t>& poly, uint8_t x) {
              uint8_t result = 0, power = 1;
              for (uint8_t coefficient : poly) {
                  result = gf_add(result, gf_mul(coefficient, power));
                  power  = gf_mul(power, x);
              }
              return result;
          }

          // syndromes: evaluate the received word at a, a^2, ... a^parity (roots of G).
          vector<uint8_t> syndrome(const vector<uint8_t>& codeword, int paritySymbols) {
              vector<uint8_t> s;
              uint8_t alpha = 2;
              for (int i = 0; i < paritySymbols; i++) {
                  s.push_back(evaluate_polynomial(codeword, alpha));
                  alpha = gf_mul(alpha, 2);
              }
              return s;
          }

          // roots of the error-locator are the error positions.
          vector<int> chien_search(const vector<uint8_t>& locator) {
              vector<int> positions;
              uint8_t alpha = 1;
              for (int i = 0; i < 255; i++) {
                  if (evaluate_polynomial(locator, alpha) == 0) positions.push_back(i);
                  alpha = gf_mul(alpha, 2);
              }
              return positions;
          }

          // next stage — find Lambda(x), its roots, then the magnitudes (still stubbed).
          vector<uint8_t> berlekamp_massey(const vector<uint8_t>& syndromes);
          vector<uint8_t> forney(const vector<uint8_t>& s,
                                 const vector<uint8_t>& locator,
                                 const vector<int>& positions);

          vector<uint8_t> decode(vector<uint8_t> codeword, int paritySymbols) {
              auto s = syndrome(codeword, paritySymbols);

              bool hasError = false;
              for (uint8_t x : s) if (x != 0) { hasError = true; break; }
              if (!hasError) return codeword;            // syndromes all zero -> clean

              // auto locator    = berlekamp_massey(s);
              // auto positions  = chien_search(locator);
              // auto magnitudes = forney(s, locator, positions);
              // for (size_t i = 0; i < positions.size(); i++)
              //     codeword[positions[i]] ^= magnitudes[i];
              return codeword;                            // detection only, for now
          }

          int main() {
              string text = "Hello";
              vector<uint8_t> message(text.begin(), text.end());
              auto codeword = systematic_encode(message, 4);
              auto s = syndrome(codeword, 4);             // {0, 0, 0, 0} when clean
              auto out = decode(codeword, 4);
          }