diff options
author | letouzey <letouzey@85f007b7-540e-0410-9357-904b9bb8a0f7> | 2012-07-30 13:45:36 +0000 |
---|---|---|
committer | letouzey <letouzey@85f007b7-540e-0410-9357-904b9bb8a0f7> | 2012-07-30 13:45:36 +0000 |
commit | c2976e6279e9aa3a661e026957b058b9a00148b3 (patch) | |
tree | b81b242ac2b4eaec8bde20dd8c53d3992103a4b6 | |
parent | 00f429861cc6bae9110d82aac2512fe5e55342bc (diff) |
Bigint : better ensure canonicity of arrays of int blocks
In the chosen representation of negative numbers, (-base) is
actually legal as first element of an int array. This situation
was wrongly handled by many functions in this module, leading to
possible non-canonical representation of a big number, and hence
faulty "equal" answers or buggy output of "to_string".
For instance, on a 32-bit machine :
open Big_int;;
to_string (sub (of_string "-10000") (of_string "-1000000"));;
(* was displaying 9810000 instead of 990000 *)
Also simplified "of_string", which now rely on "neg" for handling
negative numbers.
Btw, improved checks comparing w.r.t OCaml's Big_int instead of float.
git-svn-id: svn+ssh://scm.gforge.inria.fr/svn/coq/trunk@15665 85f007b7-540e-0410-9357-904b9bb8a0f7
-rw-r--r-- | lib/bigint.ml | 250 |
1 files changed, 160 insertions, 90 deletions
diff --git a/lib/bigint.ml b/lib/bigint.ml index 1985eceb8..4215b0ab9 100644 --- a/lib/bigint.ml +++ b/lib/bigint.ml @@ -10,26 +10,29 @@ (* Basic operations on (unbounded) integer numbers *) (***************************************************) -(* An integer is canonically represented as an array of k-digits blocs. +(* An integer is canonically represented as an array of k-digits blocs, + i.e. in base 10^k. 0 is represented by the empty array and -1 by the singleton [|-1|]. - The first bloc is in the range ]0;10^k[ for positive numbers. - The first bloc is in the range ]-10^k;-1[ for negative ones. - All other blocs are numbers in the range [0;10^k[. + The first bloc is in the range ]0;base[ for positive numbers. + The first bloc is in the range [-base;-1[ for numbers < -1. + All other blocs are numbers in the range [0;base[. - Negative numbers are represented using 2's complementation. For instance, - with 4-digits blocs, [-9655;6789] denotes -96543211 -*) + Negative numbers are represented using 2's complementation : + one unit is "borrowed" from the top block for complementing + the other blocs. For instance, with 4-digits blocs, + [|-5;6789|] denotes -43211 + since -5.10^4+6789=-((4.10^4)+(10000-6789)) = -43211 -(* The base is a power of 10 in order to facilitate the parsing and printing + The base is a power of 10 in order to facilitate the parsing and printing of numbers in digital notation. All functions, to the exception of to_string and of_string should work with an arbitrary base, even if not a power of 10. - In practice, we set k=4 so that no overflow in ocaml machine words - (i.e. the interval [-2^30;2^30-1]) occur when multiplying two - numbers less than (10^k) + In practice, we set k=4 on 32-bits machines, so that no overflow in ocaml + machine words (i.e. the interval [-2^30;2^30-1]) occur when multiplying two + numbers less than (10^k). On 64-bits machines, k=9. *) (* The main parameters *) @@ -41,6 +44,7 @@ let size = let format_size = (* How to parametrize a printf format *) if size = 4 then Printf.sprintf "%04d" + else if size = 9 then Printf.sprintf "%09d" else fun n -> let rec aux j l n = if j=size then l else aux (j+1) (string_of_int (n mod 10) :: l) (n/10) @@ -50,44 +54,81 @@ let format_size = let base = let rec exp10 = function 0 -> 1 | n -> 10 * exp10 (n-1) in exp10 size +(******************************************************************) +(* First, we represent all numbers by int arrays. + Later, we will optimize the particular case of small integers *) +(******************************************************************) + +module ArrayInt = struct + (* Basic numbers *) let zero = [||] let neg_one = [|-1|] -(* Sign of an integer *) -let is_strictly_neg n = n<>[||] && n.(0) < 0 -let is_strictly_pos n = n<>[||] && n.(0) > 0 -let is_neg_or_zero n = n=[||] or n.(0) < 0 -let is_pos_or_zero n = n=[||] or n.(0) > 0 +(* An array is canonical when + - it is empty + - it is [|-1|] + - its first bloc is in [-base;-1[U]0;base[ + and the other blocs are in [0;base[. *) + +let canonical n = + let ok x = (0 <= x && x < base) in + let rec ok_tail k = (k = 0) || (ok n.(k) && ok_tail (k-1)) in + let ok_init x = (-base <= x && x < base && x <> -1 && x <> 0) + in + (n = [||]) || (n = [|-1|]) || + (ok_init n.(0) && ok_tail (Array.length n - 1)) + +(* [normalize_pos] : removing initial blocks of 0 *) let normalize_pos n = let k = ref 0 in while !k < Array.length n & n.(!k) = 0 do incr k done; Array.sub n !k (Array.length n - !k) +(* [normalize_neg] : avoid (-1) as first bloc. + input: an array with -1 as first bloc and other blocs in [0;base[ + output: a canonical array *) + let normalize_neg n = let k = ref 1 in while !k < Array.length n & n.(!k) = base - 1 do incr k done; let n' = Array.sub n !k (Array.length n - !k) in if Array.length n' = 0 then [|-1|] else (n'.(0) <- n'.(0) - base; n') +(* [normalize] : avoid 0 and (-1) as first bloc. + input: an array with first bloc in [-base;base[ and others in [0;base[ + output: a canonical array *) + let rec normalize n = - if Array.length n = 0 then n else - if n.(0) = -1 then normalize_neg n else normalize_pos n + if Array.length n = 0 then n + else if n.(0) = -1 then normalize_neg n + else if n.(0) = 0 then normalize_pos n + else n + +(* Opposite (expects and returns canonical arrays) *) let neg m = if m = zero then zero else let n = Array.copy m in let i = ref (Array.length m - 1) in while !i > 0 & n.(!i) = 0 do decr i done; - if !i > 0 then begin + if !i = 0 then begin + n.(0) <- - n.(0); + (* n.(0) cannot be 0 since m is canonical *) + if n.(0) = -1 then normalize_neg n + else if n.(0) = base then (n.(0) <- 0; Array.append [| 1 |] n) + else n + end else begin + (* here n.(!i) <> 0, hence 0 < base - n.(!i) < base for n canonical *) n.(!i) <- base - n.(!i); decr i; while !i > 0 do n.(!i) <- base - 1 - n.(!i); decr i done; + (* since -base <= n.(0) <= base-1, hence -base <= -n.(0)-1 <= base-1 *) n.(0) <- - n.(0) - 1; - if n.(0) < -1 then (n.(0) <- n.(0) + base; Array.append [| -1 |] n) else - if n.(0) = - base then (n.(0) <- 0; Array.append [| -1 |] n) - else normalize n - end else (n.(0) <- - n.(0); n) + (* since m is canonical, m.(0)<>0 hence n.(0)<>-1, + and m=-1 is already handled above, so here m.(0)<>-1 hence n.(0)<>0 *) + n + end let push_carry r j = let j = ref j in @@ -97,10 +138,10 @@ let push_carry r j = while !j > 0 & r.(!j) >= base do r.(!j) <- r.(!j) - base; decr j; r.(!j) <- r.(!j) + 1 done; + (* here r.(0) could be in [-2*base;2*base-1] *) if r.(0) >= base then (r.(0) <- r.(0) - base; Array.append [| 1 |] r) else if r.(0) < -base then (r.(0) <- r.(0) + 2*base; Array.append [| -2 |] r) - else if r.(0) = -base then (r.(0) <- 0; Array.append [| -1 |] r) - else normalize r + else normalize r (* in case r.(0) is 0 or -1 *) let add_to r a j = if a = zero then r else begin @@ -148,6 +189,13 @@ let rec mult m n = done; normalize r +(* Comparisons *) + +let is_strictly_neg n = n<>[||] && n.(0) < 0 +let is_strictly_pos n = n<>[||] && n.(0) > 0 +let is_neg_or_zero n = n=[||] or n.(0) < 0 +let is_pos_or_zero n = n=[||] or n.(0) > 0 + let rec less_than_same_size m n i j = i < Array.length m && (m.(i) < n.(j) or (m.(i) = n.(j) && less_than_same_size m n (i+1) (j+1))) @@ -160,6 +208,8 @@ let less_than m n = is_strictly_pos n && (Array.length m < Array.length n or (Array.length m = Array.length n && less_than_same_size m n 0 0)) +(* For this equality test it is critical that n and m are canonical *) + let equal m n = (m = n) let less_than_shift_pos k m n = @@ -187,6 +237,11 @@ let sub_mult m d q k = end done +(** Euclid division m/d = (q,r) + This is the "Floor" variant, as with ocaml's / + (but not as ocaml's Big_int.quomod_big_int). + We have sign r = sign m *) + let euclid m d = let isnegm, m = if is_strictly_neg m then (-1),neg m else 1,Array.copy m in @@ -224,33 +279,21 @@ let euclid m d = (* Parsing/printing ordinary 10-based numbers *) let of_string s = - let isneg = String.length s > 1 & s.[0] = '-' in - let n = if isneg then 1 else 0 in - let d = ref n in - while !d < String.length s && s.[!d] = '0' do incr d done; - if !d = String.length s then zero else - let r = (String.length s - !d) mod size in + let len = String.length s in + let isneg = len > 1 & s.[0] = '-' in + let d = ref (if isneg then 1 else 0) in + while !d < len && s.[!d] = '0' do incr d done; + if !d = len then zero else + let r = (len - !d) mod size in let h = String.sub s (!d) r in - if !d = String.length s - 1 && isneg && h="1" then neg_one else let e = if h<>"" then 1 else 0 in - let l = (String.length s - !d) / size in - let a = Array.create (l + e + n) 0 in - if isneg then begin - a.(0) <- (-1); - let carry = ref 0 in - for i=l downto 1 do - let v = int_of_string (String.sub s ((i-1)*size + !d +r) size)+ !carry in - if v <> 0 then (a.(i+e)<- base - v; carry := 1) else carry := 0 - done; - if e=1 then a.(1) <- base - !carry - int_of_string h; - end - else begin - if e=1 then a.(0) <- int_of_string h; - for i=1 to l do - a.(i+e-1) <- int_of_string (String.sub s ((i-1)*size + !d + r) size) - done - end; - a + let l = (len - !d) / size in + let a = Array.create (l + e) 0 in + if e=1 then a.(0) <- int_of_string h; + for i=1 to l do + a.(i+e-1) <- int_of_string (String.sub s ((i-1)*size + !d + r) size) + done; + if isneg then neg a else a let to_string_pos sgn n = if Array.length n = 0 then "0" else @@ -262,25 +305,37 @@ let to_string n = if is_strictly_neg n then to_string_pos "-" (neg n) else to_string_pos "" n +end + (******************************************************************) (* Optimized operations on (unbounded) integer numbers *) (* integers smaller than base are represented as machine integers *) (******************************************************************) +open ArrayInt + type bigint = Obj.t +(* NB: n should be here in ]-base*base ; base*base] *) + +let array_of_int n = + let lo = n mod base + and hi = n / base in + if hi = base then [|1;0;lo|] + else if hi < 0 && lo <> 0 then [|hi-1;lo+base|] + else [|hi;lo|] + let ints_of_int n = - if n >= base then [| n / base; n mod base |] - else if n <= - base then [| n / base - 1; n mod base + base |] - else if n = 0 then [| |] else [| n |] + if n = 0 then [| |] + else if -base <= n && n < base then [| n |] + else array_of_int n let big_of_int n = - if n >= base then Obj.repr [| n / base; n mod base |] - else if n <= - base then Obj.repr [| n / base - 1; n mod base + base |] - else Obj.repr n + if -base <= n && n < base then Obj.repr n + else Obj.repr (array_of_int n) let big_of_ints n = - let n = normalize n in + let n = normalize n in (* TODO: using normalize here seems redundant now *) if n = zero then Obj.repr 0 else if Array.length n = 1 then Obj.repr n.(0) else Obj.repr n @@ -344,6 +399,8 @@ let is_strictly_pos n = is_strictly_pos (ints_of_z n) let is_neg_or_zero n = is_neg_or_zero (ints_of_z n) let is_pos_or_zero n = is_pos_or_zero (ints_of_z n) +let equal m n = (m = n) + (* spiwack: computes n^m *) (* The basic idea of the algorithm is that n^(2m) = (n^2)^m *) (* In practice the algorithm performs : @@ -367,43 +424,56 @@ let pow = in pow_aux one -(* Testing suite *) +(** Testing suite w.r.t. OCaml's Big_int *) + +(* +module B = struct + open Big_int + let zero = zero_big_int + let to_string = string_of_big_int + let of_string = big_int_of_string + let add = add_big_int + let opp = minus_big_int + let sub = sub_big_int + let mul = mult_big_int + let abs = abs_big_int + let sign = sign_big_int + let euclid n m = + let n' = abs n and m' = abs m in + let q',r' = quomod_big_int n' m' in + (if sign (mul n m) < 0 && sign q' <> 0 then opp q' else q'), + (if sign n < 0 then opp r' else r') +end let check () = - let numbers = [ - "1";"2";"99";"100";"101";"9999";"10000";"10001"; - "999999";"1000000";"1000001";"99999999";"100000000";"100000001"; - "1234";"5678";"12345678";"987654321"; - "-1";"-2";"-99";"-100";"-101";"-9999";"-10000";"-10001"; - "-999999";"-1000000";"-1000001";"-99999999";"-100000000";"-100000001"; - "-1234";"-5678";"-12345678";"-987654321";"0" - ] + let roots = [ 1; 100; base; 100*base; base*base ] in + let rands = [ 1234; 5678; 12345678; 987654321 ] in + let nums = (List.flatten (List.map (fun x -> [x-1;x;x+1]) roots)) @ rands in + let numbers = + List.map string_of_int nums @ + List.map (fun n -> string_of_int (-n)) nums in - let eucl n m = - let n' = abs_float n and m' = abs_float m in - let q' = floor (n' /. m') in let r' = n' -. m' *. q' in - (if n *. m < 0. & q' <> 0. then -. q' else q'), - (if n < 0. then -. r' else r') in - let round f = floor (abs_float f +. 0.5) *. (if f < 0. then -1. else 1.) in let i = ref 0 in - let compare op n n' = + let compare op x y n n' = incr i; let s = Printf.sprintf "%30s" (to_string n) in - let s' = Printf.sprintf "% 30.0f" (round n') in - if s <> s' then Printf.printf "%s: %s <> %s\n" op s s' in -List.iter (fun a -> List.iter (fun b -> - let n = of_string a and m = of_string b in - let n' = float_of_string a and m' = float_of_string b in - let a = add n m and a' = n' +. m' in - let s = sub n m and s' = n' -. m' in - let p = mult n m and p' = n' *. m' in - let q,r = try euclid n m with Division_by_zero -> zero,zero - and q',r' = eucl n' m' in - compare "+" a a'; - compare "-" s s'; - compare "*" p p'; - compare "/" q q'; - compare "%" r r') numbers) numbers; + let s' = Printf.sprintf "%30s" (B.to_string n') in + if s <> s' then Printf.printf "%s%s%s: %s <> %s\n" x op y s s' in + let test x y = + let n = of_string x and m = of_string y in + let n' = B.of_string x and m' = B.of_string y in + let a = add n m and a' = B.add n' m' in + let s = sub n m and s' = B.sub n' m' in + let p = mult n m and p' = B.mul n' m' in + let q,r = try euclid n m with Division_by_zero -> zero,zero + and q',r' = try B.euclid n' m' with Division_by_zero -> B.zero, B.zero + in + compare "+" x y a a'; + compare "-" x y s s'; + compare "*" x y p p'; + compare "/" x y q q'; + compare "%" x y r r' + in + List.iter (fun a -> List.iter (test a) numbers) numbers; Printf.printf "%i tests done\n" !i - - +*) |