Reed-solomon CODEC

In [1]:

    

#require "reedsolomon,iocaml.notebook"


    

    




	Camlp4 Parsing version 4.01.0







    




/home/andyman/.opam/4.01.0-test/lib/reedsolomon: added to search path
/home/andyman/.opam/4.01.0-test/lib/reedsolomon/reedsolomon.cma: loaded
/home/andyman/.opam/4.01.0-test/lib/iocaml: added to search path

In [2]:

    

module X = Rs.BBCTest


    

    
Out[2]:




module X :
  sig
    module Gp : sig val pp : int val pe : int end
    module G :
      sig
        type t = int
        val zero : t
        val one : t
        val ( +: ) : t -> t -> t
        val ( -: ) : t -> t -> t
        val to_string : t -> string
        val alpha : t
        val n_elems : int
        val log : t -> int
        val antilog : int -> t
        val ( *: ) : t -> t -> t
        val ( /: ) : t -> t -> t
        val ( **: ) : t -> int -> t
        val inv : t -> t
      end
    module Rp : sig val k : int val t : int val b : int end
    module R :
      sig
        type elt = int
        module M :
          sig
            type t = elt
            type matrix = t array array
            val rows : matrix -> int
            val cols : matrix -> int
            val init : int -> int -> (int -> int -> t) -> matrix
            val create : int -> int -> matrix
            val copy : matrix -> matrix
            val identity : int -> matrix
            val transpose : matrix -> matrix
            val map : (t -> t) -> matrix -> matrix
            val map2 : (t -> t -> t) -> matrix -> matrix -> matrix
            val row_vector : t array -> matrix
            val col_vector : t array -> matrix
            val ( >> ) : matrix -> matrix -> matrix
            val ( ^^ ) : matrix -> matrix -> matrix
            val sub : int -> int -> int -> int -> matrix -> matrix
            val ( +: ) : matrix -> matrix -> matrix
            val ( -: ) : matrix -> matrix -> matrix
            val ( *: ) : matrix -> matrix -> matrix
            val ( *:. ) : matrix -> t -> matrix
            val minor : int -> int -> matrix -> matrix
            val det : matrix -> t
            val adjoint_inverse : matrix -> t * matrix
            val gauss_jordan : matrix -> matrix
            val gauss_jordan_inverse : matrix -> matrix
            module Row :
              sig
                val swap : int -> int -> int -> matrix
                val mult : int -> int -> t -> matrix
                val madd : int -> int -> int -> t -> matrix
              end
          end
        module R :
          sig
            type elt = elt
            type t = elt array
            val degree : t -> int
            val zero : t
            val one : t
            val x : t
            val to_poly : elt array -> t
            val of_poly : t -> elt array
            val copy : t -> t
            type poly_format =
              Rs.BBCTest.R.R.poly_format = {
              coef : elt -> string;
              indet : int -> string;
            }
            val poly_format : poly_format
            val string_format : bool -> poly_format -> int -> elt -> string
            val to_string :
              ?down:bool -> ?str:(int -> elt -> string) -> t -> string
            val trim : t -> t
            val slice : t -> int -> t
            val ( +: ) : t -> t -> t
            val ( -: ) : t -> t -> t
            val ( *: ) : t -> t -> t
            val ( /: ) : t -> t -> t * t
            val ( *:. ) : t -> elt -> t
            val ( /:. ) : t -> elt -> t
            val ( ^: ) : t -> int -> t
            val ( **: ) : t -> int -> t
            val ext_gcd : t -> t -> t * t
            val eval : t -> elt -> elt
          end
        type poly = R.t
        type loc = int
        val root : int -> elt
        val generator : poly
        val xn : int -> poly
        val x2t : poly
        val parity : poly -> poly
        val encode : poly -> poly
        val horner : poly -> elt -> elt
        val syndromes : poly -> poly
        val key_equations : poly -> int -> M.matrix * M.matrix
        val solve_key_equations : M.matrix * M.matrix -> M.matrix
        val peterson : poly -> poly
        val euclid_inner : poly * poly -> poly * poly -> poly * poly
        val euclid : ?norm:bool -> ?lim:int -> poly -> poly * poly
        val berlekamp_massey_iter :
          poly -> int -> poly * poly * int -> poly * poly * int
        val berlekamp_massey : poly -> poly
        module Sarwate :
          sig
            val iBM : poly -> poly
            val riBM : poly -> poly * poly
            val rriBM : poly -> poly * poly
            val forney : poly -> poly -> loc -> elt
          end
        val chien : poly -> loc list
        val error_location : loc -> int
        val error_magnitude : int -> poly -> poly -> poly
        val deriv : poly -> poly
        val forney : poly -> poly -> loc -> elt
        val error : elt list -> loc list -> poly
        val correct : poly -> poly -> poly
        val decode_euclid : poly -> poly
        val decode_berlekamp_massey : poly -> poly
        val decode_peterson : poly -> poly
        val decode : poly -> poly
        val erasure_locator : int list -> poly
        val zero_erasures : poly -> int list -> poly
        val error_and_erasure :
          elt list -> loc list -> elt list -> loc list -> poly
        val decode_erasures_euclid : poly -> int list -> poly
        val decode_erasures : poly -> int list -> poly
        val decode_errors_and_erasures_euclid : poly -> int list -> poly
        val decode_errors_and_erasures_berlekamp_massey :
          poly -> int list -> poly
        val decode_errors_and_erasures : poly -> int list -> poly
      end
  end

