z.mlp

(**
   Integers.


   This file is part of the Zarith library 
   http://forge.ocamlcore.org/projects/zarith .
   It is distributed under LGPL 2 licensing, with static linking exception.
   See the LICENSE file included in the distribution.
   
   Copyright (c) 2010-2011 Antoine Miné, Abstraction project.
   Abstraction is part of the LIENS (Laboratoire d'Informatique de l'ENS),
   a joint laboratory by:
   CNRS (Centre national de la recherche scientifique, France),
   ENS (École normale supérieure, Paris, France),
   INRIA Rocquencourt (Institut national de recherche en informatique, France).

 *)

type t

exception Overflow

external init: unit -> unit = "ml_z_init"
let _ = init ()

external install_frametable: unit -> unit = install_frametable@ASM
let _ =  install_frametable ()

let _ = Callback.register_exception "ml_z_overflow" Overflow

external neg: t -> t = neg@ASM
external add: t -> t -> t = add@ASM
external sub: t -> t -> t = sub@ASM
external mul: t -> t -> t = mul@ASM
external div: t -> t -> t = div@ASM
external cdiv: t -> t -> t = "ml_z_cdiv"
external fdiv: t -> t -> t = "ml_z_fdiv"
external rem: t -> t -> t = rem@ASM
external div_rem: t -> t -> (t * t) = "ml_z_div_rem"
external succ: t -> t = succ@ASM
external pred: t -> t = pred@ASM
external abs: t -> t = abs@ASM
external logand: t -> t -> t = logand@ASM
external logor: t -> t -> t = logor@ASM
external logxor: t -> t -> t = logxor@ASM
external lognot: t -> t = lognot@ASM
external shift_left: t -> int -> t = shift_left@ASM
external shift_right: t -> int -> t = shift_right@ASM
external shift_right_trunc: t -> int -> t = shift_right_trunc@ASM
external of_int: int -> t = "ml_z_of_int" @NOALLOC
external of_int32: int32 -> t = "ml_z_of_int32"
external of_int64: int64 -> t = "ml_z_of_int64"
external of_nativeint: nativeint -> t = "ml_z_of_nativeint"
external of_float: float -> t = "ml_z_of_float"
external to_int: t -> int = "ml_z_to_int"
external to_int32: t -> int32 = "ml_z_to_int32"
external to_int64: t -> int64 = "ml_z_to_int64"
external to_nativeint: t -> nativeint = "ml_z_to_nativeint"
external format: string -> t -> string = "ml_z_format"
external of_substring_base: int -> string -> pos:int -> len:int -> t = "ml_z_of_substring_base"
external compare: t -> t -> int = "ml_z_compare" @NOALLOC
external equal: t -> t -> bool = "ml_z_equal" @NOALLOC
external sign: t -> int = "ml_z_sign" @NOALLOC
external gcd: t -> t -> t = "ml_z_gcd"
external gcdext_intern: t -> t -> (t * t * bool) = "ml_z_gcdext_intern"
external sqrt: t -> t = "ml_z_sqrt"
external sqrt_rem: t -> (t * t) = "ml_z_sqrt_rem"
external numbits: t -> int = "ml_z_numbits" @NOALLOC
external trailing_zeros: t -> int = "ml_z_trailing_zeros" @NOALLOC
external popcount: t -> int = "ml_z_popcount"
external hamdist: t -> t -> int = "ml_z_hamdist"
external size: t -> int = "ml_z_size" @NOALLOC
external fits_int: t -> bool = "ml_z_fits_int" @NOALLOC
external fits_int32: t -> bool = "ml_z_fits_int32" @NOALLOC
external fits_int64: t -> bool = "ml_z_fits_int64" @NOALLOC
external fits_nativeint: t -> bool = "ml_z_fits_nativeint" @NOALLOC
external extract: t -> int -> int -> t = "ml_z_extract"
external powm: t -> t -> t -> t = "ml_z_powm"
external pow: t -> int -> t = "ml_z_pow"
external powm_sec: t -> t -> t -> t = "ml_z_powm_sec"
external divexact: t -> t -> t = divexact@ASM
external root: t -> int -> t = "ml_z_root"
external invert: t -> t -> t = "ml_z_invert"
external perfect_power: t -> bool = "ml_z_perfect_power"
external perfect_square: t -> bool = "ml_z_perfect_square"
external probab_prime: t -> int -> int = "ml_z_probab_prime"
external nextprime: t -> t = "ml_z_nextprime"
external hash: t -> int = "ml_z_hash" @NOALLOC
external to_bits: t -> string = "ml_z_to_bits"
external of_bits: string -> t = "ml_z_of_bits"

let zero = of_int 0
let one = of_int 1
let minus_one = of_int (-1)

let min a b = if compare a b <= 0 then a else b
let max a b = if compare a b >= 0 then a else b

