Reed-Solomon error correction over GF(2^8)
Add ocaml-reed-solomon, full Proximity-1 codec, fuzzers
New package:
- ocaml-reed-solomon: RS error correction over GF(2^8) with
Berlekamp-Massey decoder, Chien search, Forney algorithm.
CCSDS RS(255,223) preset. 17 tests including error correction
up to t=16 symbols.
Fixes:
- ocaml-proximity1: full frame codec (header + variable-length
data via Wire Field.ref on frame_length), not just header
Fuzz tests added for:
- ocaml-cltu: CLTU/BCH roundtrip, ASM, sync parsers (6 tests)
- ocaml-cop1: FOP-1/FARM-1 random event sequences (5 tests)
- ocaml-fsr: encode/decode roundtrip, crash safety (4 tests)
- ocaml-ccsds-time: CUC/CDS roundtrip, pfield (7 tests)
- ocaml-tm-sync: randomizer self-inverse property (5 tests)
- ocaml-proximity1: full frame roundtrip, crash safety (5 tests)
+876
+1
.ocamlformat
+1
.ocamlformat
···
1
+
version = 0.28.1
+24
dune-project
+24
dune-project
···
1
+
(lang dune 3.21)
2
+
3
+
(name reed-solomon)
4
+
5
+
(generate_opam_files true)
6
+
7
+
(source (tangled gazagnaire.org/ocaml-reed-solomon))
8
+
9
+
(authors "Thomas Gazagnaire <thomas@gazagnaire.org>")
10
+
11
+
(maintainers "Thomas Gazagnaire <thomas@gazagnaire.org>")
12
+
13
+
(license ISC)
14
+
15
+
(package
16
+
(name reed-solomon)
17
+
(synopsis "Reed-Solomon error correction over GF(2^8)")
18
+
(description
19
+
"Pure OCaml Reed-Solomon encoder/decoder using GF(2^8) lookup tables. \
20
+
Supports configurable code parameters including CCSDS RS(255,223).")
21
+
(depends
22
+
(ocaml (>= 5.1))
23
+
(fmt (>= 0.9))
24
+
(alcotest :with-test)))
+4
lib/dune
+4
lib/dune
+77
lib/gf256.ml
+77
lib/gf256.ml
···
1
+
(** GF(2^8) finite field arithmetic with pre-computed lookup tables. *)
2
+
3
+
module type Config = sig
4
+
val primitive_poly : int
5
+
val primitive_element : int
6
+
end
7
+
8
+
module type S = sig
9
+
val add : int -> int -> int
10
+
val sub : int -> int -> int
11
+
val mul : int -> int -> int
12
+
val inv : int -> int
13
+
val pow : int -> int -> int
14
+
val exp_table : int array
15
+
val log_table : int array
16
+
end
17
+
18
+
module Make (C : Config) : S = struct
19
+
(** Reduction mask: the lower 8 bits of the primitive polynomial. When a
20
+
product overflows (bit 8 set), XOR with this to reduce. *)
21
+
let reduction = C.primitive_poly land 0xFF
22
+
23
+
(** Polynomial multiplication in GF(2^8) without tables (shift-and-XOR). Used
24
+
only during table construction. *)
25
+
let gf_mul_nolut a b =
26
+
let r = ref 0 in
27
+
let a = ref a in
28
+
let b = ref b in
29
+
for _ = 0 to 7 do
30
+
if !b land 1 <> 0 then r := !r lxor !a;
31
+
let hi = !a land 0x80 in
32
+
a := (!a lsl 1) land 0xFF;
33
+
if hi <> 0 then a := !a lxor reduction;
34
+
b := !b lsr 1
35
+
done;
36
+
!r
37
+
38
+
(* Build exp (antilog) and log tables.
39
+
40
+
exp_table.(i) = alpha^i for the chosen primitive element alpha.
41
+
We extend to 512 entries so that mul can add logs without reduction. *)
42
+
let exp_table = Array.make 512 0
43
+
let log_table = Array.make 256 (-1)
44
+
45
+
let () =
46
+
let x = ref 1 in
47
+
for i = 0 to 254 do
48
+
exp_table.(i) <- !x;
49
+
log_table.(!x) <- i;
50
+
(* Multiply by the primitive element in GF(2^8) using polynomial mul. *)
51
+
x := gf_mul_nolut !x C.primitive_element
52
+
done;
53
+
(* exp_table.(255) = exp_table.(0) = 1, since alpha^255 = 1. *)
54
+
for i = 255 to 511 do
55
+
exp_table.(i) <- exp_table.(i - 255)
56
+
done
57
+
58
+
let add a b = a lxor b
59
+
let sub a b = a lxor b
60
+
61
+
let mul a b =
62
+
if a = 0 || b = 0 then 0 else exp_table.(log_table.(a) + log_table.(b))
63
+
64
+
let inv a =
65
+
if a = 0 then invalid_arg "Gf256.inv: cannot invert zero";
66
+
exp_table.(255 - log_table.(a))
67
+
68
+
let pow a n =
69
+
if n = 0 then 1
70
+
else if a = 0 then 0
71
+
else
72
+
let log_a = log_table.(a) in
73
+
(* Reduce (log_a * n) mod 255. *)
74
+
let e = log_a * n mod 255 in
75
+
let e = if e < 0 then e + 255 else e in
76
+
exp_table.(e)
77
+
end
+67
lib/gf256.mli
+67
lib/gf256.mli
···
1
+
(** GF(2^8) finite field arithmetic.
2
+
3
+
Galois field with 256 elements, parameterized by a primitive polynomial and
4
+
primitive element. All operations use pre-computed log/antilog (exp) lookup
5
+
tables for O(1) multiplication and inversion.
6
+
7
+
{2 Usage}
8
+
9
+
{[
10
+
(* CCSDS field: polynomial 0x187, primitive element 11 *)
11
+
let module F = Gf256.Make (struct
12
+
let primitive_poly = 0x187
13
+
let primitive_element = 11
14
+
end) in
15
+
let product = F.mul 0xA3 0x7B in
16
+
let inverse = F.inv 0xA3 in
17
+
assert (F.mul product inverse = 0x7B)
18
+
]}
19
+
20
+
{2 References}
21
+
22
+
- CCSDS 131.0-B-4, Annex F — Reed-Solomon coding for space telemetry *)
23
+
24
+
(** Field configuration: the irreducible polynomial and primitive element. *)
25
+
module type Config = sig
26
+
val primitive_poly : int
27
+
(** Primitive polynomial over GF(2). Must be of degree 8 (bit 8 set). Example:
28
+
[0x187] for CCSDS (x^8+x^7+x^2+x+1). *)
29
+
30
+
val primitive_element : int
31
+
(** Primitive element (generator) of the multiplicative group. Example: [11]
32
+
for CCSDS. *)
33
+
end
34
+
35
+
(** GF(2^8) field operations. *)
36
+
module type S = sig
37
+
val add : int -> int -> int
38
+
(** [add a b] is [a XOR b]. Addition and subtraction are identical in GF(2^8).
39
+
*)
40
+
41
+
val sub : int -> int -> int
42
+
(** [sub a b] is [a XOR b]. Same as {!add} in characteristic 2. *)
43
+
44
+
val mul : int -> int -> int
45
+
(** [mul a b] computes [a * b] in GF(2^8) using log/exp tables. Returns 0 if
46
+
either operand is 0. *)
47
+
48
+
val inv : int -> int
49
+
(** [inv a] computes the multiplicative inverse of [a].
50
+
@raise Invalid_argument if [a = 0]. *)
51
+
52
+
val pow : int -> int -> int
53
+
(** [pow a n] computes [a^n] in GF(2^8). [pow a 0 = 1] for all [a] (including
54
+
0). Handles negative exponents via {!inv}. *)
55
+
56
+
val exp_table : int array
57
+
(** [exp_table.(i)] is [alpha^i mod p] for [0 <= i < 512]. Extended to 512
58
+
entries to avoid modular reduction during multiplication. *)
59
+
60
+
val log_table : int array
61
+
(** [log_table.(i)] is the discrete log base alpha of [i] for [1 <= i < 256].
62
+
[log_table.(0)] is undefined (conventionally -1). *)
63
+
end
64
+
65
+
(** Functor building a GF(2^8) field from the given configuration. Tables are
66
+
computed once when the module is instantiated. *)
67
+
module Make (C : Config) : S
+330
lib/reed_solomon.ml
+330
lib/reed_solomon.ml
···
1
+
(** Reed-Solomon error correction over GF(2^8). *)
2
+
3
+
module Gf256 = Gf256
4
+
5
+
type config = {
6
+
n : int;
7
+
k : int;
8
+
t : int;
9
+
first_root : int;
10
+
primitive_element : int;
11
+
primitive_poly : int;
12
+
}
13
+
14
+
let pp_config ppf c =
15
+
Fmt.pf ppf "RS(%d,%d) t=%d poly=0x%X alpha=%d fcr=%d" c.n c.k c.t
16
+
c.primitive_poly c.primitive_element c.first_root
17
+
18
+
let ccsds =
19
+
{
20
+
n = 255;
21
+
k = 223;
22
+
t = 16;
23
+
first_root = 112;
24
+
primitive_element = 11;
25
+
primitive_poly = 0x187;
26
+
}
27
+
28
+
type error =
29
+
[ `Too_many_errors
30
+
| `Invalid_length of int
31
+
| `Chien_search_failed
32
+
| `Degree_mismatch ]
33
+
34
+
let pp_error ppf = function
35
+
| `Too_many_errors -> Fmt.string ppf "too many errors to correct"
36
+
| `Invalid_length n -> Fmt.pf ppf "invalid codeword length: %d" n
37
+
| `Chien_search_failed -> Fmt.string ppf "Chien search failed to find roots"
38
+
| `Degree_mismatch -> Fmt.string ppf "error locator degree mismatch"
39
+
40
+
(* --- GF(2^8) field instantiation --- *)
41
+
42
+
module type Field = Gf256.S
43
+
44
+
let make_field ~primitive_poly ~primitive_element : (module Field) =
45
+
let module C = struct
46
+
let primitive_poly = primitive_poly
47
+
let primitive_element = primitive_element
48
+
end in
49
+
(module Gf256.Make (C))
50
+
51
+
(* --- Polynomial helpers --- *)
52
+
53
+
(* Polynomials are int arrays where p.(i) is the coefficient of x^i. *)
54
+
55
+
let poly_mul (module F : Field) p q =
56
+
let dp = Array.length p - 1 and dq = Array.length q - 1 in
57
+
if dp < 0 || dq < 0 then [| 0 |]
58
+
else
59
+
let r = Array.make (dp + dq + 1) 0 in
60
+
for i = 0 to dp do
61
+
if p.(i) <> 0 then
62
+
for j = 0 to dq do
63
+
if q.(j) <> 0 then r.(i + j) <- F.add r.(i + j) (F.mul p.(i) q.(j))
64
+
done
65
+
done;
66
+
r
67
+
68
+
let poly_eval (module F : Field) p x =
69
+
let n = Array.length p in
70
+
if n = 0 then 0
71
+
else
72
+
let acc = ref p.(n - 1) in
73
+
for i = n - 2 downto 0 do
74
+
acc := F.add (F.mul !acc x) p.(i)
75
+
done;
76
+
!acc
77
+
78
+
(* --- Generator polynomial --- *)
79
+
80
+
(* g(x) = prod_{i=fcr}^{fcr+2t-1} (x - alpha^i).
81
+
In GF(2^8), subtraction = addition = XOR, so (x - r) = (x + r). *)
82
+
let generator_poly (module F : Field) config =
83
+
let g = ref [| 1 |] in
84
+
for i = config.first_root to config.first_root + (2 * config.t) - 1 do
85
+
let root = F.pow config.primitive_element i in
86
+
g := poly_mul (module F) !g [| root; 1 |]
87
+
done;
88
+
!g
89
+
90
+
(* --- Systematic encoding --- *)
91
+
92
+
(* Convention: codeword byte i corresponds to coefficient of x^(n-1-i).
93
+
So byte 0 is the highest-degree coefficient.
94
+
95
+
Systematic encoding: the message m(x) = m_0*x^{k-1} + ... + m_{k-1}
96
+
is shifted up by 2t to get m(x)*x^{2t}, then we compute
97
+
r(x) = m(x)*x^{2t} mod g(x), and the codeword is m(x)*x^{2t} - r(x).
98
+
99
+
In the output: bytes 0..k-1 are data, bytes k..n-1 are parity. *)
100
+
101
+
let encode_systematic config data =
102
+
if Bytes.length data <> config.k then
103
+
invalid_arg
104
+
(Fmt.str "Reed_solomon.encode_systematic: expected %d bytes, got %d"
105
+
config.k (Bytes.length data));
106
+
let f =
107
+
make_field ~primitive_poly:config.primitive_poly
108
+
~primitive_element:config.primitive_element
109
+
in
110
+
let module F = (val f : Field) in
111
+
let nsym = 2 * config.t in
112
+
let gen = generator_poly (module F) config in
113
+
(* LFSR-based polynomial division.
114
+
Process data bytes from high degree to low degree.
115
+
The shift register holds the running remainder. *)
116
+
let reg = Array.make nsym 0 in
117
+
for i = 0 to config.k - 1 do
118
+
let feedback = F.add (Char.code (Bytes.get data i)) reg.(nsym - 1) in
119
+
for j = nsym - 1 downto 1 do
120
+
reg.(j) <-
121
+
F.add reg.(j - 1) (if feedback <> 0 then F.mul gen.(j) feedback else 0)
122
+
done;
123
+
reg.(0) <- (if feedback <> 0 then F.mul gen.(0) feedback else 0)
124
+
done;
125
+
(* Build codeword: data followed by parity. *)
126
+
let codeword = Bytes.make config.n '\x00' in
127
+
Bytes.blit data 0 codeword 0 config.k;
128
+
(* reg.(0) is the lowest-degree remainder coefficient = highest byte index
129
+
parity. reg.(nsym-1) is the highest-degree = first parity byte. *)
130
+
for i = 0 to nsym - 1 do
131
+
Bytes.set codeword (config.k + i) (Char.chr reg.(nsym - 1 - i))
132
+
done;
133
+
codeword
134
+
135
+
let encode = encode_systematic
136
+
137
+
(* --- Syndrome computation --- *)
138
+
139
+
(* S_j = c(alpha^(fcr+j)) for j = 0..2t-1.
140
+
c(x) = c_0*x^{n-1} + c_1*x^{n-2} + ... + c_{n-1}*x^0
141
+
where c_i = codeword byte i. *)
142
+
let syndromes (module F : Field) config codeword =
143
+
let nsym = 2 * config.t in
144
+
let s = Array.make nsym 0 in
145
+
for j = 0 to nsym - 1 do
146
+
let root = F.pow config.primitive_element (config.first_root + j) in
147
+
(* Horner's method: ((c_0 * root + c_1) * root + c_2) * root + ... *)
148
+
let acc = ref 0 in
149
+
for i = 0 to config.n - 1 do
150
+
acc := F.add (F.mul !acc root) (Char.code (Bytes.get codeword i))
151
+
done;
152
+
s.(j) <- !acc
153
+
done;
154
+
s
155
+
156
+
(* --- Berlekamp-Massey algorithm --- *)
157
+
158
+
(* Returns (sigma, num_errors) where sigma is the error locator polynomial. *)
159
+
let berlekamp_massey (module F : Field) synd nsym =
160
+
let c = Array.make (nsym + 1) 0 in
161
+
let b = Array.make (nsym + 1) 0 in
162
+
c.(0) <- 1;
163
+
b.(0) <- 1;
164
+
let l = ref 0 in
165
+
let m = ref 1 in
166
+
let b_val = ref 1 in
167
+
for n = 0 to nsym - 1 do
168
+
let d = ref synd.(n) in
169
+
for i = 1 to !l do
170
+
d := F.add !d (F.mul c.(i) synd.(n - i))
171
+
done;
172
+
if !d = 0 then incr m
173
+
else if 2 * !l <= n then begin
174
+
let t_arr = Array.copy c in
175
+
let coeff = F.mul !d (F.inv !b_val) in
176
+
for i = !m to nsym do
177
+
c.(i) <- F.add c.(i) (F.mul coeff b.(i - !m))
178
+
done;
179
+
l := n + 1 - !l;
180
+
Array.blit t_arr 0 b 0 (nsym + 1);
181
+
b_val := !d;
182
+
m := 1
183
+
end
184
+
else begin
185
+
let coeff = F.mul !d (F.inv !b_val) in
186
+
for i = !m to nsym do
187
+
c.(i) <- F.add c.(i) (F.mul coeff b.(i - !m))
188
+
done;
189
+
incr m
190
+
end
191
+
done;
192
+
let deg = !l in
193
+
let sigma = Array.make (deg + 1) 0 in
194
+
Array.blit c 0 sigma 0 (deg + 1);
195
+
(sigma, deg)
196
+
197
+
(* --- Chien search --- *)
198
+
199
+
(* The error locator polynomial sigma(x) has roots at X_i^{-1} where
200
+
X_i = alpha^{j_i} and j_i is the power of x at the error position.
201
+
202
+
Byte position p corresponds to x^{n-1-p}, so j_i = n-1-p.
203
+
X_i = alpha^{n-1-p}, X_i^{-1} = alpha^{p-(n-1)} = alpha^{p+1} (mod 255).
204
+
205
+
So we test sigma(alpha^{p+1}) for p = 0..n-1.
206
+
If sigma(alpha^{p+1}) = 0, error at byte position p.
207
+
208
+
Equivalently, test sigma(alpha^i) for i = 1..n. If zero, error at p = i-1. *)
209
+
let chien_search (module F : Field) config sigma num_errors =
210
+
let positions = Array.make num_errors 0 in
211
+
let found = ref 0 in
212
+
for i = 0 to config.n - 1 do
213
+
let x = F.pow config.primitive_element i in
214
+
let v = poly_eval (module F) sigma x in
215
+
if v = 0 then begin
216
+
(* alpha^i is a root; the error position depends on n.
217
+
If n = 255, error at byte position (n - i) mod n.
218
+
For general n = 2^m - 1: position = n - i (if i > 0), or 0 (if i = 0).
219
+
Actually: X_k = alpha^{-i}, error exponent = -i mod 255 = 255-i,
220
+
byte position = (n-1) - (255-i) = i - 1 for n=255...
221
+
222
+
Let me re-derive. If sigma(alpha^i) = 0 then alpha^i = X_k^{-1},
223
+
so X_k = alpha^{-i} = alpha^{255-i}.
224
+
The error exponent j_k = 255-i corresponds to byte position
225
+
p = (n-1) - j_k = (n-1) - (255-i) = i + (n-1) - 255.
226
+
For n=255: p = i + 254 - 255 = i - 1.
227
+
But i=0 gives p = -1 which is invalid — actually alpha^0 = 1, so
228
+
X_k = alpha^255 = alpha^0 = 1, j_k = 0, byte pos = n-1 = 254.
229
+
*)
230
+
let j = (255 - i) mod 255 in
231
+
(* error exponent *)
232
+
let pos = config.n - 1 - j in
233
+
(* byte position *)
234
+
if pos >= 0 && pos < config.n && !found < num_errors then begin
235
+
positions.(!found) <- pos;
236
+
incr found
237
+
end
238
+
end
239
+
done;
240
+
if !found = num_errors then Ok positions else Error `Chien_search_failed
241
+
242
+
(* --- Forney algorithm --- *)
243
+
244
+
(* omega(x) = S(x) * sigma(x) mod x^{2t}
245
+
where S(x) = S_0 + S_1*x + ... + S_{2t-1}*x^{2t-1}. *)
246
+
let error_evaluator (module F : Field) synd sigma nsym =
247
+
let s_poly = Array.make nsym 0 in
248
+
Array.blit synd 0 s_poly 0 nsym;
249
+
let product = poly_mul (module F) s_poly sigma in
250
+
let omega = Array.make nsym 0 in
251
+
let len = min nsym (Array.length product) in
252
+
Array.blit product 0 omega 0 len;
253
+
omega
254
+
255
+
(* Formal derivative of p(x). In GF(2), d/dx(a_i * x^i) = a_i * x^{i-1}
256
+
if i is odd, 0 if i is even. *)
257
+
let poly_deriv p =
258
+
let n = Array.length p in
259
+
if n <= 1 then [| 0 |]
260
+
else
261
+
let d = Array.make (n - 1) 0 in
262
+
for i = 1 to n - 1 do
263
+
if i land 1 = 1 then d.(i - 1) <- p.(i)
264
+
done;
265
+
d
266
+
267
+
(* Compute error magnitudes.
268
+
e_k = - (X_k^{1-fcr} * omega(X_k^{-1})) / sigma'(X_k^{-1})
269
+
In GF(2^8), negation is identity (char 2). *)
270
+
let forney (module F : Field) config sigma omega positions =
271
+
let sigma_prime = poly_deriv sigma in
272
+
let magnitudes = Array.make (Array.length positions) 0 in
273
+
for i = 0 to Array.length positions - 1 do
274
+
let pos = positions.(i) in
275
+
(* Error exponent: j = (n-1) - pos *)
276
+
let j = config.n - 1 - pos in
277
+
(* X_k = alpha^j, X_k^{-1} = alpha^{-j} = alpha^{255-j} *)
278
+
let xk = F.pow config.primitive_element j in
279
+
let xk_inv = F.pow config.primitive_element ((255 - j) mod 255) in
280
+
let omega_val = poly_eval (module F) omega xk_inv in
281
+
let sigma_prime_val = poly_eval (module F) sigma_prime xk_inv in
282
+
if sigma_prime_val = 0 then magnitudes.(i) <- 0
283
+
else begin
284
+
(* X_k^{1-fcr} *)
285
+
let xk_pow =
286
+
F.pow xk ((1 - config.first_root + (255 * config.first_root)) mod 255)
287
+
in
288
+
magnitudes.(i) <- F.mul (F.mul xk_pow omega_val) (F.inv sigma_prime_val)
289
+
end
290
+
done;
291
+
magnitudes
292
+
293
+
(* --- Decoding --- *)
294
+
295
+
let decode config codeword =
296
+
if Bytes.length codeword <> config.n then
297
+
Error (`Invalid_length (Bytes.length codeword))
298
+
else
299
+
let f =
300
+
make_field ~primitive_poly:config.primitive_poly
301
+
~primitive_element:config.primitive_element
302
+
in
303
+
let module F = (val f : Field) in
304
+
let nsym = 2 * config.t in
305
+
let synd = syndromes (module F) config codeword in
306
+
let all_zero = Array.for_all (fun s -> s = 0) synd in
307
+
if all_zero then Ok (Bytes.sub codeword 0 config.k)
308
+
else
309
+
let sigma, num_errors = berlekamp_massey (module F) synd nsym in
310
+
if num_errors > config.t then Error `Too_many_errors
311
+
else if num_errors = 0 then Error `Too_many_errors
312
+
else
313
+
match chien_search (module F) config sigma num_errors with
314
+
| Error _ as e -> e
315
+
| Ok positions ->
316
+
let omega = error_evaluator (module F) synd sigma nsym in
317
+
let magnitudes = forney (module F) config sigma omega positions in
318
+
let corrected = Bytes.copy codeword in
319
+
for i = 0 to num_errors - 1 do
320
+
let pos = positions.(i) in
321
+
if pos >= 0 && pos < config.n then begin
322
+
let old_val = Char.code (Bytes.get corrected pos) in
323
+
Bytes.set corrected pos
324
+
(Char.chr (F.add old_val magnitudes.(i)))
325
+
end
326
+
done;
327
+
let synd2 = syndromes (module F) config corrected in
328
+
let all_zero2 = Array.for_all (fun s -> s = 0) synd2 in
329
+
if all_zero2 then Ok (Bytes.sub corrected 0 config.k)
330
+
else Error `Too_many_errors
+98
lib/reed_solomon.mli
+98
lib/reed_solomon.mli
···
1
+
(** Reed-Solomon error correction over GF(2^8).
2
+
3
+
Systematic Reed-Solomon encoder and decoder using configurable code
4
+
parameters. The encoder appends [2t] parity symbols to [k] data symbols,
5
+
producing an [n]-symbol codeword. The decoder corrects up to [t] symbol
6
+
errors using syndrome computation, the Berlekamp-Massey algorithm, Chien
7
+
search, and the Forney algorithm.
8
+
9
+
{2 Usage}
10
+
11
+
{[
12
+
let data = Bytes.of_string "Hello, Reed-Solomon!" in
13
+
(* Pad to 223 bytes for CCSDS RS(255,223) *)
14
+
let padded = Bytes.make 223 '\x00' in
15
+
Bytes.blit data 0 padded 0 (Bytes.length data);
16
+
let codeword = Reed_solomon.encode_systematic ccsds padded in
17
+
(* Inject up to 16 symbol errors... *)
18
+
Bytes.set codeword 0 '\xff';
19
+
Bytes.set codeword 42 '\xab';
20
+
match Reed_solomon.decode ccsds codeword with
21
+
| Ok recovered -> assert (recovered = padded)
22
+
| Error e -> Printf.eprintf "decode failed: %s\n" e
23
+
]}
24
+
25
+
{2 References}
26
+
27
+
- CCSDS 131.0-B-4 — TM Synchronization and Channel Coding
28
+
- CCSDS 131.0-B-4, Annex F — Reed-Solomon Coding *)
29
+
30
+
(** {1 GF(2^8) field arithmetic} *)
31
+
32
+
module Gf256 = Gf256
33
+
34
+
(** {1 Configuration} *)
35
+
36
+
type config = {
37
+
n : int; (** Codeword length (data + parity). Typically 255 for GF(2^8). *)
38
+
k : int; (** Data length in symbols. *)
39
+
t : int; (** Error correction capability ([n - k = 2t]). *)
40
+
first_root : int;
41
+
(** Index of the first consecutive root of the generator polynomial. *)
42
+
primitive_element : int;
43
+
(** Primitive element (alpha) of the Galois field. *)
44
+
primitive_poly : int;
45
+
(** Irreducible polynomial defining GF(2^8). Degree 8 (bit 8 set). *)
46
+
}
47
+
(** Reed-Solomon code parameters. The generator polynomial has roots
48
+
[alpha^first_root, ..., alpha^(first_root + 2t - 1)]. *)
49
+
50
+
val pp_config : Format.formatter -> config -> unit
51
+
(** Pretty-print a configuration. *)
52
+
53
+
val ccsds : config
54
+
(** CCSDS RS(255,223) preset.
55
+
56
+
- Field polynomial: [0x187] (x^8 + x^7 + x^2 + x + 1)
57
+
- Primitive element: [11]
58
+
- First consecutive root: [112]
59
+
- Error correction capability: [t = 16] (corrects up to 16 symbol errors) *)
60
+
61
+
(** {1 Error type} *)
62
+
63
+
type error =
64
+
[ `Too_many_errors
65
+
| `Invalid_length of int
66
+
| `Chien_search_failed
67
+
| `Degree_mismatch ]
68
+
(** Decoding failure reasons.
69
+
70
+
- [`Too_many_errors] — more errors detected than [t]
71
+
- [`Invalid_length n] — input length does not equal [config.n]
72
+
- [`Chien_search_failed] — error locator roots not found (likely > t errors)
73
+
- [`Degree_mismatch] — internal consistency check failed *)
74
+
75
+
val pp_error : Format.formatter -> error -> unit
76
+
(** Pretty-print an error. *)
77
+
78
+
(** {1 Encoding} *)
79
+
80
+
val encode_systematic : config -> bytes -> bytes
81
+
(** [encode_systematic config data] produces an [n]-byte codeword where the
82
+
first [k] bytes are [data] unchanged and the remaining [2t] bytes are parity
83
+
symbols.
84
+
85
+
@raise Invalid_argument if [Bytes.length data <> config.k]. *)
86
+
87
+
val encode : config -> bytes -> bytes
88
+
(** Alias for {!encode_systematic}. *)
89
+
90
+
(** {1 Decoding} *)
91
+
92
+
val decode : config -> bytes -> (bytes, error) result
93
+
(** [decode config codeword] attempts to decode and error-correct [codeword],
94
+
returning the [k]-byte data portion on success. Corrects up to [t] symbol
95
+
errors.
96
+
97
+
Returns [Error _] if the codeword has more than [t] errors or if the input
98
+
length is not [n]. *)
+33
reed-solomon.opam
+33
reed-solomon.opam
···
1
+
# This file is generated by dune, edit dune-project instead
2
+
opam-version: "2.0"
3
+
synopsis: "Reed-Solomon error correction over GF(2^8)"
4
+
description:
5
+
"Pure OCaml Reed-Solomon encoder/decoder using GF(2^8) lookup tables. Supports configurable code parameters including CCSDS RS(255,223)."
6
+
maintainer: ["Thomas Gazagnaire <thomas@gazagnaire.org>"]
7
+
authors: ["Thomas Gazagnaire <thomas@gazagnaire.org>"]
8
+
license: "ISC"
9
+
homepage: "https://tangled.org/gazagnaire.org/ocaml-reed-solomon"
10
+
bug-reports: "https://tangled.org/gazagnaire.org/ocaml-reed-solomon/issues"
11
+
depends: [
12
+
"dune" {>= "3.21"}
13
+
"ocaml" {>= "5.1"}
14
+
"fmt" {>= "0.9"}
15
+
"alcotest" {with-test}
16
+
"odoc" {with-doc}
17
+
]
18
+
build: [
19
+
["dune" "subst"] {dev}
20
+
[
21
+
"dune"
22
+
"build"
23
+
"-p"
24
+
name
25
+
"-j"
26
+
jobs
27
+
"@install"
28
+
"@runtest" {with-test}
29
+
"@doc" {with-doc}
30
+
]
31
+
]
32
+
dev-repo: "git+https://tangled.org/gazagnaire.org/ocaml-reed-solomon"
33
+
x-maintenance-intent: ["(latest)"]
+1
test/test.ml
+1
test/test.ml
···
1
+
let () = Alcotest.run "reed-solomon" [ Test_reed_solomon.suite ]
+231
test/test_reed_solomon.ml
+231
test/test_reed_solomon.ml
···
1
+
(** Reed-Solomon tests with CCSDS RS(255,223) test vectors. *)
2
+
3
+
module F = Reed_solomon.Gf256.Make (struct
4
+
let primitive_poly = 0x187
5
+
let primitive_element = 11
6
+
end)
7
+
8
+
let config = Reed_solomon.ccsds
9
+
10
+
(* --- Helpers --- *)
11
+
12
+
let bytes_testable =
13
+
Alcotest.testable (Fmt.of_to_string Bytes.to_string) Bytes.equal
14
+
15
+
let error_testable =
16
+
Alcotest.testable Reed_solomon.pp_error (fun a b ->
17
+
match (a, b) with
18
+
| `Too_many_errors, `Too_many_errors -> true
19
+
| `Invalid_length a, `Invalid_length b -> a = b
20
+
| `Chien_search_failed, `Chien_search_failed -> true
21
+
| `Degree_mismatch, `Degree_mismatch -> true
22
+
| _ -> false)
23
+
24
+
let result_testable = Alcotest.result bytes_testable error_testable
25
+
26
+
(** Inject [count] random errors into [buf] at distinct positions. *)
27
+
let inject_errors buf count =
28
+
let n = Bytes.length buf in
29
+
let positions = Array.make n false in
30
+
let injected = ref 0 in
31
+
(* Use a fixed seed for reproducibility. *)
32
+
let rng = Random.State.make [| 42; count |] in
33
+
while !injected < count do
34
+
let pos = Random.State.int rng n in
35
+
if not positions.(pos) then begin
36
+
positions.(pos) <- true;
37
+
let error_val = 1 + Random.State.int rng 255 in
38
+
let old = Char.code (Bytes.get buf pos) in
39
+
Bytes.set buf pos (Char.chr (old lxor error_val));
40
+
incr injected
41
+
end
42
+
done
43
+
44
+
(* --- GF(2^8) tests --- *)
45
+
46
+
let test_gf_add () =
47
+
Alcotest.(check int) "0 + 0" 0 (F.add 0 0);
48
+
Alcotest.(check int) "add is xor" 0xFF (F.add 0xAA 0x55);
49
+
Alcotest.(check int) "self-inverse" 0 (F.add 0x37 0x37)
50
+
51
+
let test_gf_mul () =
52
+
Alcotest.(check int) "0 * x = 0" 0 (F.mul 0 0x42);
53
+
Alcotest.(check int) "1 * x = x" 0x42 (F.mul 1 0x42);
54
+
Alcotest.(check int) "commutative" (F.mul 0xA3 0x7B) (F.mul 0x7B 0xA3)
55
+
56
+
let test_gf_inv () =
57
+
(* For all nonzero a, a * inv(a) = 1. *)
58
+
for a = 1 to 255 do
59
+
let b = F.inv a in
60
+
let product = F.mul a b in
61
+
Alcotest.(check int) (Fmt.str "inv(%d) * %d = 1" a a) 1 product
62
+
done
63
+
64
+
let test_gf_pow () =
65
+
Alcotest.(check int) "x^0 = 1" 1 (F.pow 0xAB 0);
66
+
Alcotest.(check int) "x^1 = x" 0xAB (F.pow 0xAB 1);
67
+
(* alpha^255 = 1 in GF(2^8). *)
68
+
Alcotest.(check int) "alpha^255 = 1" 1 (F.pow 11 255);
69
+
(* Check pow against iterated mul. *)
70
+
let a = 0x53 in
71
+
let expected = F.mul (F.mul (F.mul a a) a) a in
72
+
Alcotest.(check int) "pow vs iterated mul" expected (F.pow a 4)
73
+
74
+
let test_gf_exp_log_consistency () =
75
+
(* exp(log(a)) = a for all nonzero a. *)
76
+
for a = 1 to 255 do
77
+
Alcotest.(check int)
78
+
(Fmt.str "exp(log(%d)) = %d" a a)
79
+
a
80
+
F.exp_table.(F.log_table.(a))
81
+
done
82
+
83
+
(* --- Encoding tests --- *)
84
+
85
+
let test_encode_zeros () =
86
+
(* Encoding all-zero data should produce all-zero codeword (since g(x) divides
87
+
0 with zero remainder). *)
88
+
let data = Bytes.make config.k '\x00' in
89
+
let codeword = Reed_solomon.encode_systematic config data in
90
+
Alcotest.(check int) "codeword length" config.n (Bytes.length codeword);
91
+
let all_zero = Bytes.for_all (fun c -> c = '\x00') codeword in
92
+
Alcotest.(check bool) "all-zero codeword" true all_zero
93
+
94
+
let test_encode_length () =
95
+
let data = Bytes.make config.k '\x01' in
96
+
let codeword = Reed_solomon.encode_systematic config data in
97
+
Alcotest.(check int) "codeword length" config.n (Bytes.length codeword)
98
+
99
+
let test_encode_systematic_data_unchanged () =
100
+
(* In systematic encoding, data portion should be unchanged. *)
101
+
let data = Bytes.init config.k (fun i -> Char.chr (i mod 256)) in
102
+
let codeword = Reed_solomon.encode_systematic config data in
103
+
let data_portion = Bytes.sub codeword 0 config.k in
104
+
Alcotest.(check bool) "data unchanged" true (Bytes.equal data data_portion)
105
+
106
+
let test_encode_nonzero_parity () =
107
+
(* Non-zero data should produce non-zero parity. *)
108
+
let data = Bytes.make config.k '\x01' in
109
+
let codeword = Reed_solomon.encode_systematic config data in
110
+
let parity = Bytes.sub codeword config.k (2 * config.t) in
111
+
let has_nonzero = not (Bytes.for_all (fun c -> c = '\x00') parity) in
112
+
Alcotest.(check bool) "nonzero parity" true has_nonzero
113
+
114
+
(* --- Decode no-error tests --- *)
115
+
116
+
let test_decode_no_errors () =
117
+
let data = Bytes.init config.k (fun i -> Char.chr (((i * 7) + 3) mod 256)) in
118
+
let codeword = Reed_solomon.encode_systematic config data in
119
+
let result = Reed_solomon.decode config codeword in
120
+
Alcotest.(check result_testable) "decode clean" (Ok data) result
121
+
122
+
(* --- Error correction tests --- *)
123
+
124
+
let test_roundtrip_1_error () =
125
+
let data = Bytes.init config.k (fun i -> Char.chr (i mod 256)) in
126
+
let codeword = Reed_solomon.encode_systematic config data in
127
+
let corrupted = Bytes.copy codeword in
128
+
inject_errors corrupted 1;
129
+
let result = Reed_solomon.decode config corrupted in
130
+
Alcotest.(check result_testable) "correct 1 error" (Ok data) result
131
+
132
+
let test_roundtrip_8_errors () =
133
+
let data = Bytes.init config.k (fun i -> Char.chr (((i * 3) + 17) mod 256)) in
134
+
let codeword = Reed_solomon.encode_systematic config data in
135
+
let corrupted = Bytes.copy codeword in
136
+
inject_errors corrupted 8;
137
+
let result = Reed_solomon.decode config corrupted in
138
+
Alcotest.(check result_testable) "correct 8 errors" (Ok data) result
139
+
140
+
let test_roundtrip_16_errors () =
141
+
let data = Bytes.init config.k (fun i -> Char.chr (((i * 11) + 5) mod 256)) in
142
+
let codeword = Reed_solomon.encode_systematic config data in
143
+
let corrupted = Bytes.copy codeword in
144
+
inject_errors corrupted 16;
145
+
let result = Reed_solomon.decode config corrupted in
146
+
Alcotest.(check result_testable) "correct 16 errors (max t)" (Ok data) result
147
+
148
+
let test_uncorrectable () =
149
+
let data = Bytes.init config.k (fun i -> Char.chr (i * 13 mod 256)) in
150
+
let codeword = Reed_solomon.encode_systematic config data in
151
+
let corrupted = Bytes.copy codeword in
152
+
inject_errors corrupted 17;
153
+
let result = Reed_solomon.decode config corrupted in
154
+
match result with
155
+
| Error _ -> () (* Any error variant is acceptable. *)
156
+
| Ok _ -> Alcotest.fail "should not decode with 17 errors"
157
+
158
+
(* --- Roundtrip with random data --- *)
159
+
160
+
let test_roundtrip_random () =
161
+
let rng = Random.State.make [| 12345 |] in
162
+
for trial = 0 to 9 do
163
+
let data =
164
+
Bytes.init config.k (fun _ -> Char.chr (Random.State.int rng 256))
165
+
in
166
+
let codeword = Reed_solomon.encode_systematic config data in
167
+
let corrupted = Bytes.copy codeword in
168
+
let num_errors = 1 + Random.State.int rng config.t in
169
+
(* Manual error injection with this trial's RNG. *)
170
+
let positions = Array.make config.n false in
171
+
let injected = ref 0 in
172
+
while !injected < num_errors do
173
+
let pos = Random.State.int rng config.n in
174
+
if not positions.(pos) then begin
175
+
positions.(pos) <- true;
176
+
let error_val = 1 + Random.State.int rng 255 in
177
+
let old = Char.code (Bytes.get corrupted pos) in
178
+
Bytes.set corrupted pos (Char.chr (old lxor error_val));
179
+
incr injected
180
+
end
181
+
done;
182
+
let result = Reed_solomon.decode config corrupted in
183
+
Alcotest.(check result_testable)
184
+
(Fmt.str "random trial %d (%d errors)" trial num_errors)
185
+
(Ok data) result
186
+
done
187
+
188
+
(* --- Invalid input tests --- *)
189
+
190
+
let test_invalid_encode_length () =
191
+
let data = Bytes.make 10 '\x00' in
192
+
Alcotest.check_raises "wrong length"
193
+
(Invalid_argument
194
+
"Reed_solomon.encode_systematic: expected 223 bytes, got 10") (fun () ->
195
+
ignore (Reed_solomon.encode_systematic config data))
196
+
197
+
let test_invalid_decode_length () =
198
+
let codeword = Bytes.make 100 '\x00' in
199
+
let result = Reed_solomon.decode config codeword in
200
+
Alcotest.(check result_testable)
201
+
"wrong codeword length"
202
+
(Error (`Invalid_length 100))
203
+
result
204
+
205
+
(* --- Suite --- *)
206
+
207
+
let suite =
208
+
( "reed-solomon",
209
+
[
210
+
Alcotest.test_case "gf256: add" `Quick test_gf_add;
211
+
Alcotest.test_case "gf256: mul" `Quick test_gf_mul;
212
+
Alcotest.test_case "gf256: inv" `Quick test_gf_inv;
213
+
Alcotest.test_case "gf256: pow" `Quick test_gf_pow;
214
+
Alcotest.test_case "gf256: exp/log" `Quick test_gf_exp_log_consistency;
215
+
Alcotest.test_case "encode: zeros" `Quick test_encode_zeros;
216
+
Alcotest.test_case "encode: length" `Quick test_encode_length;
217
+
Alcotest.test_case "encode: systematic" `Quick
218
+
test_encode_systematic_data_unchanged;
219
+
Alcotest.test_case "encode: nonzero parity" `Quick
220
+
test_encode_nonzero_parity;
221
+
Alcotest.test_case "decode: no errors" `Quick test_decode_no_errors;
222
+
Alcotest.test_case "decode: 1 error" `Quick test_roundtrip_1_error;
223
+
Alcotest.test_case "decode: 8 errors" `Quick test_roundtrip_8_errors;
224
+
Alcotest.test_case "decode: 16 errors" `Quick test_roundtrip_16_errors;
225
+
Alcotest.test_case "decode: uncorrectable" `Quick test_uncorrectable;
226
+
Alcotest.test_case "decode: random roundtrip" `Quick test_roundtrip_random;
227
+
Alcotest.test_case "encode: invalid length" `Quick
228
+
test_invalid_encode_length;
229
+
Alcotest.test_case "decode: invalid length" `Quick
230
+
test_invalid_decode_length;
231
+
] )
+7
test/test_reed_solomon.mli
+7
test/test_reed_solomon.mli
···
1
+
(** Reed-Solomon encoder/decoder tests.
2
+
3
+
Test vectors for CCSDS RS(255,223) including GF(2^8) arithmetic checks,
4
+
systematic encoding, error injection and correction up to t=16 errors, and
5
+
rejection of uncorrectable codewords. *)
6
+
7
+
val suite : string * unit Alcotest.test_case list