In [3]:

    

let pp, pe = X.Gp.pp, X.Gp.pe


    

    
Out[3]:




val pp : int = 19
val pe : int = 2

In [4]:

    

let k, t = X.Rp.k, X.Rp.t
let n = k + (2*t)


    

    
Out[4]:




val k : int = 11
val t : int = 2







    
Out[4]:




val n : int = 15

In [5]:

    

let number_of_field_elements = X.G.n_elems


    

    
Out[5]:




val number_of_field_elements : int = 16

In [6]:

    

(* utility function to display a html table *)
let table rows cols f =
    let tag x s = "<" ^ x ^ ">" ^ s ^ "</" ^ x ^ ">" in
    let concat x s = s |> Array.map (tag x) |> Array.to_list |> String.concat "" in
    let init n x f = Array.init n f |> concat x in
    Iocaml.display "text/html" 
        ("<table>" ^ (init rows "tr" (fun row -> init cols "td" (f row))) ^ "</table>")


    

    
Out[6]:




val table : int -> int -> (int -> int -> string) -> unit = <fun>

In [7]:

    

table 16 16 (fun i j -> X.G.(to_string (i +: j)))


    

    





0
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
1
0
3
2
5
4
7
6
9
8
11
10
13
12
15
14
2
3
0
1
6
7
4
5
10
11
8
9
14
15
12
13
3
2
1
0
7
6
5
4
11
10
9
8
15
14
13
12
4
5
6
7
0
1
2
3
12
13
14
15
8
9
10
11
5
4
7
6
1
0
3
2
13
12
15
14
9
8
11
10
6
7
4
5
2
3
0
1
14
15
12
13
10
11
8
9
7
6
5
4
3
2
1
0
15
14
13
12
11
10
9
8
8
9
10
11
12
13
14
15
0
1
2
3
4
5
6
7
9
8
11
10
13
12
15
14
1
0
3
2
5
4
7
6
10
11
8
9
14
15
12
13
2
3
0
1
6
7
4
5
11
10
9
8
15
14
13
12
3
2
1
0
7
6
5
4
12
13
14
15
8
9
10
11
4
5
6
7
0
1
2
3
13
12
15
14
9
8
11
10
5
4
7
6
1
0
3
2
14
15
12
13
10
11
8
9
6
7
4
5
2
3
0
1
15
14
13
12
11
10
9
8
7
6
5
4
3
2
1
0






    
Out[7]:




