/lib/bigint.ml

(************************************************************************)
(*  v      *   The Coq Proof Assistant  /  The Coq Development Team     *)
(* <O___,, *   INRIA - CNRS - LIX - LRI - PPS - Copyright 1999-2012     *)
(*   \VV/  **************************************************************)
(*    //   *      This file is distributed under the terms of the       *)
(*         *       GNU Lesser General Public License Version 2.1        *)
(************************************************************************)

(***************************************************)
(* Basic operations on (unbounded) integer numbers *)
(***************************************************)

(* 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;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 :
   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
   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 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 *)

let size =
  let rec log10 n = if n < 10 then 0 else 1 + log10 (n / 10) in
  (log10 max_int) / 2

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)
    in String.concat "" (aux 0 [] n)

(* The base is 10^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|]

(* 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 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
    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;
    (* 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
  while !j > 0 & r.(!j) < 0 do
    r.(!j) <- r.(!j) + base; decr j; r.(!j) <- r.(!j) - 1
  done;
  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 normalize r (* in case r.(0) is 0 or -1 *)

let add_to r a j =
  if a = zero then r else begin
  for i = Array.length r - 1 downto j+1 do
    r.(i) <- r.(i) + a.(i-j);
    if r.(i) >= base then (r.(i) <- r.(i) - base; r.(i-1) <- r.(i-1) + 1)
  done;
  r.(j) <- r.(j) + a.(0);
  push_carry r j
  end

let add n m =
  let d = Array.length n - Array.length m in
  if d > 0 then add_to (Array.copy n) m d else add_to (Array.copy m) n (-d)

let sub_to r a j =
  if a = zero then r else begin
  for i = Array.length r - 1 downto j+1 do
    r.(i) <- r.(i) - a.(i-j);
    if r.(i) < 0 then (r.(i) <- r.(i) + base; r.(i-1) <- r.(i-1) - 1)
  done;
  r.(j) <- r.(j) - a.(0);
  push_carry r j
  end

let sub n m =
  let d = Array.length n - Array.length m in
  if d >= 0 then sub_to (Array.copy n) m d
  else let r = neg m in add_to r n (Array.length r - Array.length n)

let rec mult m n =
  if m = zero or n = zero then zero else
  let l = Array.length m + Array.length n in
  let r = Array.create l 0 in
  for i = Array.length m - 1 downto 0 do
    for j = Array.length n - 1 downto 0 do
      let p = m.(i) * n.(j) + r.(i+j+1) in
      let (q,s) =
        if p < 0
        then (p + 1) / base - 1, (p + 1) mod base + base - 1
        else p / base, p mod base in
      r.(i+j+1) <- s;
      if q <> 0 then r.(i+j) <- r.(i+j) + q;
    done
  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)))

let less_than m n =
  if is_strictly_neg m then
    is_pos_or_zero n or Array.length m > Array.length n
      or (Array.length m = Array.length n && less_than_same_size m n 0 0)
  else
    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 =
  (Array.length m - k < Array.length n)
  or (Array.length m - k = Array.length n && less_than_same_size m n k 0)

let rec can_divide k m d i =
  (i = Array.length d) or
  (m.(k+i) > d.(i)) or
  (m.(k+i) = d.(i) && can_divide k m d (i+1))

(* For two big nums m and d and a small number q,
   computes m - d * q * base^(|m|-|d|-k) in-place (in m).
   Both m d and q are positive. *)

let sub_mult m d q k =
  if q <> 0 then
  for i = Array.length d - 1 downto 0 do
    let v = d.(i) * q in
    m.(k+i) <- m.(k+i) - v mod base;
    if m.(k+i) < 0 then (m.(k+i) <- m.(k+i) + base; m.(k+i-1) <- m.(k+i-1) -1);
    if v >= base then begin
      m.(k+i-1) <- m.(k+i-1) - v / base;
      let j = ref (i-1) in
      while m.(k + !j) < 0 do (* result is positive, hence !j remains >= 0 *)
        m.(k + !j) <- m.(k + !j) + base; decr j; m.(k + !j) <- m.(k + !j) -1
      done
    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
  let isnegd, d = if is_strictly_neg d then (-1),neg d else 1,d in
  if d = zero then raise Division_by_zero;
  let q,r =
    if less_than m d then (zero,m) else
    let ql = Array.length m - Array.length d in
    let q = Array.create (ql+1) 0 in
    let i = ref 0 in
    while not (less_than_shift_pos !i m d) do
      if m.(!i)=0 then incr i else
      if can_divide !i m d 0 then begin
        let v =
          if Array.length d > 1 && d.(0) <> m.(!i) then
            (m.(!i) * base + m.(!i+1)) / (d.(0) * base + d.(1) + 1)
          else
            m.(!i) / d.(0) in
        q.(!i) <- q.(!i) + v;
	sub_mult m d v !i
      end else begin
        let v = (m.(!i) * base + m.(!i+1)) / (d.(0) + 1) in
        q.(!i) <- q.(!i) + v / base;
	sub_mult m d (v / base) !i;
        q.(!i+1) <- q.(!i+1) + v mod base;
	if q.(!i+1) >= base then
	  (q.(!i+1) <- q.(!i+1)-base; q.(!i) <- q.(!i)+1);
	sub_mult m d (v mod base) (!i+1)
      end
    done;
    (normalize q, normalize m) in
  (if isnegd * isnegm = -1 then neg q else q),
  (if isnegm = -1 then neg r else r)

