ocaml/otherlibs/num/nat.ml

589 lines
22 KiB
OCaml

(***********************************************************************)
(* *)
(* 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 Library General Public License, with *)
(* the special exception on linking described in file ../../LICENSE. *)
(* *)
(***********************************************************************)
open Int_misc
type nat;;
external create_nat: int -> nat = "create_nat"
external set_to_zero_nat: nat -> int -> int -> unit = "set_to_zero_nat"
external blit_nat: nat -> int -> nat -> int -> int -> unit = "blit_nat"
external set_digit_nat: nat -> int -> int -> unit = "set_digit_nat"
external nth_digit_nat: nat -> int -> int = "nth_digit_nat"
external set_digit_nat_native: nat -> int -> nativeint -> unit
= "set_digit_nat_native"
external nth_digit_nat_native: nat -> int -> nativeint = "nth_digit_nat_native"
external num_digits_nat: nat -> int -> int -> int = "num_digits_nat"
external num_leading_zero_bits_in_digit: nat -> int -> int
= "num_leading_zero_bits_in_digit"
external is_digit_int: nat -> int -> bool = "is_digit_int"
external is_digit_zero: nat -> int -> bool = "is_digit_zero"
external is_digit_normalized: nat -> int -> bool = "is_digit_normalized"
external is_digit_odd: nat -> int -> bool = "is_digit_odd"
external incr_nat: nat -> int -> int -> int -> int = "incr_nat"
external add_nat: nat -> int -> int -> nat -> int -> int -> int -> int
= "add_nat" "add_nat_native"
external complement_nat: nat -> int -> int -> unit = "complement_nat"
external decr_nat: nat -> int -> int -> int -> int = "decr_nat"
external sub_nat: nat -> int -> int -> nat -> int -> int -> int -> int
= "sub_nat" "sub_nat_native"
external mult_digit_nat:
nat -> int -> int -> nat -> int -> int -> nat -> int -> int
= "mult_digit_nat" "mult_digit_nat_native"
external mult_nat:
nat -> int -> int -> nat -> int -> int -> nat -> int -> int -> int
= "mult_nat" "mult_nat_native"
external square_nat: nat -> int -> int -> nat -> int -> int -> int
= "square_nat" "square_nat_native"
external shift_left_nat: nat -> int -> int -> nat -> int -> int -> unit
= "shift_left_nat" "shift_left_nat_native"
external div_digit_nat:
nat -> int -> nat -> int -> nat -> int -> int -> nat -> int -> unit
= "div_digit_nat" "div_digit_nat_native"
external div_nat: nat -> int -> int -> nat -> int -> int -> unit
= "div_nat" "div_nat_native"
external shift_right_nat: nat -> int -> int -> nat -> int -> int -> unit
= "shift_right_nat" "shift_right_nat_native"
external compare_digits_nat: nat -> int -> nat -> int -> int
= "compare_digits_nat"
external compare_nat: nat -> int -> int -> nat -> int -> int -> int
= "compare_nat" "compare_nat_native"
external land_digit_nat: nat -> int -> nat -> int -> unit = "land_digit_nat"
external lor_digit_nat: nat -> int -> nat -> int -> unit = "lor_digit_nat"
external lxor_digit_nat: nat -> int -> nat -> int -> unit = "lxor_digit_nat"
external initialize_nat: unit -> unit = "initialize_nat"
let _ = initialize_nat()
let length_nat (n : nat) = Obj.size (Obj.repr n) - 1
let length_of_digit = Sys.word_size;;
let make_nat len =
if len < 0 then invalid_arg "make_nat" else
let res = create_nat len in set_to_zero_nat res 0 len; res
(* Nat temporaries *)
let a_2 = make_nat 2
and a_1 = make_nat 1
and b_2 = make_nat 2
let copy_nat nat off_set length =
let res = create_nat (length) in
blit_nat res 0 nat off_set length;
res
let is_zero_nat n off len =
compare_nat (make_nat 1) 0 1 n off (num_digits_nat n off len) = 0
let is_nat_int nat off len =
num_digits_nat nat off len = 1 && is_digit_int nat off
let sys_int_of_nat nat off len =
if is_nat_int nat off len
then nth_digit_nat nat off
else failwith "int_of_nat"
let int_of_nat nat =
sys_int_of_nat nat 0 (length_nat nat)
let nat_of_int i =
if i < 0 then invalid_arg "nat_of_int" else
let res = make_nat 1 in
if i = 0 then res else begin set_digit_nat res 0 i; res end
let eq_nat nat1 off1 len1 nat2 off2 len2 =
compare_nat nat1 off1 (num_digits_nat nat1 off1 len1)
nat2 off2 (num_digits_nat nat2 off2 len2) = 0
and le_nat nat1 off1 len1 nat2 off2 len2 =
compare_nat nat1 off1 (num_digits_nat nat1 off1 len1)
nat2 off2 (num_digits_nat nat2 off2 len2) <= 0
and lt_nat nat1 off1 len1 nat2 off2 len2 =
compare_nat nat1 off1 (num_digits_nat nat1 off1 len1)
nat2 off2 (num_digits_nat nat2 off2 len2) < 0
and ge_nat nat1 off1 len1 nat2 off2 len2 =
compare_nat nat1 off1 (num_digits_nat nat1 off1 len1)
nat2 off2 (num_digits_nat nat2 off2 len2) >= 0
and gt_nat nat1 off1 len1 nat2 off2 len2 =
compare_nat nat1 off1 (num_digits_nat nat1 off1 len1)
nat2 off2 (num_digits_nat nat2 off2 len2) > 0
(* XL: now implemented in C for better performance.
The code below doesn't handle carries correctly.
Fortunately, the carry is never used. *)
(***
let square_nat nat1 off1 len1 nat2 off2 len2 =
let c = ref 0
and trash = make_nat 1 in
(* Double product *)
for i = 0 to len2 - 2 do
c := !c + mult_digit_nat
nat1
(succ (off1 + 2 * i))
(2 * (pred (len2 - i)))
nat2
(succ (off2 + i))
(pred (len2 - i))
nat2
(off2 + i)
done;
shift_left_nat nat1 0 len1 trash 0 1;
(* Square of digit *)
for i = 0 to len2 - 1 do
c := !c + mult_digit_nat
nat1
(off1 + 2 * i)
(len1 - 2 * i)
nat2
(off2 + i)
1
nat2
(off2 + i)
done;
!c
***)
let gcd_int_nat i nat off len =
if i = 0 then 1 else
if is_nat_int nat off len then begin
set_digit_nat nat off (gcd_int (nth_digit_nat nat off) i); 0
end else begin
let len_copy = succ len in
let copy = create_nat len_copy
and quotient = create_nat 1
and remainder = create_nat 1 in
blit_nat copy 0 nat off len;
set_digit_nat copy len 0;
div_digit_nat quotient 0 remainder 0 copy 0 len_copy (nat_of_int i) 0;
set_digit_nat nat off (gcd_int (nth_digit_nat remainder 0) i);
0
end
let exchange r1 r2 =
let old1 = !r1 in r1 := !r2; r2 := old1
let gcd_nat nat1 off1 len1 nat2 off2 len2 =
if is_zero_nat nat1 off1 len1 then begin
blit_nat nat1 off1 nat2 off2 len2; len2
end else begin
let copy1 = ref (create_nat (succ len1))
and copy2 = ref (create_nat (succ len2)) in
blit_nat !copy1 0 nat1 off1 len1;
blit_nat !copy2 0 nat2 off2 len2;
set_digit_nat !copy1 len1 0;
set_digit_nat !copy2 len2 0;
if lt_nat !copy1 0 len1 !copy2 0 len2
then exchange copy1 copy2;
let real_len1 =
ref (num_digits_nat !copy1 0 (length_nat !copy1))
and real_len2 =
ref (num_digits_nat !copy2 0 (length_nat !copy2)) in
while not (is_zero_nat !copy2 0 !real_len2) do
set_digit_nat !copy1 !real_len1 0;
div_nat !copy1 0 (succ !real_len1) !copy2 0 !real_len2;
exchange copy1 copy2;
real_len1 := !real_len2;
real_len2 := num_digits_nat !copy2 0 !real_len2
done;
blit_nat nat1 off1 !copy1 0 !real_len1;
!real_len1
end
(* Integer square root using newton method (nearest integer by default) *)
(* Theorem: the sequence x_{n+1} = ( x_n + a/x_n )/2 converges toward
the integer square root (by default) of a for any starting value x_0
strictly greater than the square root of a except if a + 1 is a
perfect square. In this situation, the sequence alternates between
the excess and default integer square root. In any case, the last
strictly decreasing term is the expected result *)
let sqrt_nat rad off len =
let len = num_digits_nat rad off len in
(* Working copy of radicand *)
let len_parity = len mod 2 in
let rad_len = len + 1 + len_parity in
let rad =
let res = create_nat rad_len in
blit_nat res 0 rad off len;
set_digit_nat res len 0;
set_digit_nat res (rad_len - 1) 0;
res in
let cand_len = (len + 1) / 2 in (* ceiling len / 2 *)
let cand_rest = rad_len - cand_len in
(* Candidate square root cand = "|FFFF .... |" *)
let cand = make_nat cand_len in
(* Improve starting square root:
We compute nbb, the number of significant bits of the first digit of the
candidate
(half of the number of significant bits in the first two digits
of the radicand extended to an even length).
shift_cand is word_size - nbb *)
let shift_cand =
((num_leading_zero_bits_in_digit rad (len-1)) +
Sys.word_size * len_parity) / 2 in
(* All radicand bits are zeroed, we give back 0. *)
if shift_cand = Sys.word_size then cand else
begin
complement_nat cand 0 cand_len;
shift_right_nat cand 0 1 a_1 0 shift_cand;
let next_cand = create_nat rad_len in
(* Repeat until *)
let rec loop () =
(* next_cand := rad *)
blit_nat next_cand 0 rad 0 rad_len;
(* next_cand <- next_cand / cand *)
div_nat next_cand 0 rad_len cand 0 cand_len;
(* next_cand (strong weight) <- next_cand (strong weight) + cand,
i.e. next_cand <- cand + rad / cand *)
ignore (add_nat next_cand cand_len cand_rest cand 0 cand_len 0);
(* next_cand <- next_cand / 2 *)
shift_right_nat next_cand cand_len cand_rest a_1 0 1;
if lt_nat next_cand cand_len cand_rest cand 0 cand_len then
begin (* cand <- next_cand *)
blit_nat cand 0 next_cand cand_len cand_len; loop ()
end
else cand in
loop ()
end;;
let power_base_max = make_nat 2;;
match length_of_digit with
| 64 ->
set_digit_nat power_base_max 0 (Int64.to_int 1000000000000000000L);
ignore
(mult_digit_nat power_base_max 0 2
power_base_max 0 1 (nat_of_int 9) 0)
| 32 -> set_digit_nat power_base_max 0 1000000000
| _ -> assert false
;;
let pmax =
match length_of_digit with
| 64 -> 19
| 32 -> 9
| _ -> assert false
;;
let max_superscript_10_power_in_int =
match length_of_digit with
| 64 -> 18
| 32 -> 9
| _ -> assert false
;;
let max_power_10_power_in_int =
match length_of_digit with
| 64 -> nat_of_int (Int64.to_int 1000000000000000000L)
| 32 -> nat_of_int 1000000000
| _ -> assert false
;;
let raw_string_of_digit nat off =
if is_nat_int nat off 1
then begin string_of_int (nth_digit_nat nat off) end
else begin
blit_nat b_2 0 nat off 1;
div_digit_nat a_2 0 a_1 0 b_2 0 2 max_power_10_power_in_int 0;
let leading_digits = nth_digit_nat a_2 0
and s1 = string_of_int (nth_digit_nat a_1 0) in
let len = String.length s1 in
if leading_digits < 10 then begin
let result = Bytes.make (max_superscript_10_power_in_int+1) '0' in
Bytes.set result 0 (Char.chr (48 + leading_digits));
String.blit s1 0 result (Bytes.length result - len) len;
Bytes.to_string result
end else begin
let result = Bytes.make (max_superscript_10_power_in_int+2) '0' in
String.blit (string_of_int leading_digits) 0 result 0 2;
String.blit s1 0 result (Bytes.length result - len) len;
Bytes.to_string result
end
end
(* XL: suppression de string_of_digit et de sys_string_of_digit.
La copie est de toute facon faite dans string_of_nat, qui est le
seul point d entree public dans ce code.
| Deletion of string_of_digit and sys_string_of_digit.
The copy is already done in string_of_nat which is the only
public entry point in this code
*)
(******
let sys_string_of_digit nat off =
let s = raw_string_of_digit nat off in
let result = String.create (String.length s) in
String.blit s 0 result 0 (String.length s);
s
let string_of_digit nat =
sys_string_of_digit nat 0
*******)
let digits = "0123456789ABCDEF"
(*
make_power_base affecte power_base des puissances successives de base a
partir de la puissance 1-ieme.
A la fin de la boucle i-1 est la plus grande puissance de la base qui tient
sur un seul digit et j est la plus grande puissance de la base qui tient
sur un int.
This function returns [(pmax, pint)] where:
[pmax] is the index of the digit of [power_base] that contains the
the maximum power of [base] that fits in a digit. This is also one
less than the exponent of that power.
[pint] is the exponent of the maximum power of [base] that fits in an [int].
*)
let make_power_base base power_base =
let i = ref 0
and j = ref 0 in
set_digit_nat power_base 0 base;
while incr i; is_digit_zero power_base !i do
ignore
(mult_digit_nat power_base !i 2
power_base (pred !i) 1
power_base 0)
done;
while !j < !i - 1 && is_digit_int power_base !j do incr j done;
(!i - 2, !j)
(*
int_to_string places the representation of the integer int in base 'base'
in the string s by starting from the end position pos and going towards
the start, for 'times' places and updates the value of pos.
*)
let int_to_string int s pos_ref base times =
let i = ref int
and j = ref times in
while ((!i != 0) || (!j != 0)) && (!pos_ref != -1) do
Bytes.set s !pos_ref (String.get digits (!i mod base));
decr pos_ref;
decr j;
i := !i / base
done
let power_base_int base i =
if i = 0 || base = 1 then
nat_of_int 1
else if base = 0 then
nat_of_int 0
else if i < 0 then
invalid_arg "power_base_int"
else begin
let power_base = make_nat (succ length_of_digit) in
let (pmax, pint) = make_power_base base power_base in
let n = i / (succ pmax)
and rem = i mod (succ pmax) in
if n > 0 then begin
let newn =
if i = biggest_int then n else (succ n) in
let res = make_nat newn
and res2 = make_nat newn
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 newn in
let len2 = min n (2 * len) in
let succ_len2 = succ len2 in
ignore (square_nat res2 0 len2 res 0 len);
if n land (1 lsl i) > 0 then begin
set_to_zero_nat res 0 len;
ignore
(mult_digit_nat res 0 succ_len2
res2 0 len2 power_base pmax)
end else
blit_nat res 0 res2 0 len2;
set_to_zero_nat res2 0 len2
done;
if rem > 0 then begin
ignore
(mult_digit_nat res2 0 newn
res 0 n power_base (pred rem));
res2
end else res
end else
copy_nat power_base (pred rem) 1
end
(* the ith element (i >= 2) of num_digits_max_vector is :
| |
| biggest_string_length * log (i) |
| ------------------------------- | + 1
| length_of_digit * log (2) |
-- --
*)
(* XL: ai specialise le code d origine a length_of_digit = 32.
| the original code have been specialized to a length_of_digit = 32. *)
(* Now deleted (useless?) *)
(******
let num_digits_max_vector =
[|0; 0; 1024; 1623; 2048; 2378; 2647; 2875; 3072; 3246; 3402;
3543; 3671; 3789; 3899; 4001; 4096|]
let num_digits_max_vector =
match length_of_digit with
16 -> [|0; 0; 2048; 3246; 4096; 4755; 5294; 5749; 6144; 6492; 6803;
7085; 7342; 7578; 7797; 8001; 8192|]
(* If really exotic machines !!!!
| 17 -> [|0; 0; 1928; 3055; 3855; 4476; 4983; 5411; 5783; 6110; 6403;
6668; 6910; 7133; 7339; 7530; 7710|]
| 18 -> [|0; 0; 1821; 2886; 3641; 4227; 4706; 5111; 5461; 5771; 6047;
6298; 6526; 6736; 6931; 7112; 7282|]
| 19 -> [|0; 0; 1725; 2734; 3449; 4005; 4458; 4842; 5174; 5467; 5729;
5966; 6183; 6382; 6566; 6738; 6898|]
| 20 -> [|0; 0; 1639; 2597; 3277; 3804; 4235; 4600; 4915; 5194; 5443;
5668; 5874; 6063; 6238; 6401; 6553|]
| 21 -> [|0; 0; 1561; 2473; 3121; 3623; 4034; 4381; 4681; 4946; 5183;
5398; 5594; 5774; 5941; 6096; 6241|]
| 22 -> [|0; 0; 1490; 2361; 2979; 3459; 3850; 4182; 4468; 4722; 4948;
5153; 5340; 5512; 5671; 5819; 5958|]
| 23 -> [|0; 0; 1425; 2258; 2850; 3308; 3683; 4000; 4274; 4516; 4733;
4929; 5108; 5272; 5424; 5566; 5699|]
| 24 -> [|0; 0; 1366; 2164; 2731; 3170; 3530; 3833; 4096; 4328; 4536;
4723; 4895; 5052; 5198; 5334; 5461|]
| 25 -> [|0; 0; 1311; 2078; 2622; 3044; 3388; 3680; 3932; 4155; 4354;
4534; 4699; 4850; 4990; 5121; 5243|]
| 26 -> [|0; 0; 1261; 1998; 2521; 2927; 3258; 3538; 3781; 3995; 4187;
4360; 4518; 4664; 4798; 4924; 5041|]
| 27 -> [|0; 0; 1214; 1924; 2428; 2818; 3137; 3407; 3641; 3847; 4032;
4199; 4351; 4491; 4621; 4742; 4855|]
| 28 -> [|0; 0; 1171; 1855; 2341; 2718; 3025; 3286; 3511; 3710; 3888;
4049; 4196; 4331; 4456; 4572; 4681|]
| 29 -> [|0; 0; 1130; 1791; 2260; 2624; 2921; 3172; 3390; 3582; 3754;
3909; 4051; 4181; 4302; 4415; 4520|]
| 30 -> [|0; 0; 1093; 1732; 2185; 2536; 2824; 3067; 3277; 3463; 3629;
3779; 3916; 4042; 4159; 4267; 4369|]
| 31 -> [|0; 0; 1057; 1676; 2114; 2455; 2733; 2968; 3171; 3351; 3512;
3657; 3790; 3912; 4025; 4130; 4228|]
*)
| 32 -> [|0; 0; 1024; 1623; 2048; 2378; 2647; 2875; 3072; 3246; 3402;
3543; 3671; 3789; 3899; 4001; 4096|]
| n -> failwith "num_digits_max_vector"
******)
let unadjusted_string_of_nat nat off len_nat =
let len = num_digits_nat nat off len_nat in
if len = 1 then
raw_string_of_digit nat off
else
let len_copy = ref (succ len) in
let copy1 = create_nat !len_copy
and copy2 = make_nat !len_copy
and rest_digit = make_nat 2 in
if len > biggest_int / (succ pmax)
then failwith "number too long"
else let len_s = (succ pmax) * len in
let s = Bytes.make len_s '0'
and pos_ref = ref len_s in
len_copy := pred !len_copy;
blit_nat copy1 0 nat off len;
set_digit_nat copy1 len 0;
while not (is_zero_nat copy1 0 !len_copy) do
div_digit_nat copy2 0
rest_digit 0
copy1 0 (succ !len_copy)
power_base_max 0;
let str = raw_string_of_digit rest_digit 0 in
String.blit str 0
s (!pos_ref - String.length str)
(String.length str);
pos_ref := !pos_ref - pmax;
len_copy := num_digits_nat copy2 0 !len_copy;
blit_nat copy1 0 copy2 0 !len_copy;
set_digit_nat copy1 !len_copy 0
done;
Bytes.unsafe_to_string s
let string_of_nat nat =
let s = unadjusted_string_of_nat nat 0 (length_nat nat)
and index = ref 0 in
begin try
for i = 0 to String.length s - 2 do
if String.get s i <> '0' then (index:= i; raise Exit)
done
with Exit -> ()
end;
String.sub s !index (String.length s - !index)
let base_digit_of_char c base =
let n = Char.code c in
if n >= 48 && n <= 47 + min base 10 then n - 48
else if n >= 65 && n <= 65 + base - 11 then n - 55
else if n >= 97 && n <= 97 + base - 11 then n - 87
else failwith "invalid digit"
(*
The substring (s, off, len) represents a nat in base 'base' which is
determined here
*)
let sys_nat_of_string base s off len =
let power_base = make_nat (succ length_of_digit) in
let (pmax, pint) = make_power_base base power_base in
let new_len = ref (1 + len / (pmax + 1))
and current_len = ref 1 in
let possible_len = ref (min 2 !new_len) in
let nat1 = make_nat !new_len
and nat2 = make_nat !new_len
and digits_read = ref 0
and bound = off + len - 1
and int = ref 0 in
for i = off to bound do
(*
we read (at most) pint digits, we transform them in a int
and integrate it to the number
*)
let c = String.get s i in
begin match c with
' ' | '\t' | '\n' | '\r' | '\\' -> ()
| '_' when i > off -> ()
| _ -> int := !int * base + base_digit_of_char c base;
incr digits_read
end;
if (!digits_read = pint || i = bound) && not (!digits_read = 0) then
begin
set_digit_nat nat1 0 !int;
let erase_len = if !new_len = !current_len then !current_len - 1
else !current_len in
for j = 1 to erase_len do
set_digit_nat nat1 j 0
done;
ignore
(mult_digit_nat nat1 0 !possible_len
nat2 0 !current_len power_base (pred !digits_read));
blit_nat nat2 0 nat1 0 !possible_len;
current_len := num_digits_nat nat1 0 !possible_len;
possible_len := min !new_len (succ !current_len);
int := 0;
digits_read := 0
end
done;
(*
We reframe nat
*)
let nat = create_nat !current_len in
blit_nat nat 0 nat1 0 !current_len;
nat
let nat_of_string s = sys_nat_of_string 10 s 0 (String.length s)
let float_of_nat nat = float_of_string(string_of_nat nat)