- : unit = ()

In [8]:

    

table 16 16 (fun i j -> X.G.(to_string (i *: j)))


    

    





0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
0
2
4
6
8
10
12
14
3
1
7
5
11
9
15
13
0
3
6
5
12
15
10
9
11
8
13
14
7
4
1
2
0
4
8
12
3
7
11
15
6
2
14
10
5
1
13
9
0
5
10
15
7
2
13
8
14
11
4
1
9
12
3
6
0
6
12
10
11
13
7
1
5
3
9
15
14
8
2
4
0
7
14
9
15
8
1
6
13
10
3
4
2
5
12
11
0
8
3
11
6
14
5
13
12
4
15
7
10
2
9
1
0
9
1
8
2
11
3
10
4
13
5
12
6
15
7
14
0
10
7
13
14
4
9
3
15
5
8
2
1
11
6
12
0
11
5
14
10
1
15
4
7
12
2
9
13
6
8
3
0
12
11
7
5
9
14
2
10
6
1
13
15
3
4
8
0
13
9
4
1
12
8
5
2
15
11
6
3
14
10
7
0
14
15
1
13
3
2
12
9
7
6
8
4
10
11
5
0
15
13
2
9
6
4
11
1
14
12
3
8
7
5
10






    
Out[8]:




- : unit = ()

In [9]:

    

table 1 16 (fun _ i -> X.G.(to_string (inv i)))


    

    





0
1
9
14
13
11
7
6
15
2
12
5
10
4
3
8






    
Out[9]:




- : unit = ()

In [10]:

    

module IPoly = Poly.Make(Ops.Int)


    

    
Out[10]:




module IPoly :
  sig
    type elt = Ops.Int.t
    type t = Ops.Int.t array
    val degree : t -> int
    val zero : t
    val one : t
    val x : t
    val to_poly : elt array -> t
    val of_poly : t -> elt array
    val copy : t -> t
    type poly_format =
      Poly.Make(Ops.Int).poly_format = {
      coef : elt -> string;
      indet : int -> string;
    }
    val poly_format : poly_format
    val string_format : bool -> poly_format -> int -> elt -> string
    val to_string : ?down:bool -> ?str:(int -> elt -> string) -> t -> string
    val trim : t -> t
    val slice : t -> int -> t
    val ( +: ) : t -> t -> t
    val ( -: ) : t -> t -> t
    val ( *: ) : t -> t -> t
    val ( /: ) : t -> t -> t * t
    val ( *:. ) : t -> elt -> t
    val ( /:. ) : t -> elt -> t
    val ( ^: ) : t -> int -> t
    val ( **: ) : t -> int -> t
    val ext_gcd : t -> t -> t * t
    val eval : t -> elt -> elt
  end

In [11]:

    

(* display poly *)
let display_poly poly = 
    (* modifiy printing a little to suit mathjax *)
    let str = X.R.R.(string_format true 
        { poly_format with 
            indet = function 0 -> "" | 1 -> "x" | _ as n -> "x^{" ^ string_of_int n ^ "}"; }) in
    Iocaml.display "text/html" ("$$" ^ X.R.R.to_string ~str:str poly ^ "$$")


    

    
Out[11]:




val display_poly : X.R.R.t -> unit = <fun>

In [12]:

    

let f,g = [| 7; 3 |], [| 4; 0; 2 |]


    

    
Out[12]:




val f : int array = [|7; 3|]
val g : int array = [|4; 0; 2|]

In [13]:

    

display_poly f;;
display_poly g;;


    

    





$$3.x + 7$$






    





$$2.x^{2} + 4$$






    
Out[13]:




- : unit = ()







    
Out[13]:




- : unit = ()

In [14]:

    

display_poly IPoly.(f *: g)


    

    





$$6.x^{3} + 14.x^{2} + 12.x + 28$$






    
Out[14]:




- : unit = ()

In [15]:

    

