(* Field *) (* $Id: field.ml,v 1.1 2003/01/12 15:33:00 berke Exp $ *) (* By Berke Durak *) open Printf (* Operations in GF(2). These are really trivial. *) type poly_gf2 = bool array let zero_gf2 = false let one_gf2 = true let is_zero_gf2 = not let is_one_gf2 x = x let add_gf2 x y = if x then not y else y let sub_gf2 = add_gf2 let neg_gf2 x = x let mul_gf2 x y = if x then y else x let inv_gf2 x = if not x then raise Division_by_zero else x let div_gf2 x y = if not y then raise Division_by_zero else x let next_gf2 x = not x (* Polynomial operations over an arbitrary *field*. *) let is_poly_zero_gf2 p = if Array.length p = 0 then true else try begin for i = 0 to Array.length p - 1 do if not (is_zero_gf2 p.(i)) then raise Not_found; done; true end with Not_found -> false | x -> raise x let print_coef_gf2 c = if c then print_string "1" else print_string "0" let print_monome_gf2 d c = if c then match d with 0 -> print_coef_gf2 c | 1 -> begin if not (is_one_gf2 c) then print_coef_gf2 c; print_string "X" end | x -> begin if not (is_one_gf2 c) then print_coef_gf2 c; print_string "X^"; print_int x end else print_coef_gf2 zero_gf2 let print_poly_gf2 p = let l = Array.length p and f = ref false in begin for i = 0 to l - 1 do if not (is_zero_gf2 p.(i)) then begin if !f then print_string " + " else f := true; print_monome_gf2 i p.(i) end done; if not !f then print_coef_gf2 zero_gf2 end type poly_degree = Minus_infinity | Poly_degree of int;; let degree_poly_gf2 p = let rec loop i = if i < 0 then Minus_infinity else if not (is_zero_gf2 p.(i)) then Poly_degree i else loop (i - 1) in loop (Array.length p - 1) let zero_poly_gf2 = [| |] let rec add_poly_gf2 p1 p2 = let d1 = degree_poly_gf2 p1 and d2 = degree_poly_gf2 p2 in match d1,d2 with Minus_infinity, _ -> p2 | _, Minus_infinity -> p1 | Poly_degree d1, Poly_degree d2 -> if d2 < d1 then add_poly_gf2 p2 p1 else let d3 = max d1 d2 in let p3 = Array.make (d3 + 1) false in begin for i = 0 to d1 do p3.(i) <- add_gf2 p1.(i) p2.(i) done; for i = (d1 + 1) to d2 do p3.(i) <- p2.(i) done; p3 end let are_equal_poly_gf2 p1 p2 = let d1 = degree_poly_gf2 p1 and d2 = degree_poly_gf2 p2 in match d1, d2 with Minus_infinity, Minus_infinity -> true | Minus_infinity, _ -> false | _, Minus_infinity -> false | Poly_degree d1, Poly_degree d2 -> if d1 = d2 then let rec loop i = if i > d1 then true else if p1.(i) != p2.(i) then false else loop (i + 1) in loop 0 else false let rec sub_poly_gf2 p1 p2 = let d1 = degree_poly_gf2 p1 and d2 = degree_poly_gf2 p2 in match d1,d2 with Minus_infinity, _ -> p2 | _, Minus_infinity -> p1 | Poly_degree d1, Poly_degree d2 -> if d2 < d1 then add_poly_gf2 p2 p1 else let d3 = max d1 d2 in let p3 = Array.make (d3 + 1) false in begin for i = 0 to d1 do p3.(i) <- sub_gf2 p1.(i) p2.(i) done; for i = (d1 + 1) to d2 do p3.(i) <- p2.(i) done; p3 end let set_coef_poly_gf2 p1 i c = let n = Array.length p1 in if i >= n then let p2 = Array.create (i + i / 4 + 1) false in begin for i = 0 to n - 1 do p2.(i) <- p1.(i) done; p2.(i) <- c; p2 end else begin p1.(i) <- c; p1 end let get_coef_poly_gf2 p1 i = let n = Array.length p1 in if i >= n then zero_gf2 else p1.(i) let next_poly_gf2 p = let rec loop p i = let c = get_coef_poly_gf2 p i in let c_next = next_gf2 c in let p = set_coef_poly_gf2 p i c_next in if is_zero_gf2 c_next then loop p (i + 1) else p in loop p 0 let print_poly_compact_gf2 p n = for i = 0 to n do print_coef_gf2 (get_coef_poly_gf2 p i) done (* multiply p2 by cX^n *) let rec shift_mul_poly_gf2 p1 n c = if is_zero_gf2 c then zero_poly_gf2 else match degree_poly_gf2 p1 with Minus_infinity -> zero_poly_gf2 | Poly_degree d1 -> let p2 = Array.create (d1 + n + 1) false in begin for i = 0 to d1 do p2.(i + n) <- mul_gf2 c p1.(i) done; p2 end let rec mul_poly_gf2 p1 p2 = let d1 = degree_poly_gf2 p1 and d2 = degree_poly_gf2 p2 in match d1,d2 with Minus_infinity, _ -> zero_poly_gf2 | _, Minus_infinity -> zero_poly_gf2 | Poly_degree d1, Poly_degree d2 -> if d2 < d1 then mul_poly_gf2 p2 p1 else let d3 = d1 + d2 in let p3 = Array.make (d3 + 1) false in begin for i = 0 to d1 do for j = 0 to d2 do p3.(i + j) <- add_gf2 p3.(i + j) (mul_gf2 p1.(i) p2.(j)) done done; p3 end (* compute p1 div p2, p1 Mod p2 *) let rec div_poly_gf2 p1 p2 = let d1 = degree_poly_gf2 p1 and d2 = degree_poly_gf2 p2 in match d1,d2 with Minus_infinity, _ -> (zero_poly_gf2, zero_poly_gf2) | _, Minus_infinity -> raise Division_by_zero | Poly_degree d1, Poly_degree d2 -> (* p3 = p1 div p2, of degree deg(p1) - deg(p2) *) (* p4 = p1 Mod p2 of degree < deg(p2) *) if d1 < d2 then (zero_poly_gf2, p1) else (* d1 <= d2 *) (* d1 is the degree of p1 *) (* p1 is the current partial remainder *) (* p3 is the current partial quotient *) let rec loop p1 d1 p3 = if d1 < d2 then (p3, p1) else let c = div_gf2 p2.(d2) p1.(d1) in let p4 = shift_mul_poly_gf2 p2 (d1 - d2) c in let p1 = sub_poly_gf2 p1 p4 in let p3 = set_coef_poly_gf2 p3 (d1 - d2) c in match degree_poly_gf2 p1 with Minus_infinity -> (p3, p1) | Poly_degree d1 -> loop p1 d1 p3 in loop p1 d1 zero_poly_gf2 (* X^3 *) (* 1+X^4+X^7+X^8 *) let int_of_poly_gf2 p = match degree_poly_gf2 p with Poly_degree d -> begin let rec loop i x = if i < 0 then x else loop (i - 1) (x * 2 + (if get_coef_poly_gf2 p i then 1 else 0)) in loop d 0 end | Minus_infinity -> 0 let check_div_poly p1 p2 = let q, r = div_poly_gf2 p1 p2 in begin print_string "("; print_poly_gf2 q; print_string ","; print_poly_gf2 r; print_string "); q * p2 + r = "; print_poly_gf2 (add_poly_gf2 (mul_poly_gf2 p2 q) r); print_string "\n"; end let zozo x = let p1 = [| false; false; false; true |] and p2 = [| true; false; false; false; true; false; false; true; true |] in begin print_string "p1 = "; print_poly_gf2 p1; print_string "\n"; print_string "p2 = "; print_poly_gf2 p2; print_string "\n"; print_string "p1 + p2 = "; print_poly_gf2 (add_poly_gf2 p1 p2); print_string "\n"; print_string "p1 * p2 = "; print_poly_gf2 (mul_poly_gf2 p1 p2); print_string "\n"; print_string "(p1 div p2, p1 Mod p2) = "; check_div_poly p1 p2; print_string "(p2 div p1, p2 Mod p1) = "; check_div_poly p2 p1; end;; let string_to_poly_gf2 s = let l = String.length s and p = ref zero_poly_gf2 in begin for i = 0 to l - 1 do let c = s.[i] in if c = '1' then p := set_coef_poly_gf2 !p i one_gf2 done; !p end exception Not_primitive let display_generator_table_gf2 pm g = begin print_string "/* table generator = "; print_poly_gf2 g; print_string " primitive = "; print_poly_gf2 pm; print_string "*/\n\n"; print_string "unsigned char table[] = {"; match degree_poly_gf2 pm with Poly_degree dm -> let x = ref (set_coef_poly_gf2 zero_poly_gf2 0 one_gf2) and order = ref ((1 lsl dm) - 1) in begin for j = 0 to (1 lsl dm) - 2 do let (q,r) = div_poly_gf2 !x pm in begin if j mod 10 = 0 then print_string "\n\t"; if degree_poly_gf2 r = Minus_infinity then raise Not_primitive; printf "0x%02x, " (int_of_poly_gf2 r); x := mul_poly_gf2 !x g end done; print_string "};\n"; end | Minus_infinity -> raise Division_by_zero end ;; let display_generators_table_gf2 pm = begin print_string "Table for "; print_poly_gf2 pm; print_string "\n"; match degree_poly_gf2 pm with Poly_degree dm -> let p = ref (next_poly_gf2 zero_poly_gf2) and one = next_poly_gf2 zero_poly_gf2 in for i = 1 to (1 lsl dm) - 1 do begin if 0 = i land 0x3 then begin prerr_string "Processed "; prerr_int i; prerr_newline () end; print_string "\t["; print_poly_compact_gf2 !p (dm - 1); print_string "] "; let g = ref (set_coef_poly_gf2 zero_poly_gf2 0 one_gf2) and order = ref ((1 lsl dm) - 1) in begin try for j = 0 to (1 lsl dm) - 1 do let (q,r) = div_poly_gf2 !g pm in begin if degree_poly_gf2 r = Minus_infinity then raise Not_primitive; print_poly_compact_gf2 r (dm - 1); print_string " "; if j > 0 && are_equal_poly_gf2 r one then begin order := j; raise Not_found end; g := mul_poly_gf2 !g !p end done; raise Not_found with Not_found -> begin print_string " "; print_int !order; print_string "\n" end | x -> raise x end; p := next_poly_gf2 !p end done | _ -> raise Division_by_zero end;; (* display_generators_table_gf2 [| true; true; false; false; true |];; *) (* 1 + X^4 *) (* zozo 1 *) (* let main argv = let pm = if Array.length argv <= 1 then begin print_string "Enter polynomial over GF(2) in compact notation: "; flush stdout; let s = read_line () in string_to_poly_gf2 s end else string_to_poly_gf2 argv.(1) in begin print_poly_gf2 pm; print_string "\n"; display_generators_table_gf2 pm; end;; *) let main argv = let (pm, g) = if Array.length argv <= 2 then begin print_string "Enter primitive polynomial over GF(2) in compact notation: "; flush stdout; let s1 = read_line () in begin print_string "Enter generator polynomial: "; flush stdout; let s2 = read_line () in (string_to_poly_gf2 s1, string_to_poly_gf2 s2) end end else (string_to_poly_gf2 argv.(1), string_to_poly_gf2 argv.(2)) in begin print_poly_gf2 pm; print_string "\n"; display_generator_table_gf2 pm g; end;; main Sys.argv;;