(* Parsing/printing ordinary 10-based numbers *)

let of_string s =
  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
  let e = if h<>"" then 1 else 0 in
  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
     sgn ^
     String.concat ""
      (string_of_int n.(0) :: List.map format_size (List.tl (Array.to_list n)))

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

(* Since base is the largest power of 10 such that base*base <= max_int,
   we have max_int < 100*base*base : any int can be represented
   by at most three blocs *)

let small n = (-base <= n) && (n < base)

let mkarray n =
  (* n isn't small, this case is handled separately below *)
  let lo = n mod base
  and hi = n / base in
  let t = if small hi then [|hi;lo|] else [|hi/base;hi mod base;lo|]
  in
  for i = Array.length t -1 downto 1 do
    if t.(i) < 0 then (t.(i) <- t.(i) + base; t.(i-1) <- t.(i-1) -1)
  done;
  t

let ints_of_int n =
  if n = 0 then [| |]
  else if small n then [| n |]
  else mkarray n

let of_int n =
  if small n then Obj.repr n else Obj.repr (mkarray n)

let of_ints n =
  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

let coerce_to_int = (Obj.magic : Obj.t -> int)
let coerce_to_ints = (Obj.magic : Obj.t -> int array)

let to_ints n =
  if Obj.is_int n then ints_of_int (coerce_to_int n)
  else coerce_to_ints n

let int_of_ints =
  let maxi = mkarray max_int and mini = mkarray min_int in
  fun t ->
    let l = Array.length t in
    if (l > 3) || (l = 3 && (less_than maxi t || less_than t mini))
    then failwith "Bigint.to_int: too large";
    let sum = ref 0 in
    let pow = ref 1 in
    for i = l-1 downto 0 do
      sum := !sum + t.(i) * !pow;
      pow := !pow*base;
    done;
    !sum

let to_int n =
  if Obj.is_int n then coerce_to_int n
  else int_of_ints (coerce_to_ints n)

let app_pair f (m, n) =
  (f m, f n)

let add m n =
  if Obj.is_int m & Obj.is_int n
  then of_int (coerce_to_int m + coerce_to_int n)
  else of_ints (add (to_ints m) (to_ints n))

let sub m n =
  if Obj.is_int m & Obj.is_int n
  then of_int (coerce_to_int m - coerce_to_int n)
  else of_ints (sub (to_ints m) (to_ints n))

let mult m n =
  if Obj.is_int m & Obj.is_int n
  then of_int (coerce_to_int m * coerce_to_int n)
  else of_ints (mult (to_ints m) (to_ints n))

let euclid m n =
  if Obj.is_int m & Obj.is_int n
  then app_pair of_int
    (coerce_to_int m / coerce_to_int n, coerce_to_int m mod coerce_to_int n)
  else app_pair of_ints (euclid (to_ints m) (to_ints n))

let less_than m n =
  if Obj.is_int m & Obj.is_int n
  then coerce_to_int m < coerce_to_int n
  else less_than (to_ints m) (to_ints n)

let neg n =
  if Obj.is_int n then of_int (- (coerce_to_int n))
  else of_ints (neg (to_ints n))

let of_string m = of_ints (of_string m)
let to_string m = to_string (to_ints m)

let zero = of_int 0
let one = of_int 1
let two = of_int 2
let sub_1 n = sub n one
let add_1 n = add n one
let mult_2 n = add n n

let div2_with_rest n =
  let (q,b) = euclid n two in
  (q, b = one)

let is_strictly_neg n = is_strictly_neg (to_ints n)
let is_strictly_pos n = is_strictly_pos (to_ints n)
let is_neg_or_zero n = is_neg_or_zero (to_ints n)
let is_pos_or_zero n = is_pos_or_zero (to_ints 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 :
    k*n^0 = k
    k*n^(2m) = k*(n*n)^m
    k*n^(2m+1) = (n*k)*(n*n)^m *)
let pow  =
  let rec pow_aux odd_rest n m = (* odd_rest is the k from above *)
    if m<=0 then
      odd_rest
    else
      let quo = m lsr 1 (* i.e. m/2 *)
      and odd = (m land 1) <> 0 in
      pow_aux
	(if odd then mult n odd_rest else odd_rest)
	(mult n n)
	quo
  in
  pow_aux one

(** 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 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 i = ref 0 in
  let compare op x y n n' =
    incr i;
    let s = Printf.sprintf "%30s" (to_string n) in
    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
*)