Added code to compute the probability of non-trivial predictions; I also fixed the...

diff --git a/polys_over_GF2.ml b/polys_over_GF2.ml
index bbe6b41..b215fe4 100644 (file)
--- a/polys_over_GF2.ml
+++ b/polys_over_GF2.ml
@@ -104,11 +104,12 @@ 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. *)
+   found. At the moment, we can only generate polynomials up to degree
+   [62]. Larger degrees will upset the random number generator. *)
 
 let random deg =
-  let l = 1 lsl deg in
-  Polynom (Z.of_int (Random.int l))
+  let l = Int64.shift_left 1L deg in
+  Polynom (Z.of_int64 (Random.int64 l))
 
 let random_nonzero deg =
   let rec loop () =

diff --git a/test.ml b/test.ml
index 38ca6ce..2db4727 100644 (file)
--- a/test.ml
+++ b/test.ml
@@ -119,6 +119,30 @@ let test deg n =
   Printf.printf "### %d tests passed.\n" n
 
 
+(* Figure out how large the probability might be that [gcd a b]
+   divided by [gcd a b c] is non-trivial. Here [size] is the maximal
+   degree of [a], [b] and [c]; and [times] is the number of samples to
+   be drawn. *)
+
+let prob size times =
+  let rec loop total n =
+    if n <= 0
+    then float total /. float times
+    else
+      let a = R.random size
+      and b = R.random size
+      and c = R.random_nonzero size
+      in
+      let g = R.gcd a b
+      in
+      let w = R.gcd g c
+      in
+      loop (if R.eql (R.quo g w) R.one then total else total+1) (n-1)
+  in
+  Printf.printf "size = %d, times = %d, nontrivial GCD = %.2f%%.\n"
+    size times (100. *. loop 0 times)
+    
+    
 (* Parsing command line arguments (just for fun). *)
 
 
@@ -142,6 +166,10 @@ let args =
     "--test",
     Arg.Unit (set_experiment test),
     "Test polynomial arithmetic.";
+
+    "--probability",
+    Arg.Unit (set_experiment prob),
+    "How often is GCD(a,b)/GCD(a,b,p) non-trivial in GF(2)[x]?";
     
     "--size",
     Arg.Int set_size,