What one needs to know in order to implement O(N*log(N)) Reed-Solomon error-correcting codes
Prerequsites:
- Complex arithmetic, in particular add/sub/mul, powers, roots and polar form
- Finite fields, with the same topics in mind
Topics:
- Discrete Fourier transform
- Fast Fourier transform and in particular Cooley–Tukey FFT algorithm as O(N*log(N)) algorithm implementing DFT
- Fast Number-Theoretic Transform as modified FFT employing the same add/sub/mul operations and unity roots, but in Galois Field
Once you grasped all these topics, you can grab some FFT implementation and convert it to NTT.
If source vector A = (A0, A1...A[n-1])
represents polynomial A(x) = A0 + A1*x + ... + A[n-1]*x**(n-1)
,
then NTT(A,a) computes polynomial values at the points 1, a, a**2 ... a**(n-1)
,
representing n
roots of power n
of 1 (if a
is a primitive root of power n
of 1).
It's equivalent to multiplying vector A
to Vandermonde Matrix defined by these points.
So, NTT converts polynomial coefficients to polynomial values at these specific points, and iNTT (inverse NTT) reverses this transformation.
iNTT is almost the same as NTT, i.e. iNTT(A,a) = NTT(A,1/a)/n
(with an element-wise division by n) reverses the effect of NTT(A,a)
.
This means that NTT(A,a)
provides the fast way to convert polynomial coefficients A
into values at a**i
points and vice versa.
The idea behind Reed-Solomon n+k
encoding (i.e. with n
source words and k
ecc words) is simple:
multiply source vector A
to some fixed matrix M[n*k]
and get resulting vector C
.
Decoding needs a mutiplication of combined vector (A,C) to inverse of combined matrix (I,M)
where I[n*n]
is identity matrix.
Of course, the decoding required only when we lost some values in A
, so real decoding combines alive values in A
and C
to vector of size n
(if more values are alive, they are just ignored), then deletes from (I,M)
combined matrix the k
lines corresponding to deleted values
and then mutiplies partial (A,C)
vector of size n
to inverse of partial (I,M)
matrix of size n*n
.
This operation computes the original A
vector. And of course, all computations are performed in some Galois Field to ensure
that all the intermediate results still fit to some 8-32 bit word.
Therefore, all that we need are guarantees that any n
lines of the (I,M)
combined matrix form a invertible matrix.
If that's true, then we have an MDS code - we can lose any k
words
and restore the original data from n
remaining ones.
There are so many ways to ensure that matrices are invertible:
- We can use non-systematic code, i.e. code where all
n+k
codes are computed rather than includen
original andk
computed codes. Non-systematic codes are ok, for example, for data transfer over noisy channels. So we just multiply sourceA
vector by Vandermonde(n+k)*n
matrix generated fromn+k
differenta[i]
numbers. It's guaranteed that anyn
differenta[i]
numbers form an invertible Vandermonde matrix, so we can restore from anyn
remaining words after a loss. - Plank proposed to start with Vandermonde
(n+k)*n
matrix and then apply the Gaussian elimination in order to convert it to some(I,M)
matrix. As far as we perform this operation only once per a lot of parity computations, we can ignore the time required by this operation. - PAR2 format employs
(I,V)
encoding matrix, i.e. it employs Vandermondek*n
matrix to computek
ecc words while employing the systematic code. Despite of special form ofa[i]
used in their Vandermonde matrix, the restoration matrix is sometimes non-invertible. But it seems to be a good compromise between the speed/complexity of computations and recovery strength. - Persicum proposed us to use the Cauchy matrix of the form
M[i,j] = 1/(n+i+j)
(and it probably was empoyed by RAR5). He said that it guarantees invertability of anyn*n
submatrix of(I,M)
, but i have no idea whether it's true.
The first and most obvious idea to perform RS encoding in O(N*log(N)) time is to compute non-systematic code:
extend n
source words with zeroes to n+k
, NTT them and send n+k
words produced to the channel. We are done!
This can be written as C = V(a)*ext(A)
where ext(A)
is a source vector A
extended with zeroes to n+k
elements,
C
is encoded vector to be sent to the channel, V(a)
is Vandermonde (n+k)*(n+k)
matix produced by powers of a
,
i.e. (1, a, a**2 ... a**(n+k-1))
and a
is a primitive root of 1 of power n+k
.
Converting this to systematic code means the equation (A,C) = V(a)*X
, so A
becomes a part of encoded word
and we need to compute only remaining C
part, and X
is some unknown vector we have to compute.