table 2 2 (fun i j -> Galois.Primitive.GF2.(to_string (i +: j)));;
table 2 2 (fun i j -> Galois.Primitive.GF2.(to_string (i *: j)));;


    

    





0
1
1
0






    





0
0
0
1






    
Out[15]:




- : unit = ()







    
Out[15]:




- : unit = ()

In [16]:

    

module GFPoly = Galois.Extension.Make(struct 
    module Poly = Poly.Make(Galois.Primitive.GF2)
    let pp = [| 1; 1; 0; 0; 1 |]
end)


    

    
Out[16]:




module GFPoly :
  sig
    type t = Galois.Primitive.GF2.t array
    val zero : t
    val one : t
    val ( +: ) : t -> t -> t
    val ( -: ) : t -> t -> t
    val ( *: ) : t -> t -> t
    val ( /: ) : t -> t -> t
    val to_string : t -> string
  end

In [17]:

    

let x : GFPoly.t = [|0; 1; 0; 0|]


    

    
Out[17]:




val x : GFPoly.t = [|0; 1; 0; 0|]

In [18]:

    

GFPoly.to_string x


    

    
Out[18]:




- : string = "x"

In [19]:

    

GFPoly.(to_string (x *: x))


    

    
Out[19]:




- : string = "x^2"

In [20]:

    

let rec pow = GFPoly.(function
    | 0 -> one
    | 1 -> x
    | _ as n -> x *: pow (n-1))


    

    
Out[20]:




val pow : int -> GFPoly.t = <fun>

In [21]:

    

table 15 3 (fun i -> function 
    | 0 -> "$x^{" ^ string_of_int i ^ "}$"
    | 1 -> "$" ^ GFPoly.to_string (pow i) ^ "$"
    | _ -> X.G.(to_string (X.Gp.pe **: i)))


    

    





$x^{0}$
$1$
1
$x^{1}$
$x$
2
$x^{2}$
$x^2$
4
$x^{3}$
$x^3$
8
$x^{4}$
$x + 1$
3
$x^{5}$
$x^2 + x$
6
$x^{6}$
$x^3 + x^2$
12
$x^{7}$
$x^3 + x + 1$
11
$x^{8}$
$x^2 + 1$
5
$x^{9}$
$x^3 + x$
10
$x^{10}$
$x^2 + x + 1$
7
$x^{11}$
$x^3 + x^2 + x$
14
$x^{12}$
$x^3 + x^2 + x + 1$
15
$x^{13}$
$x^3 + x^2 + 1$
13
$x^{14}$
$x^3 + 1$
9






    
Out[21]:




- : unit = ()

In [22]:

    

let message = [| 11; 10; 9; 8; 7; 6; 5; 4; 3; 2; 1 |]
let () = display_poly message


    

    





$$x^{10} + 2.x^{9} + 3.x^{8} + 4.x^{7} + 5.x^{6} + 6.x^{5} + 7.x^{4} + 8.x^{3} + 9.x^{2} + 10.x + 11$$






    
Out[22]:




val message : int array = [|11; 10; 9; 8; 7; 6; 5; 4; 3; 2; 1|]

In [23]:

    

let codeword = X.R.encode message
let () = display_poly codeword


    

    





$$x^{14} + 2.x^{13} + 3.x^{12} + 4.x^{11} + 5.x^{10} + 6.x^{9} + 7.x^{8} + 8.x^{7} + 9.x^{6} + 10.x^{5} + 11.x^{4} + 3.x^{3} + 3.x^{2} + 12.x + 12$$






    
Out[23]:




val codeword : X.R.poly = [|12; 12; 3; 3; 11; 10; 9; 8; 7; 6; 5; 4; 3; 2; 1|]

In [24]:

    

let myencode message = 
    let open X.R.R in
    let (%:) a b = snd (a /: b) in (* mod; the remainder of division *)
    let message = message *: (x **: (2*t)) in (* shift message *)
    message +: (message %: X.R.generator)


    

    