let leq a b = compare a b <= 0
let geq a b = compare a b >= 0
let lt a b = compare a b < 0
let gt a b = compare a b > 0

let to_string = format "%d"

let of_string s = of_substring_base 0 s ~pos:0 ~len:(String.length s)
let of_substring = of_substring_base 0
let of_string_base base s = of_substring_base base s ~pos:0 ~len:(String.length s)

let ediv_rem a b =
  (* we have a = q * b + r, but [Big_int]'s remainder satisfies 0 <= r < |b|,
     while [Z]'s remainder satisfies -|b| < r < |b| and sign(r) = sign(a)
   *)
   let q,r = div_rem a b in
   if sign r >= 0 then (q,r) else
   if sign b >= 0 then (pred q, add r b)
   else (succ q, sub r b)

let ediv a b =
   if sign b >= 0 then fdiv a b else cdiv a b

let erem a b =
   let r = rem a b in
   if sign r >= 0 then r else add r (abs b)

let gcdext u v =
  let g,s,z = gcdext_intern u v in
  if z then g, s, div (sub g (mul u s)) v
  else g, div (sub g (mul v s)) u, s

let lcm u v = 
  let g = gcd u v in
  abs (mul (divexact u g) v)

external testbit_internal: t -> int -> bool = "ml_z_testbit" @NOALLOC
let testbit x n =
  if n >= 0 then testbit_internal x n else invalid_arg "Z.testbit"
(* The test [n >= 0] is done in Caml rather than in the C stub code
   so that the latter raises no exceptions and can be declared "noalloc". *)

let is_odd x = testbit_internal x 0
let is_even x  = not (testbit_internal x 0)

let signed_extract x o l =
  if o < 0 then invalid_arg "Z.signed_extract: negative bit offset";
  if l < 1 then invalid_arg "Z.signed_extract: non-positive bit length";
  if testbit x (o + l - 1)
  then lognot (extract (lognot x) o l)
  else extract x o l

let log2 x =
  if sign x > 0 then (numbits x) - 1 else invalid_arg "Z.log2"
let log2up x =
  if sign x > 0 then numbits (pred x) else invalid_arg "Z.log2up"

(* Consider a real number [r] such that
   - the integral part of [r] is the bigint [x]
   - 2^54 <= |x| < 2^63
   - the fractional part of [r] is 0 if [exact = true], 
     nonzero if [exact = false].
   Then, the following function returns [r] correctly rounded
   according to the current rounding mode of the processor.
   This is an instance of the "round to odd" technique formalized in
   "When double rounding is odd" by S. Boldo and G. Melquiond.
   The claim above is lemma Fappli_IEEE_extra.round_odd_fix
   from the CompCert Coq development. *)

let round_to_float x exact =
  let m = to_int64 x in
  (* Unless the fractional part is exactly 0, round m to an odd integer *)
  let m = if exact then m else Int64.logor m 1L in
  (* Then convert m to float, with the current rounding mode. *)
  Int64.to_float m

let to_float x =
  if Obj.is_int (Obj.repr x) then
    (* Fast path *)
    float_of_int (Obj.magic x : int)
  else begin
    let n = numbits x in
    if n <= 63 then
      Int64.to_float (to_int64 x)
    else begin
      let n = n - 55 in
      (* Extract top 55 bits of x *)
      let top = shift_right x n in
      (* Check if the other bits are all zero *)
      let exact = equal x (shift_left top n) in
      (* Round to float and apply exponent *)
      ldexp (round_to_float top exact) n
    end
  end

let print x = print_string (to_string x)
let output chan x = output_string chan (to_string x)
let sprint () x = to_string x
let bprint b x = Buffer.add_string b (to_string x)
let pp_print f x = Format.pp_print_string f (to_string x)

external (~-): t -> t = neg@ASM
let (~+) x = x
external (+): t -> t -> t = add@ASM
external (-): t -> t -> t = sub@ASM
external ( * ): t -> t -> t = mul@ASM
external (/): t -> t -> t = div@ASM
external (/>): t -> t -> t = "ml_z_cdiv"
external (/<): t -> t -> t = "ml_z_fdiv"
external (/|): t -> t -> t = divexact@ASM
external (mod): t -> t -> t = rem@ASM
external (land): t -> t -> t = logand@ASM
external (lor): t -> t -> t = logor@ASM
external (lxor): t -> t -> t = logxor@ASM
external (~!): t -> t = lognot@ASM
external (lsl): t -> int -> t = shift_left@ASM
external (asr): t -> int -> t = shift_right@ASM
external (~$): int -> t = "ml_z_of_int" @NOALLOC
external ( ** ): t -> int -> t = "ml_z_pow"
let (=) = equal
let (<) = lt
let (>) = gt
let (<=) = leq
let (>=) = geq
let (<>) a b = not (equal a b)

let version = @VERSION