The current state of the program computes the number of common factors.
--- /dev/null
+OASISFormat: 0.4
+Name: lu
+Version: 0.1
+Synopsis: Simple LDU implementation in OCaml.
+Authors: Johannes Middeke
+License: GPL
+Plugins: META (0.4)
+
+Executable "lu"
+ Path: .
+ BuildTools: ocamlbuild
+ MainIs: test.ml
+ CompiledObject: best
+ BuildDepends: zarith, tools
--- /dev/null
+let data = [
+10,9.50000,3.46667;
+10,9.96667,3.86667;
+10,8.43333,3.30000;
+10,10.36667,4.16667;
+10,10.20000,4.06667;
+10,10.03333,3.90000;
+10,7.63333,3.46667;
+10,9.50000,3.20000;
+10,7.93333,3.56667;
+10,9.43333,3.53333;
+15,13.36667,4.86667;
+15,14.43333,5.70000;
+15,14.40000,5.50000;
+15,14.93333,5.66667;
+15,16.76667,5.73333;
+15,15.26667,5.93333;
+15,15.06667,5.03333;
+15,17.26667,6.26667;
+15,15.86667,6.23333;
+15,14.76667,5.96667;
+20,19.73333,8.30000;
+20,21.60000,8.20000;
+20,20.00000,7.30000;
+20,20.80000,8.76667;
+20,23.23333,8.13333;
+20,21.60000,7.86667;
+20,23.50000,7.70000;
+20,23.36667,8.26667;
+20,20.96667,8.06667;
+20,23.40000,9.33333;
+25,26.73333,9.93333;
+25,30.40000,10.76667;
+25,27.33333,10.16667;
+25,29.80000,11.76667;
+25,28.83333,10.40000;
+25,28.06667,9.93333;
+25,28.50000,9.96667;
+25,28.36667,9.50000;
+25,30.70000,10.66667;
+25,29.83333,11.33333;
+30,31.83333,12.10000;
+30,36.26667,12.16667;
+30,33.30000,12.33333;
+30,31.53333,12.00000;
+30,36.16667,13.06667;
+30,37.43333,13.60000;
+30,35.70000,12.93333;
+30,35.73333,12.66667;
+30,31.53333,12.10000;
+30,36.73333,14.20000;
+35,42.43333,16.20000;
+35,43.90000,15.70000;
+35,41.33333,15.03333;
+35,43.36667,16.53333;
+35,41.06667,13.96667;
+35,39.60000,14.13333;
+35,42.13333,14.66667;
+35,38.13333,15.43333;
+35,42.90000,14.70000;
+35,42.06667,15.26667;
+5,3.93333,1.66667;
+5,3.13333,1.53333;
+5,3.43333,1.50000;
+5,2.76667,1.30000;
+5,4.00000,1.70000;
+5,3.06667,1.43333;
+5,3.93333,1.66667;
+5,4.03333,1.70000;
+5,3.30000,1.43333;
+5,3.46667,1.70000;
+]
--- /dev/null
+#use "topfind";;
+#require "tools";;
+#use "data.ml";;
+
+let data2 =
+ let replace (b,c) (x,y) = (b@x, c@y) in
+ List.map (fun (n,total,predict) -> (n, ([predict],[total]))) data
+ |> List.fold_left (ExtList.add_assoc ~replace) []
+ |> ExtList.sort_by ~compare fst
+
+let mean, deviation =
+ let mean (a, (b,c)) = (a, (Statistics.mean b, Statistics.mean c))
+ and dev (a, (b,c)) =
+ (a, (Statistics.standard_deviation b, Statistics.standard_deviation c))
+ in
+ List.map mean data2, List.map dev data2
+
+let ratio =
+ List.map (fun (a,(b,c)) -> (a, b/.c)) mean
+
+let mean_ratio, deviation_ratio =
+ let data = List.map snd ratio in
+ Statistics.mean data,
+ Statistics.standard_deviation data
+
+let () =
+ let (@@) f g x = f (g x) in
+ List.map (string_of_float @@ snd) ratio
+ |> String.concat ", "
+ |> Printf.printf "%%%% Ratios: %s.\n";
+ Printf.printf "\\gdef\\MeanRatio{%.2f}\n" (100. *. mean_ratio);
+ Printf.printf "\\gdef\\DeviationRatio{%.2f}\n"
+ (100. *. deviation_ratio)
+
+let tikz_mean =
+ let point n x =
+ Printf.sprintf "(%d,%.3f) node[circle,fill,inner sep=1pt] {}" n x
+ in
+ match mean with
+ | [] -> ()
+ | (n,(_p,a)) :: mean ->
+ begin
+ Printf.printf "\\draw[total] %s" (point n a);
+
+ ExtList.each mean (fun (n,(_p,a)) ->
+ Printf.printf "\n -- %s" (point n a));
+
+
+ Printf.printf "node[right] {total};\n"
+ end
+
+
+
+let tikz_mean_predicted =
+ let point n x =
+ Printf.sprintf "(%d,%.3f) node[circle,fill,inner sep=1pt] {}" n x
+ in
+ match mean with
+ | [] -> ()
+ | (n,(p,_a)) :: mean ->
+ begin
+ Printf.printf "\\draw[predicted] %s" (point n p);
+
+ ExtList.each mean (fun (n,(p,_a)) ->
+ Printf.printf "\n -- %s" (point n p));
+
+
+ Printf.printf "node[right] {predicted};\n"
+ end
--- /dev/null
+#!/bin/bash
+
+while true
+do
+ bash prep.sh
+ ocaml makeHTML.ml > Bareiss.html
+ scp Bareiss.html jmiddeke@www:/home/www/people/jmiddeke/
+ sleep 1m
+done
--- /dev/null
+type t = Z.t
+
+let eql = Z.equal
+
+let add = Z.add
+and sub = Z.sub
+and neg = Z.neg
+and nul = Z.zero
+
+let mul = Z.mul
+and div _ _ = failwith "Division not implemented."
+and inv _ = failwith "Inversion not implemented."
+and one = Z.one
+
+let exact_div = Z.divexact
+
+let show = Z.to_string
+
+let of_int = Z.of_int
+
+let size = Z.size
+
+let gcd a b =
+ if Z.equal a Z.zero
+ then b
+ else
+ if Z.equal b Z.zero
+ then a
+ else Z.gcd a b
+
+let rem = Z.rem
+
--- /dev/null
+#use "topfind";;
+#require "tools";;
+#use "data.ml";;
+
+
+
+(* The raw input data is a list consisting of tuples of the form
+ [(n,total,prediction)] where each [n] can occur several times. We
+ first need to group triples with the same [n] together creating a
+ list with entries [(n, totals, predictions)] where each [n] is
+ unique and [totals] and [predictions] are lists containing all the
+ values of [total] and [prediction] for this [n]. *)
+
+let data2 =
+ let replace (b,c) (x,y) = (b@x, c@y) in
+ List.map (fun (n,total,predict) -> (n, ([predict],[total]))) data
+ |> List.fold_left (ExtList.add_assoc ~replace) []
+ |> ExtList.sort_by ~compare fst
+
+(* From the list with triples [(n,totals,predictions)], we create two
+ new lists [mean] and [deviation] where [mean] contains triples
+ [(n,total_mean,prediction_mean)] such that [total_mean] is the mean
+ value of the list [totals] for this [n] and similar for the other
+ lists. *)
+
+let mean, deviation =
+ let mean (a, (b,c)) = (a, (Statistics.mean b, Statistics.mean c))
+ and dev (a, (b,c)) =
+ (a, (Statistics.standard_deviation b, Statistics.standard_deviation c))
+ in
+ List.map mean data2, List.map dev data2
+
+(* The next value tells us for each [n] how many factors we found. *)
+
+let ratio =
+ List.map (fun (a,(b,c)) -> (a, b/.c)) mean
+
+(* The total amount of factors found by the predicition for all
+ possible [n]. *)
+
+let mean_ratio, deviation_ratio =
+ let data = List.map snd ratio in
+ Statistics.mean data,
+ Statistics.standard_deviation data
+
+(* *)
+
+
+let svg_of_coord (x,y) = (5.*.x +. 25., 425. -. 5.*.y)
+
+let coord ff (x,y) =
+ let (x,y) = svg_of_coord (x,y) in
+ Printf.fprintf ff "%.2f,%.2f" x y
+
+let svg_poly_line colour marker = function
+ | [] -> ()
+ | x :: xs ->
+ Printf.printf "<path d=\"M%a" coord x;
+ List.iter (fun y -> Printf.printf " L%a" coord y) xs;
+ Printf.printf "\" style=\"stroke:%s; fill:none;\" " colour;
+ Printf.printf " marker-start=\"url(#circle-%s)\"" marker;
+ Printf.printf " marker-mid=\"url(#circle-%s)\"" marker;
+ Printf.printf " marker-end=\"url(#circle-%s)\" />" marker;
+ let z = ExtList.last (x :: xs) in
+ let (x,y) = svg_of_coord z in
+ Printf.printf "<text x=\"%.2f\" y=\"%.2f\" " (x +. 5.) (y +. 3.);
+ Printf.printf "fill=\"%s\" " colour;
+ Printf.printf "style=\"font-size:14px\" ";
+ Printf.printf "text-anchor=\"start\"> %s </text>" marker
+
+
+
+
+let x_tick x =
+ let (x1,y1) = svg_of_coord (float x, 0.)
+ and (x2,y2) = svg_of_coord (float x,-1.)
+ in
+ Printf.printf "<path d=\"M%.2f,%.2f L%.2f,%.2f\" " x1 y1 x2 y2;
+ Printf.printf "stroke=\"black\" stroke-width=\"1\" />\n";
+ Printf.printf "<text x=\"%.2f\" y=\"%.2f\" " x2 (y2 +. 12.);
+ Printf.printf "style=\"font-size:12px\" text-anchor=\"middle\">\n";
+ Printf.printf " %d\n" x;
+ Printf.printf "</text>\n"
+
+let rec x_ticks ?(by=5) from until =
+ if from <= until
+ then (x_tick from; x_ticks ~by (from+by) until)
+
+let y_tick y =
+ let (x1,y1) = svg_of_coord ( 0., float y)
+ and (x2,y2) = svg_of_coord (-1., float y)
+ in
+ Printf.printf "<path d=\"M%.2f,%.2f L%.2f,%.2f\" " x1 y1 x2 y2;
+ Printf.printf "stroke=\"black\" stroke-width=\"1\" />\n";
+ Printf.printf "<text x=\"%.2f\" y=\"%.2f\" " (x2 -. 1.) (y2 +. 5.);
+ Printf.printf "style=\"font-size:12px\" text-anchor=\"end\">\n";
+ Printf.printf " %d\n" y;
+ Printf.printf "</text>\n"
+
+let rec y_ticks ?(by=5) from until =
+ if from <= until
+ then (y_tick from; y_ticks ~by (from+by) until)
+
+
+let axes () =
+ Printf.printf "<path d=\"M25,25 L25,425 L425,425\"\n";
+ Printf.printf " stroke=\"#000\"\n";
+ Printf.printf " stroke-width=\"2\"\n";
+ Printf.printf " fill=\"none\"\n";
+ Printf.printf " marker-start=\"url(#arrow-begin)\" \n";
+ Printf.printf " marker-end=\"url(#arrow-end)\" \n";
+ Printf.printf " />\n";
+ Printf.printf "\n";
+ Printf.printf "<text x=\"435\" y=\"430\"";
+ Printf.printf " style=\"font-size:14px;font-style:italic;\"";
+ Printf.printf " text-anchor=\"start\">\n";
+ Printf.printf " n\n";
+ Printf.printf "</text>\n";
+ Printf.printf "<text x=\"35\" y=\"25\" style=\"font-size:14px\"";
+ Printf.printf " text-anchor=\"start\">\n";
+ Printf.printf " # factors\n";
+ Printf.printf "</text>\n";
+ Printf.printf "\n";
+ x_ticks ~by:5 5 75;
+ y_ticks ~by:5 5 75
+
+let head total predict =
+ Printf.printf "<!DOCTYPE html>\n";
+ Printf.printf "<html lang=\"en-CA\">\n";
+ Printf.printf "<head>\n";
+ Printf.printf "<title>Bareiss & Polynomials over GF(2)[x]</title>\n";
+ Printf.printf "\n";
+ Printf.printf "<link href='http://fonts.googleapis.com/css?family=Lato' \n";
+ Printf.printf "rel='stylesheet' type='text/css' />\n";
+ Printf.printf "<style>\n";
+ Printf.printf "body { font-family: 'Lato'; font-size: 14px; }\n";
+ Printf.printf "div#main { margin: 25px auto; width: 450px; }\n";
+ Printf.printf "h1 { font-size: 18px; }\n";
+ Printf.printf "</style>\n";
+ Printf.printf "\n";
+ Printf.printf "</head>\n";
+ Printf.printf "\n";
+ Printf.printf "<body>\n";
+ Printf.printf "\n";
+ Printf.printf "<div id=\"main\">\n";
+ Printf.printf "<h1>Common Factors in Bareiss' Algorithm over GF(2)[<i>x</i>]</h1>\n";
+ Printf.printf "\n";
+ Printf.printf "\n";
+ Printf.printf "<svg width=\"450\" height=\"450\">\n";
+ Printf.printf "\n";
+ Printf.printf "<defs>\n";
+ Printf.printf "<marker id=\"arrow-end\"\n";
+ Printf.printf "markerWidth=\"7\" markerHeight=\"7\"\n";
+ Printf.printf "refX=\"2\" refY=\"4\"\n";
+ Printf.printf "orient=\"auto\">\n";
+ Printf.printf "<path d=\"M2,2 L2,6 L5,4 L2,2\" style=\"fill: #000000;\" />\n";
+ Printf.printf "</marker>\n";
+ Printf.printf "\n";
+ Printf.printf "<marker id=\"arrow-begin\"\n";
+ Printf.printf "markerWidth=\"7\" markerHeight=\"7\"\n";
+ Printf.printf "refX=\"5\" refY=\"4\"\n";
+ Printf.printf "orient=\"auto\">\n";
+ Printf.printf "<path d=\"M5,2 L5,6 L2,4 L5,2\" style=\"fill: #000000;\" />\n";
+ Printf.printf "</marker>\n";
+ Printf.printf "\n";
+ Printf.printf "\n";
+ Printf.printf "<marker id=\"circle-total\"\n";
+ Printf.printf "markerWidth=\"4\" markerHeight=\"4\"\n";
+ Printf.printf "refx=\"2\" refy=\"2\">\n";
+ Printf.printf "<circle cx=\"2\" cy=\"2\" r=\"2\" stroke=\"none\" ";
+ Printf.printf "fill=\"%s\"/>\n" total;
+ Printf.printf "</marker>\n";
+ Printf.printf "\n";
+ Printf.printf "<marker id=\"circle-prediction\"\n";
+ Printf.printf "markerWidth=\"4\" markerHeight=\"4\"\n";
+ Printf.printf "refx=\"2\" refy=\"2\">\n";
+ Printf.printf "<circle cx=\"2\" cy=\"2\" r=\"2\" stroke=\"none\" ";
+ Printf.printf "fill=\"%s\"/>\n" predict;
+ Printf.printf "</marker>\n";
+ Printf.printf "\n";
+ Printf.printf "</defs>\n";
+ Printf.printf "\n"
+
+let total = "#6720A1"
+and predict = "#B31E61"
+
+let total2 = "#553B75"
+and predict2 = "#BD5500"
+
+let foot () =
+ Printf.printf "</svg>\n";
+ Printf.printf "<p>\n";
+ Printf.printf "Number of common row factors of <i>U</i> in the\n";
+ Printf.printf "<i>LD<sup>-1</sup>U</i> factorisation of <i>A</i>. Input were\n";
+ Printf.printf "300 random <i>n</i>-by-<i>n</i> matrices <i>A</i> over\n";
+ Printf.printf "GF(2)[<i>x</i>] with entries of degrees at most 10. Row\n";
+ Printf.printf "factors ignore the last row since it contains just\n";
+ Printf.printf "det <i>A</i>. On average, we can predict %.2f%% of the\n"
+ (100. *. mean_ratio);
+ Printf.printf "common factors (with a devation of %.2f%%).\n"
+ (100. *. deviation_ratio);
+ Printf.printf "</p>\n";
+ Printf.printf "</div>\n";
+ Printf.printf "\n";
+ Printf.printf "\n";
+ Printf.printf "</body>\n";
+ Printf.printf "</html>\n";
+ Printf.printf "\n"
+
+
+let main =
+ head total predict;
+ axes ();
+ mean
+ |> List.map (fun (n,(_,t)) -> (float n, t))
+ |> svg_poly_line total "total";
+ mean
+ |> List.map (fun (n,(p,_)) -> (float n, p))
+ |> svg_poly_line predict "prediction";
+ foot ()
--- /dev/null
+(* Implementation of reduction of a matrix into upper triangular form
+ using Gaussian elimination, cross-multiplication, or Bareiss
+ one-step algorithm. This program does only compute the resulting
+ matrix but not the transformation matrices. *)
+
+
+(* Define a ring with exact division. *)
+
+module type RING = sig
+ type t
+
+ val eql : t -> t -> bool
+
+ val add : t -> t -> t
+ val sub : t -> t -> t
+ val neg : t -> t
+ val nul : t
+
+ val mul : t -> t -> t
+ val div : t -> t -> t
+ val inv : t -> t
+ val one : t
+
+ val exact_div : t -> t -> t
+
+ val show : t -> string
+
+ val of_int : int -> t
+
+ val size : t -> int (* Used for pivoting *)
+ val gcd : t -> t -> t (* Used for row GCDs *)
+ end
+
+
+(* The reductions are in this functor. *)
+
+module Matrix (R : RING) = struct
+ (* We implement matrices in an imperative (ie, destructive) style
+ using two-dimensional arrays. *)
+
+ type t = R.t array array
+
+ (* Pretty print a matrix in LaTeX format (assuming that the
+ "amsmath" package is used). *)
+
+ let as_latex a =
+ "\\begin{pmatrix}\n"
+ ^ (ExtArray.map_matrix R.show a
+ |> ExtArray.join_matrix ~colsep:" & " ~rowsep:" \\\\\n" ~align:`Left)
+ ^ "\n\\end{pmatrix}"
+ let as_maple a =
+ "<"
+ ^ (ExtArray.map_matrix R.show a
+ |> Array.map (ExtArray.join ~separator:", ")
+ |> ExtArray.join ~separator:"; ")
+ ^ ">"
+ let as_ascii a =
+ "[ "^
+ (ExtArray.map_matrix R.show a
+ |> ExtArray.join_matrix ~colsep:" " ~rowsep:" ]\n[ " ~align:`Right)
+ ^" ]"
+
+ (* Next are elementary row and column operations as needed for the
+ reductions. The matrix input is always called [a], and the
+ parameters [p] and [q] give a range of row (or columns) on which
+ the transformation acts (inclding [p] and excluding [q]). This
+ allows us, eg, to omit a number of pointless additions as the
+ lower rows become filled with zeroes during the reductions.
+
+ The function [addrow] adds [s] times the [j]-th row to the [i]-th
+ row. With [mulrow] we multiply the [i]-th row by [s]. Finally,
+ [swaprow] and [swapcol] exchange the [i]-th row (or column) and
+ the [j]-th row (or column). *)
+
+ let addrow a i j s p q =
+ let ai = a.(i) and aj = a.(j) in
+ for x = p to q-1 do
+ ai.(x) <- R.add ai.(x) (R.mul s aj.(x))
+ done
+ and multrow a i s p q =
+ let ai = a.(i) in
+ for x = p to q-1 do
+ ai.(x) <- R.mul s ai.(x)
+ done
+ and swaprow a i j (_p : int) (_q : int) =
+ let tmp = a.(i) in
+ a.(i) <- a.(j);
+ a.(j) <- tmp
+ and swapcol a i j p q =
+ for x = p to q-1 do
+ let tmp = a.(x).(i) in
+ a.(x).(i) <- a.(x).(j);
+ a.(x).(j) <- tmp
+ done
+
+ (* Currently, we always use the first possible pivot. In order to
+ abort the search as soon as a pivot is found, we use an exception
+ to get out of the loop in the function [find_pivot]
+ prematurely. Given a matrix [a] with [m] rows and [n] columns as
+ well as an index [k], this function will find a pivot in the
+ submatrix of [a] containing the rows [k] to [m] and the columns
+ [k] to [n]. It returns the pivot with its position. The function
+ [pivot] does the same search but already applies row and column
+ transformations which bring the pivot to position [(k,k)]. It
+ returns just the pivot. Both functions return [None] if no pivot
+ was found. *)
+
+
+ let find_pivot a k m n =
+ let rec loop candidate i j =
+ if j >= n
+ then Option.map (fun (p,i,j,_) -> (p,i,j)) candidate
+ else
+ if i >= m
+ then loop candidate k (j+1)
+ else
+ let c = a.(i).(j) in
+ if R.eql c R.nul
+ then loop candidate (i+1) j
+ else
+ let c_size = R.size c in
+ match candidate with
+ | Some (p,_,_,p_size) ->
+ if p_size > c_size
+ then loop (Some (c,i,j,c_size)) (i+1) j
+ else loop candidate (i+1) j
+ | _ -> loop (Some (c,i,j,c_size)) (i+1) j
+ in
+ loop None k k
+
+
+ let pivot p l a q k m n = match find_pivot a k m n with
+ | None -> None
+ | Some (pivot,i,j) ->
+ if i <> k then begin
+ swapcol p i k 0 m;
+ swaprow l i k 0 k;
+ swaprow a i k k m
+ end;
+ if j <> k then begin
+ swapcol a j k 0 n;
+ swaprow q j k 0 n
+ end;
+ Some pivot
+
+
+ (* Create a deep copy of a matrix, ie, a matrix containing the same
+ entries but not sharing any (pointers to) rows. *)
+
+ let deepcopy a = ExtArray.map_matrix (fun x -> x) a
+
+
+
+
+ (* The first method we implement is cross-multiplication. An
+ auxillary function [crossmult] does the elimination step. Its
+ inputs are a matrix [a] with [m] rows and [n] columns, the
+ current row/column index [k], the pivot [p]. The function
+ [crossmultiplication] has as input a matrix [a] and will return
+ the upper triangular reduction of [a]. *)
+
+ let crossmult l a k p m n =
+ l.(k).(k) <- p;
+ for i = k+1 to m-1 do
+ l.(i).(k) <- a.(i).(k);
+ multrow a i p k m
+ done;
+ for i = k+1 to m-1 do
+ addrow a i k (R.neg l.(i).(k)) k m
+ done
+
+ (* Bareiss algorithm extends cross-multiplication by an exact
+ division step. This step is carried out by the function [divide]
+ which needs the [m] by [n] matrix [a], the last pivot [oldp] and
+ [k] where [(k,k)] is the position of the current pivot. The
+ function [bareiss] applies the algorithm to a matrix [a]. *)
+
+ let divide a k oldp m n =
+ for i = k+1 to m-1 do
+ for j = k+1 to n-1 do
+ a.(i).(j) <- R.exact_div a.(i).(j) oldp
+ done
+ done
+
+
+ (* Create an [n]-by-[n] identity matrix or an [m]-by-[n] zero
+ matrix. *)
+
+ let identity n =
+ ExtArray.init_matrix n n (fun i j -> if i = j then R.one else R.nul)
+ and zero m n =
+ ExtArray.init_matrix m n (fun _ _ -> R.nul)
+
+
+ (* The main function for Bareiss' algorithm. *)
+
+ let bareiss a =
+ let a = deepcopy a in
+ let (m,n) = ExtArray.length_matrix a in
+ let p = identity m
+ and l = zero m m
+ and d = identity m
+ and q = identity n
+ in
+ let rec loop oldp k =
+ if k < m && k < n then
+ match pivot p l a q k m n with
+ | None -> (p,l,d,a,q)
+ | Some pivot ->
+ crossmult l a k pivot m n;
+ divide a k oldp m n;
+ d.(k).(k) <- R.mul oldp pivot;
+ loop pivot (k+1)
+ else (p,l,d,a,q)
+ in
+ loop R.one 0
+
+
+ (* Compute the row GCDs of a matrix. Result is an array. *)
+
+ let rowgcds = Array.map (Array.fold_left R.gcd R.nul)
+
+
+ (* Predict common factors in [U] using the new prediction method and
+ reading the values from [L]. We cannot predict anything for the
+ first row with this method and we leave out the (trivial) factor
+ in the last row. *)
+
+ let prediction l =
+ let (m,_n) = ExtArray.length_matrix l in
+ Array.init m (fun k ->
+ if k = 0 || k = m-1
+ then R.one
+ else
+ let g = R.gcd l.(k-1).(k-1) l.(k).(k-1)
+ in
+ if k = 1
+ then g
+ else R.exact_div g (R.gcd g l.(k-2).(k-2)))
+end
+
+
+
+
+
+
+
+
--- /dev/null
+(* This module contains an implementation of polynomials over GF(2)
+ using integers (ie, bit-vectors) to represent the polynomials. Note
+ that this restricts the degrees to at most 63. *)
+
+(* As mentioned above, we use integers to represent polynomials. The
+ [n]-th bit will correspond to the [n]-th coefficient. We define
+ some constants, too. *)
+
+type t = Polynom of Z.t
+
+let nul = Polynom (Z.of_int 0b0)
+and one = Polynom (Z.of_int 0b1)
+and x = Polynom (Z.of_int 0b10)
+
+(* Comparing polynomials is simply comparison of integers. *)
+
+let eql (Polynom a) (Polynom b) = Z.equal a b
+
+(* This function converts polynomials to strings using Maple
+ notation. *)
+
+let show (Polynom a) =
+ let rec loop accu a n =
+ if Z.equal a Z.zero
+ then
+ if accu = []
+ then "0"
+ else String.concat " + " accu
+ else
+ let accu' = match (Z.to_int (Z.logand a Z.one), n) with
+ | (1,0) -> "1" :: accu
+ | (1,1) -> "x" :: accu
+ | (1,n) -> ("x^" ^ string_of_int n) :: accu
+ | _ -> accu
+ in
+ loop accu' (Z.shift_right a 1) (n+1)
+ in
+ loop [] a 0
+
+
+(* Use the conversion above, we can define a printer for polynomials
+ at the toplevel. (Install with [#install_printer fprintf;;].) *)
+
+let fprintf ff a = Format.fprintf ff "%s" (show a)
+
+(* Addition (and subtraction) of polynomials is simply addition of
+ bit-vectors. *)
+
+let add (Polynom a) (Polynom b) = Polynom (Z.logxor a b)
+let sub = add and neg a = a
+
+(* Multiplication uses the naive algorithm. Based on this, we define
+ exponentiation using a library function. *)
+
+let mul (Polynom a) (Polynom b) =
+ let rec loop c b =
+ if Z.equal b Z.zero
+ then nul
+ else
+ let q = loop (Z.shift_left c 1) (Z.shift_right b 1) in
+ if Z.to_int (Z.logand b Z.one) = 1
+ then add (Polynom c) q
+ else q
+ in
+ loop a b
+
+let expt = Numbers.expt ~one ~times:mul
+
+
+(* The degree is simply the number of digits. *)
+
+let deg (Polynom a) =
+ let rec loop a =
+ if Z.equal a Z.zero
+ then -1
+ else 1 + loop (Z.shift_right a 1)
+ in
+ loop a
+
+
+(* The quotient and remainder after division of [a] by [b]. This uses
+ the naive algorithm. *)
+
+
+let quorem (Polynom a) (Polynom b) =
+ if Z.equal b Z.zero
+ then failwith "Division by zero."
+ else
+ let n = deg (Polynom b) in
+ let rec loop c =
+ let m = deg (Polynom c) - n in
+ if m < 0
+ then (Z.zero, c)
+ else
+ let (q,r) = loop (Z.logxor c (Z.shift_left b m))
+ in
+ (Z.logxor (Z.shift_left Z.one m) q, r)
+ in
+ let (q,r) = loop a in (Polynom q, Polynom r)
+
+let quo x y = fst (quorem x y)
+and rem x y = snd (quorem x y)
+
+
+(* Create random polynomials of at most a certain degree. A second
+ function creates random polynomials until an non-zero specimen is
+ found. *)
+
+let random deg =
+ let l = 1 lsl deg in
+ Polynom (Z.of_int (Random.int l))
+
+let random_nonzero deg =
+ let rec loop () =
+ let r = random deg in
+ if eql r nul then loop () else r
+ in loop ()
+
+
+(* The following functions will be used during Bareiss' algorithm. *)
+
+let exact_div a b =
+ let (q,r) = quorem a b in
+ if eql r nul then q else failwith "Division is not exact."
+
+let inv a = if eql a one then one else failwith "Not invertible."
+
+let div a b = if eql b one then a else failwith "Division by non-unit."
+
+let of_int a = Polynom (Z.of_int (a land 1))
+
+let size = deg
+
+
+(* Greatest common divisor. *)
+
+let rec gcd a b = if eql b nul then a else gcd b (rem a b)
+
+
+
--- /dev/null
+#!/bin/bash
+
+echo 'let data = ['
+
+cat DATA.* \
+ | sed 's/[a-z =]//g' \
+ | cut -d, -f1,3,5 \
+ | sed 's/\.$/;/' \
+
+echo ']'
--- /dev/null
+#!/bin/sh
+
+for N in `seq 5 5 80`
+do
+ for I in `seq 1 10`
+ do
+ nohup sh -c "./test.native $N 30 | maple -q > DATA.$N.$I" &
+ done
+done
+
+
+
--- /dev/null
+open Matrix
+
+module R = Polys_over_GF2
+module M = Matrix(R)
+module N = Matrix(Integers)
+
+let () = Random.self_init ()
+
+
+let testGF2 n =
+ let a = ExtArray.init_matrix n n (fun _ _ -> R.random 5) in
+ let (p,l,d,u,q) = M.bareiss a in
+ let g = M.rowgcds u and w = M.prediction l in
+ (* Print factors and prediction: *)
+ Printf.printf "G := %s:\n" (M.as_maple [|g|]);
+ Printf.printf "W := %s:\n" (M.as_maple [|w|]);
+ (* Check prediction: *)
+ assert(ExtArray.map2 R.rem g w
+ |> Array.map (R.eql R.nul)
+ |> Array.fold_left ( && ) true);
+ (* Print the results: *)
+ Printf.printf "A := %s:\n" (M.as_maple a);
+ Printf.printf "Pr := %s:\n" (M.as_maple p);
+ Printf.printf "LL := %s:\n" (M.as_maple l);
+ Printf.printf "DD := %s:\n" (M.as_maple d);
+ Printf.printf "UU := %s:\n" (M.as_maple u);
+ Printf.printf "Pc := %s:\n" (M.as_maple q);
+ (* Check decomposition: *)
+ Printf.printf "Dinv := Inverse(DD) mod 2:\n";
+ Printf.printf "A - Pr.LL.Dinv.UU.Pc:\n";
+ Printf.printf "simplify(%%) mod 2:\n";
+ Printf.printf "if not Equal(%%, ZeroMatrix(%d,%d)) then\n" n n;
+ Printf.printf " printf(\"Faulty Decomposition!\\n\"):\n";
+ Printf.printf " lprint(A):\n";
+ Printf.printf " `quit`(1):\n";
+ Printf.printf "end if:\n";
+ (* Total amount of factors (leave out last row): *)
+ Printf.printf "map(f -> Factors(f) mod 2, G[1..-2]):\n";
+ Printf.printf "map(f -> map(e -> e[2], f[2]), %%):\n";
+ Printf.printf "map(f -> foldl(`+`, 0, op(f)), %%):\n";
+ Printf.printf "foldl(`+`, 0, op(convert(%%,list))):\n";
+ Printf.printf "total := total + %%:\n";
+ (* Predicted amount of factors: *)
+ Printf.printf "map(f -> Factors(f) mod 2, W[1..-2]):\n";
+ Printf.printf "map(f -> map(e -> e[2], f[2]), %%):\n";
+ Printf.printf "map(f -> foldl(`+`, 0, op(f)), %%):\n";
+ Printf.printf "foldl(`+`, 0, op(convert(%%,list))):\n";
+ Printf.printf "predicted := predicted + %%:\n"
+
+let factorsGF2 n times =
+ Printf.printf "restart:\n";
+ Printf.printf "with(LinearAlgebra):\n";
+ Printf.printf "total := 0:\n";
+ Printf.printf "predicted := 0:\n";
+ for i = 1 to times do
+ testGF2 n
+ done;
+ Printf.printf "printf(\"n = %d, " n;
+ Printf.printf "total = %%d, average total = %%.5f, ";
+ Printf.printf "prediction = %%d, average prediction = %%.5f.\\n\", ";
+ Printf.printf "total, total/%d, predicted, predicted/%d);\n" times times
+
+
+
+
+let testInts () =
+ let b =
+ ExtArray.init_matrix 15 15 (fun _ _ -> 5 - Random.int 11)
+ |> ExtArray.map_matrix Integers.of_int
+ in
+ let (p,l,d,u,q) = N.bareiss b in
+ let g = N.rowgcds u
+ and w = N.prediction l
+ in
+
+ Printf.printf "G := %s:\n%!" (N.as_maple [|g|]);
+ Printf.printf "W := %s:\n%!" (N.as_maple [|w|]);
+
+ assert(ExtArray.map2 Integers.rem g w
+ |> Array.map (Integers.eql Integers.nul)
+ |> Array.fold_left ( && ) true);
+
+ Printf.printf "A - Pr.LL.DD^(-1).UU.Pc;\n"
+
+
+
+
+(* We test the arithmetic in this file by generating random
+ polynomials of a certain degree and checking whether the ring
+ axioms hold. We also check division with remainder. *)
+
+let test ?(deg=10) n =
+ let (+.) = R.add
+ and l = R.one
+ and ( *.) = R.mul
+ and (==) = R.eql
+ and o = R.nul
+ in
+ for i = 0 to n do
+ let a = R.random deg
+ and b = R.random_nonzero deg
+ and c = R.random deg in
+ assert (a *. (b *. c) == (a *. b) *. c);
+ assert ((a *. b) == (b *. a));
+ assert (l *. a == a);
+ assert (let (q,r) = R.quorem a b in a == q *. b +. r);
+ assert (let g = R.gcd a c in
+ g == o || begin
+ let (qa,ra) = R.quorem a g
+ and (qc,rc) = R.quorem c g
+ in
+ ra == o && rc == o && qa *. g == a && qc *. g == c
+ end)
+ done;
+ Printf.printf "### %d tests passed.\n" n
+
+
+let main =
+ if not !Sys.interactive then begin
+ if Array.length Sys.argv < 2
+ then failwith "Not enough arguments!"
+ else
+ let n = Sys.argv.(1) |> int_of_string
+ and times = Sys.argv.(2) |> int_of_string
+ in
+ factorsGF2 n times
+ end
+
--- /dev/null
+#!/bin/bash
+
+bash ./prep.sh > data.ml
+ocaml draw.ml > drawing.tex
+pdflatex graph.tex
+pdflatex graph.tex
projects / jmiddeke / LU.git / commitdiff