Out[24]:




val myencode : X.R.R.t -> X.R.R.t = <fun>

In [25]:

    

let mycodeword = myencode message
let () = display_poly mycodeword


    

    





$$x^{14} + 2.x^{13} + 3.x^{12} + 4.x^{11} + 5.x^{10} + 6.x^{9} + 7.x^{8} + 8.x^{7} + 9.x^{6} + 10.x^{5} + 11.x^{4} + 3.x^{3} + 3.x^{2} + 12.x + 12$$






    
Out[25]:




val mycodeword : X.R.R.t =
  [|12; 12; 3; 3; 11; 10; 9; 8; 7; 6; 5; 4; 3; 2; 1|]

In [26]:

    

let error = [| 0; 0; 2; 0; 0; 0; 0; 0; 0; 13; 0; 0; 0; 0; 0 |]
let () = display_poly error


    

    





$$13.x^{9} + 2.x^{2}$$






    
Out[26]:




val error : int array = [|0; 0; 2; 0; 0; 0; 0; 0; 0; 13; 0; 0; 0; 0; 0|]

In [27]:

    

let received = X.R.R.(codeword +: error)
let () = display_poly received


    

    





$$x^{14} + 2.x^{13} + 3.x^{12} + 4.x^{11} + 5.x^{10} + 11.x^{9} + 7.x^{8} + 8.x^{7} + 9.x^{6} + 10.x^{5} + 11.x^{4} + 3.x^{3} + x^{2} + 12.x + 12$$






    
Out[27]:




val received : X.R.R.t = [|12; 12; 1; 3; 11; 10; 9; 8; 7; 11; 5; 4; 3; 2; 1|]

In [28]:

    

let corrected = X.R.decode received
let () = display_poly corrected
let okay = corrected = codeword


    

    





$$x^{14} + 2.x^{13} + 3.x^{12} + 4.x^{11} + 5.x^{10} + 6.x^{9} + 7.x^{8} + 8.x^{7} + 9.x^{6} + 10.x^{5} + 11.x^{4} + 3.x^{3} + 3.x^{2} + 12.x + 12$$






    
Out[28]:




val corrected : X.R.poly =
  [|12; 12; 3; 3; 11; 10; 9; 8; 7; 6; 5; 4; 3; 2; 1|]







    
Out[28]:




val okay : bool = true

In [29]:

    

let syndromes = Array.init (2*t) X.R.(fun i -> horner received (root i))


    

    
Out[29]:




val syndromes : X.R.elt array = [|15; 3; 4; 12|]

In [30]:

    

let v, l = X.R.euclid syndromes


    

    
Out[30]:




val v : X.R.poly = [|14; 3|]
val l : X.R.poly = [|9; 7; 7|]

In [31]:

    

let el = X.R.chien l
let el' = List.map X.R.error_location el (* compute actual locations *)


    

    
Out[31]:




val el : X.R.loc list = [6; 13]







    
Out[31]:




val el' : int list = [9; 2]

In [32]:

    

let ev = List.map (X.R.forney v l) el


    

    
Out[32]:




val ev : X.R.elt list = [13; 2]

In [33]:

    

let e = X.R.error ev el


    

    
Out[33]:




val e : X.R.poly = [|0; 0; 2; 0; 0; 0; 0; 0; 0; 13|]

In [34]:

    

let corrected = X.R.R.(received +: e)
let () = display_poly corrected


    

    





$$x^{14} + 2.x^{13} + 3.x^{12} + 4.x^{11} + 5.x^{10} + 6.x^{9} + 7.x^{8} + 8.x^{7} + 9.x^{6} + 10.x^{5} + 11.x^{4} + 3.x^{3} + 3.x^{2} + 12.x + 12$$






    
Out[34]:




val corrected : X.R.R.t = [|12; 12; 3; 3; 11; 10; 9; 8; 7; 6; 5; 4; 3; 2; 1|]