From b84af069df5c28f54a0343e099a8b989c9ded616 Mon Sep 17 00:00:00 2001 From: Johannes Middeke Date: Wed, 4 May 2016 05:59:35 +0200 Subject: [PATCH 1/1] Initial commit. The current state of the program computes the number of common factors. --- _oasis | 14 +++ data.ml | 72 ++++++++++++++++ draw.ml | 69 +++++++++++++++ html.sh | 9 ++ integers.ml | 32 +++++++ makeHTML.ml | 220 ++++++++++++++++++++++++++++++++++++++++++++++++ matrix.ml | 248 ++++++++++++++++++++++++++++++++++++++++++++++++++++++ polys_over_GF2.ml | 140 ++++++++++++++++++++++++++++++ prep.sh | 10 +++ run.sh | 12 +++ test.ml | 128 ++++++++++++++++++++++++++++ visualise.sh | 6 ++ 12 files changed, 960 insertions(+) create mode 100644 _oasis create mode 100644 data.ml create mode 100644 draw.ml create mode 100644 html.sh create mode 100644 integers.ml create mode 100644 makeHTML.ml create mode 100644 matrix.ml create mode 100644 polys_over_GF2.ml create mode 100644 prep.sh create mode 100644 run.sh create mode 100644 test.ml create mode 100644 visualise.sh diff --git a/_oasis b/_oasis new file mode 100644 index 0000000..d473356 --- /dev/null +++ b/_oasis @@ -0,0 +1,14 @@ +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 diff --git a/data.ml b/data.ml new file mode 100644 index 0000000..f61ea38 --- /dev/null +++ b/data.ml @@ -0,0 +1,72 @@ +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; +] diff --git a/draw.ml b/draw.ml new file mode 100644 index 0000000..e0ea276 --- /dev/null +++ b/draw.ml @@ -0,0 +1,69 @@ +#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 diff --git a/html.sh b/html.sh new file mode 100644 index 0000000..983128e --- /dev/null +++ b/html.sh @@ -0,0 +1,9 @@ +#!/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 diff --git a/integers.ml b/integers.ml new file mode 100644 index 0000000..283215e --- /dev/null +++ b/integers.ml @@ -0,0 +1,32 @@ +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 + diff --git a/makeHTML.ml b/makeHTML.ml new file mode 100644 index 0000000..5465a50 --- /dev/null +++ b/makeHTML.ml @@ -0,0 +1,220 @@ +#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 " 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 " %s " 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 "\n"; + Printf.printf "\n"; + Printf.printf " %d\n" x; + Printf.printf "\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 "\n"; + Printf.printf "\n"; + Printf.printf " %d\n" y; + Printf.printf "\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 "\n"; + Printf.printf "\n"; + Printf.printf "\n"; + Printf.printf " n\n"; + Printf.printf "\n"; + Printf.printf "\n"; + Printf.printf " # factors\n"; + Printf.printf "\n"; + Printf.printf "\n"; + x_ticks ~by:5 5 75; + y_ticks ~by:5 5 75 + +let head total predict = + Printf.printf "\n"; + Printf.printf "\n"; + Printf.printf "\n"; + Printf.printf "Bareiss & Polynomials over GF(2)[x]\n"; + Printf.printf "\n"; + Printf.printf "\n"; + Printf.printf "\n"; + Printf.printf "\n"; + Printf.printf "\n"; + Printf.printf "\n"; + Printf.printf "\n"; + Printf.printf "\n"; + Printf.printf "
\n"; + Printf.printf "

Common Factors in Bareiss' Algorithm over GF(2)[x]

\n"; + Printf.printf "\n"; + Printf.printf "\n"; + Printf.printf "\n"; + Printf.printf "\n"; + Printf.printf "\n"; + Printf.printf "\n"; + Printf.printf "\n"; + Printf.printf "\n"; + Printf.printf "\n"; + Printf.printf "\n"; + Printf.printf "\n"; + Printf.printf "\n"; + Printf.printf "\n"; + Printf.printf "\n"; + Printf.printf "\n"; + Printf.printf "\n" total; + Printf.printf "\n"; + Printf.printf "\n"; + Printf.printf "\n"; + Printf.printf "\n" predict; + Printf.printf "\n"; + Printf.printf "\n"; + Printf.printf "\n"; + Printf.printf "\n" + +let total = "#6720A1" +and predict = "#B31E61" + +let total2 = "#553B75" +and predict2 = "#BD5500" + +let foot () = + Printf.printf "\n"; + Printf.printf "

\n"; + Printf.printf "Number of common row factors of U in the\n"; + Printf.printf "LD-1U factorisation of A. Input were\n"; + Printf.printf "300 random n-by-n matrices A over\n"; + Printf.printf "GF(2)[x] with entries of degrees at most 10. Row\n"; + Printf.printf "factors ignore the last row since it contains just\n"; + Printf.printf "det A. 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 "

\n"; + Printf.printf "
\n"; + Printf.printf "\n"; + Printf.printf "\n"; + Printf.printf "\n"; + Printf.printf "\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 () diff --git a/matrix.ml b/matrix.ml new file mode 100644 index 0000000..2fca419 --- /dev/null +++ b/matrix.ml @@ -0,0 +1,248 @@ +(* 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 + + + + + + + + diff --git a/polys_over_GF2.ml b/polys_over_GF2.ml new file mode 100644 index 0000000..bbe6b41 --- /dev/null +++ b/polys_over_GF2.ml @@ -0,0 +1,140 @@ +(* 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) + + + diff --git a/prep.sh b/prep.sh new file mode 100644 index 0000000..6331e71 --- /dev/null +++ b/prep.sh @@ -0,0 +1,10 @@ +#!/bin/bash + +echo 'let data = [' + +cat DATA.* \ + | sed 's/[a-z =]//g' \ + | cut -d, -f1,3,5 \ + | sed 's/\.$/;/' \ + +echo ']' diff --git a/run.sh b/run.sh new file mode 100644 index 0000000..cd78b61 --- /dev/null +++ b/run.sh @@ -0,0 +1,12 @@ +#!/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 + + + diff --git a/test.ml b/test.ml new file mode 100644 index 0000000..631ae20 --- /dev/null +++ b/test.ml @@ -0,0 +1,128 @@ +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 + diff --git a/visualise.sh b/visualise.sh new file mode 100644 index 0000000..16ac011 --- /dev/null +++ b/visualise.sh @@ -0,0 +1,6 @@ +#!/bin/bash + +bash ./prep.sh > data.ml +ocaml draw.ml > drawing.tex +pdflatex graph.tex +pdflatex graph.tex -- 2.1.4