diff options
Diffstat (limited to 'otherlibs/num/big_int.ml')
| -rw-r--r-- | otherlibs/num/big_int.ml | 898 |
1 files changed, 0 insertions, 898 deletions
diff --git a/otherlibs/num/big_int.ml b/otherlibs/num/big_int.ml deleted file mode 100644 index 45cea9ca78..0000000000 --- a/otherlibs/num/big_int.ml +++ /dev/null @@ -1,898 +0,0 @@ -(**************************************************************************) -(* *) -(* OCaml *) -(* *) -(* Valerie Menissier-Morain, projet Cristal, INRIA Rocquencourt *) -(* *) -(* Copyright 1996 Institut National de Recherche en Informatique et *) -(* en Automatique. *) -(* *) -(* All rights reserved. This file is distributed under the terms of *) -(* the GNU Lesser General Public License version 2.1, with the *) -(* special exception on linking described in the file LICENSE. *) -(* *) -(**************************************************************************) - -open Int_misc -open Nat - -type big_int = - { sign : int; - abs_value : nat } - -let create_big_int sign nat = - if sign = 1 || sign = -1 || - (sign = 0 && - is_zero_nat nat 0 (num_digits_nat nat 0 (length_nat nat))) - then { sign = sign; - abs_value = nat } - else invalid_arg "create_big_int" - -(* Sign of a big_int *) -let sign_big_int bi = bi.sign - -let zero_big_int = - { sign = 0; - abs_value = make_nat 1 } - -let unit_big_int = - { sign = 1; - abs_value = nat_of_int 1 } - -(* Number of digits in a big_int *) -let num_digits_big_int bi = - num_digits_nat (bi.abs_value) 0 (length_nat bi.abs_value) - -(* Number of bits in a big_int *) -let num_bits_big_int bi = - let nd = num_digits_nat (bi.abs_value) 0 (length_nat bi.abs_value) in - (* nd = 1 if bi = 0 *) - let lz = num_leading_zero_bits_in_digit bi.abs_value (nd - 1) in - (* lz = length_of_digit if bi = 0 *) - nd * length_of_digit - lz - (* = 0 if bi = 0 *) - -(* Opposite of a big_int *) -let minus_big_int bi = - { sign = - bi.sign; - abs_value = copy_nat (bi.abs_value) 0 (num_digits_big_int bi)} - -(* Absolute value of a big_int *) -let abs_big_int bi = - { sign = if bi.sign = 0 then 0 else 1; - abs_value = copy_nat (bi.abs_value) 0 (num_digits_big_int bi)} - -(* Comparison operators on big_int *) - -(* - compare_big_int (bi, bi2) = sign of (bi-bi2) - i.e. 1 if bi > bi2 - 0 if bi = bi2 - -1 if bi < bi2 -*) -let compare_big_int bi1 bi2 = - if bi1.sign = 0 && bi2.sign = 0 then 0 - else if bi1.sign < bi2.sign then -1 - else if bi1.sign > bi2.sign then 1 - else if bi1.sign = 1 then - compare_nat (bi1.abs_value) 0 (num_digits_big_int bi1) - (bi2.abs_value) 0 (num_digits_big_int bi2) - else - compare_nat (bi2.abs_value) 0 (num_digits_big_int bi2) - (bi1.abs_value) 0 (num_digits_big_int bi1) - -let eq_big_int bi1 bi2 = compare_big_int bi1 bi2 = 0 -and le_big_int bi1 bi2 = compare_big_int bi1 bi2 <= 0 -and ge_big_int bi1 bi2 = compare_big_int bi1 bi2 >= 0 -and lt_big_int bi1 bi2 = compare_big_int bi1 bi2 < 0 -and gt_big_int bi1 bi2 = compare_big_int bi1 bi2 > 0 - -let max_big_int bi1 bi2 = if lt_big_int bi1 bi2 then bi2 else bi1 -and min_big_int bi1 bi2 = if gt_big_int bi1 bi2 then bi2 else bi1 - -(* Operations on big_int *) - -let pred_big_int bi = - match bi.sign with - 0 -> { sign = -1; abs_value = nat_of_int 1} - | 1 -> let size_bi = num_digits_big_int bi in - let copy_bi = copy_nat (bi.abs_value) 0 size_bi in - ignore (decr_nat copy_bi 0 size_bi 0); - { sign = if is_zero_nat copy_bi 0 size_bi then 0 else 1; - abs_value = copy_bi } - | _ -> let size_bi = num_digits_big_int bi in - let size_res = succ (size_bi) in - let copy_bi = create_nat (size_res) in - blit_nat copy_bi 0 (bi.abs_value) 0 size_bi; - set_digit_nat copy_bi size_bi 0; - ignore (incr_nat copy_bi 0 size_res 1); - { sign = -1; - abs_value = copy_bi } - -let succ_big_int bi = - match bi.sign with - 0 -> {sign = 1; abs_value = nat_of_int 1} - | -1 -> let size_bi = num_digits_big_int bi in - let copy_bi = copy_nat (bi.abs_value) 0 size_bi in - ignore (decr_nat copy_bi 0 size_bi 0); - { sign = if is_zero_nat copy_bi 0 size_bi then 0 else -1; - abs_value = copy_bi } - | _ -> let size_bi = num_digits_big_int bi in - let size_res = succ (size_bi) in - let copy_bi = create_nat (size_res) in - blit_nat copy_bi 0 (bi.abs_value) 0 size_bi; - set_digit_nat copy_bi size_bi 0; - ignore (incr_nat copy_bi 0 size_res 1); - { sign = 1; - abs_value = copy_bi } - -let add_big_int bi1 bi2 = - let size_bi1 = num_digits_big_int bi1 - and size_bi2 = num_digits_big_int bi2 in - if bi1.sign = bi2.sign - then (* Add absolute values if signs are the same *) - { sign = bi1.sign; - abs_value = - match compare_nat (bi1.abs_value) 0 size_bi1 - (bi2.abs_value) 0 size_bi2 with - -1 -> let res = create_nat (succ size_bi2) in - (blit_nat res 0 (bi2.abs_value) 0 size_bi2; - set_digit_nat res size_bi2 0; - ignore - (add_nat res 0 (succ size_bi2) - (bi1.abs_value) 0 size_bi1 0); - res) - |_ -> let res = create_nat (succ size_bi1) in - (blit_nat res 0 (bi1.abs_value) 0 size_bi1; - set_digit_nat res size_bi1 0; - ignore (add_nat res 0 (succ size_bi1) - (bi2.abs_value) 0 size_bi2 0); - res)} - - else (* Subtract absolute values if signs are different *) - match compare_nat (bi1.abs_value) 0 size_bi1 - (bi2.abs_value) 0 size_bi2 with - 0 -> zero_big_int - | 1 -> { sign = bi1.sign; - abs_value = - let res = copy_nat (bi1.abs_value) 0 size_bi1 in - (ignore (sub_nat res 0 size_bi1 - (bi2.abs_value) 0 size_bi2 1); - res) } - | _ -> { sign = bi2.sign; - abs_value = - let res = copy_nat (bi2.abs_value) 0 size_bi2 in - (ignore (sub_nat res 0 size_bi2 - (bi1.abs_value) 0 size_bi1 1); - res) } - -(* Coercion with int type *) -let big_int_of_int i = - { sign = sign_int i; - abs_value = - let res = (create_nat 1) - in (if i = monster_int - then (set_digit_nat res 0 biggest_int; - ignore (incr_nat res 0 1 1)) - else set_digit_nat res 0 (abs i)); - res } - -let add_int_big_int i bi = add_big_int (big_int_of_int i) bi - -let sub_big_int bi1 bi2 = add_big_int bi1 (minus_big_int bi2) - -(* Returns i * bi *) -let mult_int_big_int i bi = - let size_bi = num_digits_big_int bi in - let size_res = succ size_bi in - if i = monster_int - then let res = create_nat size_res in - blit_nat res 0 (bi.abs_value) 0 size_bi; - set_digit_nat res size_bi 0; - ignore (mult_digit_nat res 0 size_res (bi.abs_value) 0 size_bi - (nat_of_int biggest_int) 0); - { sign = - (sign_big_int bi); - abs_value = res } - else let res = make_nat (size_res) in - ignore (mult_digit_nat res 0 size_res (bi.abs_value) 0 size_bi - (nat_of_int (abs i)) 0); - { sign = (sign_int i) * (sign_big_int bi); - abs_value = res } - -let mult_big_int bi1 bi2 = - let size_bi1 = num_digits_big_int bi1 - and size_bi2 = num_digits_big_int bi2 in - let size_res = size_bi1 + size_bi2 in - let res = make_nat (size_res) in - { sign = bi1.sign * bi2.sign; - abs_value = - if size_bi2 > size_bi1 - then (ignore (mult_nat res 0 size_res (bi2.abs_value) 0 size_bi2 - (bi1.abs_value) 0 size_bi1);res) - else (ignore (mult_nat res 0 size_res (bi1.abs_value) 0 size_bi1 - (bi2.abs_value) 0 size_bi2);res) } - -(* (quotient, remainder ) of the euclidian division of 2 big_int *) -let quomod_big_int bi1 bi2 = - if bi2.sign = 0 then raise Division_by_zero - else - let size_bi1 = num_digits_big_int bi1 - and size_bi2 = num_digits_big_int bi2 in - match compare_nat (bi1.abs_value) 0 size_bi1 - (bi2.abs_value) 0 size_bi2 with - -1 -> (* 1/2 -> 0, remains 1, -1/2 -> -1, remains 1 *) - (* 1/-2 -> 0, remains 1, -1/-2 -> 1, remains 1 *) - if bi1.sign >= 0 then - (big_int_of_int 0, bi1) - else if bi2.sign >= 0 then - (big_int_of_int(-1), add_big_int bi2 bi1) - else - (big_int_of_int 1, sub_big_int bi1 bi2) - | 0 -> (big_int_of_int (bi1.sign * bi2.sign), zero_big_int) - | _ -> let bi1_negatif = bi1.sign = -1 in - let size_q = - if bi1_negatif - then succ (max (succ (size_bi1 - size_bi2)) 1) - else max (succ (size_bi1 - size_bi2)) 1 - and size_r = succ (max size_bi1 size_bi2) - (* r is long enough to contain both quotient and remainder *) - (* of the euclidian division *) - in - (* set up quotient, remainder *) - let q = create_nat size_q - and r = create_nat size_r in - blit_nat r 0 (bi1.abs_value) 0 size_bi1; - set_to_zero_nat r size_bi1 (size_r - size_bi1); - - (* do the division of |bi1| by |bi2| - - at the beginning, r contains |bi1| - - at the end, r contains - * in the size_bi2 least significant digits, the remainder - * in the size_r-size_bi2 most significant digits, the quotient - note the conditions for application of div_nat are verified here - *) - div_nat r 0 size_r (bi2.abs_value) 0 size_bi2; - - (* separate quotient and remainder *) - blit_nat q 0 r size_bi2 (size_r - size_bi2); - let not_null_mod = not (is_zero_nat r 0 size_bi2) in - - (* correct the signs, adjusting the quotient and remainder *) - if bi1_negatif && not_null_mod - then - (* bi1<0, r>0, noting r for (r, size_bi2) the remainder, *) - (* we have |bi1|=q * |bi2| + r with 0 < r < |bi2|, *) - (* thus -bi1 = q * |bi2| + r *) - (* and bi1 = (-q) * |bi2| + (-r) with -|bi2| < (-r) < 0 *) - (* thus bi1 = -(q+1) * |bi2| + (|bi2|-r) *) - (* with 0 < (|bi2|-r) < |bi2| *) - (* so the quotient has for sign the opposite of the bi2'one *) - (* and for value q+1 *) - (* and the remainder is strictly positive *) - (* has for value |bi2|-r *) - (let new_r = copy_nat (bi2.abs_value) 0 size_bi2 in - (* new_r contains (r, size_bi2) the remainder *) - { sign = - bi2.sign; - abs_value = (set_digit_nat q (pred size_q) 0; - ignore (incr_nat q 0 size_q 1); q) }, - { sign = 1; - abs_value = - (ignore (sub_nat new_r 0 size_bi2 r 0 size_bi2 1); - new_r) }) - else - (if bi1_negatif then set_digit_nat q (pred size_q) 0; - { sign = if is_zero_nat q 0 size_q - then 0 - else bi1.sign * bi2.sign; - abs_value = q }, - { sign = if not_null_mod then 1 else 0; - abs_value = copy_nat r 0 size_bi2 }) - -let div_big_int bi1 bi2 = fst (quomod_big_int bi1 bi2) -and mod_big_int bi1 bi2 = snd (quomod_big_int bi1 bi2) - -let gcd_big_int bi1 bi2 = - let size_bi1 = num_digits_big_int bi1 - and size_bi2 = num_digits_big_int bi2 in - if is_zero_nat (bi1.abs_value) 0 size_bi1 then abs_big_int bi2 - else if is_zero_nat (bi2.abs_value) 0 size_bi2 then - { sign = 1; - abs_value = bi1.abs_value } - else - { sign = 1; - abs_value = - match compare_nat (bi1.abs_value) 0 size_bi1 - (bi2.abs_value) 0 size_bi2 with - 0 -> bi1.abs_value - | 1 -> - let res = copy_nat (bi1.abs_value) 0 size_bi1 in - let len = - gcd_nat res 0 size_bi1 (bi2.abs_value) 0 size_bi2 in - copy_nat res 0 len - | _ -> - let res = copy_nat (bi2.abs_value) 0 size_bi2 in - let len = - gcd_nat res 0 size_bi2 (bi1.abs_value) 0 size_bi1 in - copy_nat res 0 len - } - -(* Coercion operators *) - -let monster_big_int = big_int_of_int monster_int;; - -let monster_nat = monster_big_int.abs_value;; - -let is_int_big_int bi = - num_digits_big_int bi == 1 && - match compare_nat bi.abs_value 0 1 monster_nat 0 1 with - | 0 -> bi.sign == -1 - | -1 -> true - | _ -> false;; - -let int_of_big_int bi = - try let n = int_of_nat bi.abs_value in - if bi.sign = -1 then - n else n - with Failure _ -> - if eq_big_int bi monster_big_int then monster_int - else failwith "int_of_big_int";; - -let int_of_big_int_opt bi = - try Some (int_of_big_int bi) with Failure _ -> None - -let big_int_of_nativeint i = - if i = 0n then - zero_big_int - else if i > 0n then begin - let res = create_nat 1 in - set_digit_nat_native res 0 i; - { sign = 1; abs_value = res } - end else begin - let res = create_nat 1 in - set_digit_nat_native res 0 (Nativeint.neg i); - { sign = -1; abs_value = res } - end - -let nativeint_of_big_int bi = - if num_digits_big_int bi > 1 then failwith "nativeint_of_big_int"; - let i = nth_digit_nat_native bi.abs_value 0 in - if bi.sign >= 0 then - if i >= 0n then i else failwith "nativeint_of_big_int" - else - if i >= 0n || i = Nativeint.min_int - then Nativeint.neg i - else failwith "nativeint_of_big_int" - -let nativeint_of_big_int_opt bi = - try Some (nativeint_of_big_int bi) with Failure _ -> None - -let big_int_of_int32 i = big_int_of_nativeint (Nativeint.of_int32 i) - -let int32_of_big_int bi = - let i = nativeint_of_big_int bi in - if i <= 0x7FFF_FFFFn && i >= -0x8000_0000n - then Nativeint.to_int32 i - else failwith "int32_of_big_int" - -let int32_of_big_int_opt bi = - try Some (int32_of_big_int bi) with Failure _ -> None - -let big_int_of_int64 i = - if Sys.word_size = 64 then - big_int_of_nativeint (Int64.to_nativeint i) - else begin - let (sg, absi) = - if i = 0L then (0, 0L) - else if i > 0L then (1, i) - else (-1, Int64.neg i) in - let res = create_nat 2 in - set_digit_nat_native res 0 (Int64.to_nativeint absi); - set_digit_nat_native res 1 (Int64.to_nativeint (Int64.shift_right absi 32)); - { sign = sg; abs_value = res } - end - -let int64_of_big_int bi = - if Sys.word_size = 64 then - Int64.of_nativeint (nativeint_of_big_int bi) - else begin - let i = - match num_digits_big_int bi with - | 1 -> Int64.logand - (Int64.of_nativeint (nth_digit_nat_native bi.abs_value 0)) - 0xFFFFFFFFL - | 2 -> Int64.logor - (Int64.logand - (Int64.of_nativeint (nth_digit_nat_native bi.abs_value 0)) - 0xFFFFFFFFL) - (Int64.shift_left - (Int64.of_nativeint (nth_digit_nat_native bi.abs_value 1)) - 32) - | _ -> failwith "int64_of_big_int" in - if bi.sign >= 0 then - if i >= 0L then i else failwith "int64_of_big_int" - else - if i >= 0L || i = Int64.min_int - then Int64.neg i - else failwith "int64_of_big_int" - end - -let int64_of_big_int_opt bi = - try Some (int64_of_big_int bi) with Failure _ -> None - -(* Coercion with nat type *) -let nat_of_big_int bi = - if bi.sign = -1 - then failwith "nat_of_big_int" - else copy_nat (bi.abs_value) 0 (num_digits_big_int bi) - -let sys_big_int_of_nat nat off len = - let length = num_digits_nat nat off len in - { sign = if is_zero_nat nat off length then 0 else 1; - abs_value = copy_nat nat off length } - -let big_int_of_nat nat = - sys_big_int_of_nat nat 0 (length_nat nat) - -(* Coercion with string type *) - -let string_of_big_int bi = - if bi.sign = -1 - then "-" ^ string_of_nat bi.abs_value - else string_of_nat bi.abs_value - - -let sys_big_int_of_string_aux s ofs len sgn base = - if len < 1 then failwith "sys_big_int_of_string"; - let n = sys_nat_of_string base s ofs len in - if is_zero_nat n 0 (length_nat n) then zero_big_int - else {sign = sgn; abs_value = n} -;; - -let sys_big_int_of_string_base s ofs len sgn = - if len < 1 then failwith "sys_big_int_of_string"; - if len < 2 then sys_big_int_of_string_aux s ofs len sgn 10 - else - match (s.[ofs], s.[ofs+1]) with - | ('0', 'x') | ('0', 'X') -> - sys_big_int_of_string_aux s (ofs+2) (len-2) sgn 16 - | ('0', 'o') | ('0', 'O') -> - sys_big_int_of_string_aux s (ofs+2) (len-2) sgn 8 - | ('0', 'b') | ('0', 'B') -> - sys_big_int_of_string_aux s (ofs+2) (len-2) sgn 2 - | _ -> sys_big_int_of_string_aux s ofs len sgn 10 -;; - -let sys_big_int_of_string s ofs len = - if len < 1 then failwith "sys_big_int_of_string"; - match s.[ofs] with - | '-' -> sys_big_int_of_string_base s (ofs+1) (len-1) (-1) - | '+' -> sys_big_int_of_string_base s (ofs+1) (len-1) 1 - | _ -> sys_big_int_of_string_base s ofs len 1 -;; - -let big_int_of_string s = - sys_big_int_of_string s 0 (String.length s) - -let big_int_of_string_opt s = - try Some (big_int_of_string s) with Failure _ -> None - -let power_base_nat base nat off len = - if base = 0 then nat_of_int 0 else - if is_zero_nat nat off len || base = 1 then nat_of_int 1 else - let power_base = make_nat (succ length_of_digit) in - let (pmax, _pint) = make_power_base base power_base in - let (n, rem) = - let (x, y) = quomod_big_int (sys_big_int_of_nat nat off len) - (big_int_of_int (succ pmax)) in - (int_of_big_int x, int_of_big_int y) in - if n = 0 then copy_nat power_base (pred rem) 1 else - begin - let res = make_nat n - and res2 = make_nat (succ n) - and l = num_bits_int n - 2 in - blit_nat res 0 power_base pmax 1; - for i = l downto 0 do - let len = num_digits_nat res 0 n in - let len2 = min n (2 * len) in - let succ_len2 = succ len2 in - ignore (square_nat res2 0 len2 res 0 len); - begin - if n land (1 lsl i) > 0 - then (set_to_zero_nat res 0 len; - ignore (mult_digit_nat res 0 succ_len2 - res2 0 len2 power_base pmax)) - else blit_nat res 0 res2 0 len2 - end; - set_to_zero_nat res2 0 len2 - done; - if rem > 0 - then (ignore (mult_digit_nat res2 0 (succ n) - res 0 n power_base (pred rem)); - res2) - else res - end - -let power_int_positive_int i n = - match sign_int n with - 0 -> unit_big_int - | -1 -> invalid_arg "power_int_positive_int" - | _ -> let nat = power_base_int (abs i) n in - { sign = if i >= 0 - then sign_int i - else if n land 1 = 0 - then 1 - else -1; - abs_value = nat} - -let power_big_int_positive_int bi n = - match sign_int n with - 0 -> unit_big_int - | -1 -> invalid_arg "power_big_int_positive_int" - | _ -> let bi_len = num_digits_big_int bi in - let res_len = bi_len * n in - let res = make_nat res_len - and res2 = make_nat res_len - and l = num_bits_int n - 2 in - blit_nat res 0 bi.abs_value 0 bi_len; - for i = l downto 0 do - let len = num_digits_nat res 0 res_len in - let len2 = min res_len (2 * len) in - set_to_zero_nat res2 0 len2; - ignore (square_nat res2 0 len2 res 0 len); - if n land (1 lsl i) > 0 then begin - let lenp = min res_len (len2 + bi_len) in - set_to_zero_nat res 0 lenp; - ignore(mult_nat res 0 lenp res2 0 len2 (bi.abs_value) 0 bi_len) - end else begin - blit_nat res 0 res2 0 len2 - end - done; - {sign = if bi.sign >= 0 then bi.sign - else if n land 1 = 0 then 1 else -1; - abs_value = res} - -let power_int_positive_big_int i bi = - match sign_big_int bi with - 0 -> unit_big_int - | -1 -> invalid_arg "power_int_positive_big_int" - | _ -> let nat = power_base_nat - (abs i) (bi.abs_value) 0 (num_digits_big_int bi) in - { sign = if i >= 0 - then sign_int i - else if is_digit_odd (bi.abs_value) 0 - then -1 - else 1; - abs_value = nat } - -let power_big_int_positive_big_int bi1 bi2 = - match sign_big_int bi2 with - 0 -> unit_big_int - | -1 -> invalid_arg "power_big_int_positive_big_int" - | _ -> try - power_big_int_positive_int bi1 (int_of_big_int bi2) - with Failure _ -> - try - power_int_positive_big_int (int_of_big_int bi1) bi2 - with Failure _ -> - raise Out_of_memory - (* If neither bi1 nor bi2 is a small integer, bi1^bi2 is not - representable. Indeed, on a 32-bit platform, - |bi1| >= 2 and |bi2| >= 2^30, hence bi1^bi2 has at least - 2^30 bits = 2^27 bytes, greater than the max size of - allocated blocks. On a 64-bit platform, - |bi1| >= 2 and |bi2| >= 2^62, hence bi1^bi2 has at least - 2^62 bits = 2^59 bytes, greater than the max size of - allocated blocks. *) - -(* base_power_big_int compute bi*base^n *) -let base_power_big_int base n bi = - match sign_int n with - 0 -> bi - | -1 -> let nat = power_base_int base (-n) in - let len_nat = num_digits_nat nat 0 (length_nat nat) - and len_bi = num_digits_big_int bi in - if len_bi < len_nat then - invalid_arg "base_power_big_int" - else if len_bi = len_nat && - compare_digits_nat (bi.abs_value) len_bi nat len_nat = -1 - then invalid_arg "base_power_big_int" - else - let copy = create_nat (succ len_bi) in - blit_nat copy 0 (bi.abs_value) 0 len_bi; - set_digit_nat copy len_bi 0; - div_nat copy 0 (succ len_bi) - nat 0 len_nat; - if not (is_zero_nat copy 0 len_nat) - then invalid_arg "base_power_big_int" - else { sign = bi.sign; - abs_value = copy_nat copy len_nat 1 } - | _ -> let nat = power_base_int base n in - let len_nat = num_digits_nat nat 0 (length_nat nat) - and len_bi = num_digits_big_int bi in - let new_len = len_bi + len_nat in - let res = make_nat new_len in - ignore - (if len_bi > len_nat - then mult_nat res 0 new_len - (bi.abs_value) 0 len_bi - nat 0 len_nat - else mult_nat res 0 new_len - nat 0 len_nat - (bi.abs_value) 0 len_bi) - ; if is_zero_nat res 0 new_len - then zero_big_int - else create_big_int (bi.sign) res - -(* Other functions needed *) - -(* Integer part of the square root of a big_int *) -let sqrt_big_int bi = - match bi.sign with - | 0 -> zero_big_int - | -1 -> invalid_arg "sqrt_big_int" - | _ -> {sign = 1; - abs_value = sqrt_nat (bi.abs_value) 0 (num_digits_big_int bi)} - -let square_big_int bi = - if bi.sign == 0 then zero_big_int else - let len_bi = num_digits_big_int bi in - let len_res = 2 * len_bi in - let res = make_nat len_res in - ignore (square_nat res 0 len_res (bi.abs_value) 0 len_bi); - {sign = 1; abs_value = res} - -(* round off of the futur last digit (of the integer represented by the string - argument of the function) that is now the previous one. - if s contains an integer of the form (10^n)-1 - then s <- only 0 digits and the result_int is true - else s <- the round number and the result_int is false *) -let round_futur_last_digit s off_set length = - let l = pred (length + off_set) in - if Char.code(Bytes.get s l) >= Char.code '5' - then - let rec round_rec l = - if l < off_set then true else begin - let current_char = Bytes.get s l in - if current_char = '9' then - (Bytes.set s l '0'; round_rec (pred l)) - else - (Bytes.set s l (Char.chr (succ (Char.code current_char))); - false) - end - in round_rec (pred l) - else false - - -(* Approximation with floating decimal point a` la approx_ratio_exp *) -let approx_big_int prec bi = - let len_bi = num_digits_big_int bi in - let n = - max 0 - (int_of_big_int ( - add_int_big_int - (-prec) - (div_big_int (mult_big_int (big_int_of_int (pred len_bi)) - (big_int_of_string "963295986")) - (big_int_of_string "100000000")))) in - let s = - Bytes.unsafe_of_string - (string_of_big_int (div_big_int bi (power_int_positive_int 10 n))) - in - let (sign, off) = - if Bytes.get s 0 = '-' - then ("-", 1) - else ("", 0) in - if (round_futur_last_digit s off (succ prec)) - then (sign^"1."^(String.make prec '0')^"e"^ - (string_of_int (n + 1 - off + Bytes.length s))) - else (sign^(Bytes.sub_string s off 1)^"."^ - (Bytes.sub_string s (succ off) (pred prec)) - ^"e"^(string_of_int (n - succ off + Bytes.length s))) - -(* Logical operations *) - -(* Shift left by N bits *) - -let shift_left_big_int bi n = - if n < 0 then invalid_arg "shift_left_big_int" - else if n = 0 then bi - else if bi.sign = 0 then bi - else begin - let size_bi = num_digits_big_int bi in - let size_res = size_bi + ((n + length_of_digit - 1) / length_of_digit) in - let res = create_nat size_res in - let ndigits = n / length_of_digit in - set_to_zero_nat res 0 ndigits; - blit_nat res ndigits bi.abs_value 0 size_bi; - let nbits = n mod length_of_digit in - if nbits > 0 then - shift_left_nat res ndigits size_bi res (ndigits + size_bi) nbits; - { sign = bi.sign; abs_value = res } - end - -(* Shift right by N bits (rounds toward zero) *) - -let shift_right_towards_zero_big_int bi n = - if n < 0 then invalid_arg "shift_right_towards_zero_big_int" - else if n = 0 then bi - else if bi.sign = 0 then bi - else begin - let size_bi = num_digits_big_int bi in - let ndigits = n / length_of_digit in - let nbits = n mod length_of_digit in - if ndigits >= size_bi then zero_big_int else begin - let size_res = size_bi - ndigits in - let res = create_nat size_res in - blit_nat res 0 bi.abs_value ndigits size_res; - if nbits > 0 then begin - let tmp = create_nat 1 in - shift_right_nat res 0 size_res tmp 0 nbits - end; - if is_zero_nat res 0 size_res - then zero_big_int - else { sign = bi.sign; abs_value = res } - end - end - -(* Compute 2^n - 1 *) - -let two_power_m1_big_int n = - if n < 0 then invalid_arg "two_power_m1_big_int" - else if n = 0 then zero_big_int - else begin - let idx = n / length_of_digit in - let size_res = idx + 1 in - let res = make_nat size_res in - set_digit_nat_native res idx - (Nativeint.shift_left 1n (n mod length_of_digit)); - ignore (decr_nat res 0 size_res 0); - { sign = 1; abs_value = res } - end - -(* Shift right by N bits (rounds toward minus infinity) *) - -let shift_right_big_int bi n = - if n < 0 then invalid_arg "shift_right_big_int" - else if bi.sign >= 0 then shift_right_towards_zero_big_int bi n - else - shift_right_towards_zero_big_int (sub_big_int bi (two_power_m1_big_int n)) n - -(* Extract N bits starting at ofs. - Treats bi in two's complement. - Result is always positive. *) - -let extract_big_int bi ofs n = - if ofs < 0 || n < 0 then invalid_arg "extract_big_int" - else if bi.sign = 0 then bi - else begin - let size_bi = num_digits_big_int bi in - let size_res = (n + length_of_digit - 1) / length_of_digit in - let ndigits = ofs / length_of_digit in - let nbits = ofs mod length_of_digit in - let res = make_nat size_res in - if ndigits < size_bi then - blit_nat res 0 bi.abs_value ndigits (min size_res (size_bi - ndigits)); - if bi.sign < 0 then begin - (* Two's complement *) - complement_nat res 0 size_res; - (* PR#6010: need to increment res iff digits 0...ndigits-1 of bi are 0. - In this case, digits 0...ndigits-1 of not(bi) are all 0xFF...FF, - and adding 1 to them produces a carry out at ndigits. *) - let rec carry_incr i = - i >= ndigits || i >= size_bi || - (is_digit_zero bi.abs_value i && carry_incr (i + 1)) in - if carry_incr 0 then ignore (incr_nat res 0 size_res 1) - end; - if nbits > 0 then begin - let tmp = create_nat 1 in - shift_right_nat res 0 size_res tmp 0 nbits - end; - let n' = n mod length_of_digit in - if n' > 0 then begin - let tmp = create_nat 1 in - set_digit_nat_native tmp 0 - (Nativeint.shift_right_logical (-1n) (length_of_digit - n')); - land_digit_nat res (size_res - 1) tmp 0 - end; - if is_zero_nat res 0 size_res - then zero_big_int - else { sign = 1; abs_value = res } - end - -(* Bitwise logical operations. Arguments must be >= 0. *) - -let and_big_int a b = - if a.sign < 0 || b.sign < 0 then invalid_arg "and_big_int" - else if a.sign = 0 || b.sign = 0 then zero_big_int - else begin - let size_a = num_digits_big_int a - and size_b = num_digits_big_int b in - let size_res = min size_a size_b in - let res = create_nat size_res in - blit_nat res 0 a.abs_value 0 size_res; - for i = 0 to size_res - 1 do - land_digit_nat res i b.abs_value i - done; - if is_zero_nat res 0 size_res - then zero_big_int - else { sign = 1; abs_value = res } - end - -let or_big_int a b = - if a.sign < 0 || b.sign < 0 then invalid_arg "or_big_int" - else if a.sign = 0 then b - else if b.sign = 0 then a - else begin - let size_a = num_digits_big_int a - and size_b = num_digits_big_int b in - let size_res = max size_a size_b in - let res = create_nat size_res in - let or_aux a' b' size_b' = - blit_nat res 0 a'.abs_value 0 size_res; - for i = 0 to size_b' - 1 do - lor_digit_nat res i b'.abs_value i - done in - if size_a >= size_b - then or_aux a b size_b - else or_aux b a size_a; - if is_zero_nat res 0 size_res - then zero_big_int - else { sign = 1; abs_value = res } - end - -let xor_big_int a b = - if a.sign < 0 || b.sign < 0 then invalid_arg "xor_big_int" - else if a.sign = 0 then b - else if b.sign = 0 then a - else begin - let size_a = num_digits_big_int a - and size_b = num_digits_big_int b in - let size_res = max size_a size_b in - let res = create_nat size_res in - let xor_aux a' b' size_b' = - blit_nat res 0 a'.abs_value 0 size_res; - for i = 0 to size_b' - 1 do - lxor_digit_nat res i b'.abs_value i - done in - if size_a >= size_b - then xor_aux a b size_b - else xor_aux b a size_a; - if is_zero_nat res 0 size_res - then zero_big_int - else { sign = 1; abs_value = res } - end - -(* Coercion with float type *) - -(* 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 to - the nearest double-precision floating-point number. - 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_big_int_to_float x exact = - assert (let n = num_bits_big_int x in 55 <= n && n <= 63); - let m = int64_of_big_int 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 normal rounding mode. *) - Int64.to_float m - -let float_of_big_int x = - let n = num_bits_big_int x in - if n <= 63 then - Int64.to_float (int64_of_big_int x) - else begin - let n = n - 55 in - (* Extract top 55 bits of x *) - let top = shift_right_big_int x n in - (* Check if the other bits are all zero *) - let exact = eq_big_int x (shift_left_big_int top n) in - (* Round to float and apply exponent *) - ldexp (round_big_int_to_float top exact) n - end |