diff options
Diffstat (limited to 'otherlibs')
28 files changed, 0 insertions, 5637 deletions
diff --git a/otherlibs/num/.depend b/otherlibs/num/.depend deleted file mode 100644 index dfd9e70fc3..0000000000 --- a/otherlibs/num/.depend +++ /dev/null @@ -1,40 +0,0 @@ -bng.$(O): bng.c bng.h ../../byterun/caml/config.h ../../byterun/caml/m.h \ - ../../byterun/caml/s.h bng_digit.c -bng_amd64.$(O): bng_amd64.c -bng_arm64.$(O): bng_arm64.c -bng_digit.$(O): bng_digit.c -bng_ia32.$(O): bng_ia32.c -bng_ppc.$(O): bng_ppc.c -bng_sparc.$(O): bng_sparc.c -nat_stubs.$(O): nat_stubs.c ../../byterun/caml/alloc.h \ - ../../byterun/caml/misc.h ../../byterun/caml/config.h \ - ../../byterun/caml/m.h ../../byterun/caml/s.h \ - ../../byterun/caml/mlvalues.h ../../byterun/caml/config.h \ - ../../byterun/caml/custom.h ../../byterun/caml/intext.h \ - ../../byterun/caml/io.h ../../byterun/caml/fail.h \ - ../../byterun/caml/hash.h ../../byterun/caml/memory.h \ - ../../byterun/caml/gc.h ../../byterun/caml/major_gc.h \ - ../../byterun/caml/freelist.h ../../byterun/caml/minor_gc.h \ - ../../byterun/caml/address_class.h ../../byterun/caml/mlvalues.h bng.h \ - nat.h -arith_flags.cmo : arith_flags.cmi -arith_flags.cmx : arith_flags.cmi -arith_flags.cmi : -arith_status.cmo : arith_flags.cmi arith_status.cmi -arith_status.cmx : arith_flags.cmx arith_status.cmi -arith_status.cmi : -big_int.cmo : nat.cmi int_misc.cmi big_int.cmi -big_int.cmx : nat.cmx int_misc.cmx big_int.cmi -big_int.cmi : nat.cmi -int_misc.cmo : int_misc.cmi -int_misc.cmx : int_misc.cmi -int_misc.cmi : -nat.cmo : int_misc.cmi nat.cmi -nat.cmx : int_misc.cmx nat.cmi -nat.cmi : -num.cmo : ratio.cmi nat.cmi int_misc.cmi big_int.cmi arith_flags.cmi num.cmi -num.cmx : ratio.cmx nat.cmx int_misc.cmx big_int.cmx arith_flags.cmx num.cmi -num.cmi : ratio.cmi nat.cmi big_int.cmi -ratio.cmo : nat.cmi int_misc.cmi big_int.cmi arith_flags.cmi ratio.cmi -ratio.cmx : nat.cmx int_misc.cmx big_int.cmx arith_flags.cmx ratio.cmi -ratio.cmi : nat.cmi big_int.cmi diff --git a/otherlibs/num/Makefile b/otherlibs/num/Makefile deleted file mode 100644 index a43b58a004..0000000000 --- a/otherlibs/num/Makefile +++ /dev/null @@ -1,42 +0,0 @@ -#************************************************************************** -#* * -#* OCaml * -#* * -#* Xavier Leroy, projet Cristal, INRIA Rocquencourt * -#* * -#* Copyright 1999 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. * -#* * -#************************************************************************** - -# Makefile for the "num" (exact rational arithmetic) library - -LIBNAME=nums -EXTRACFLAGS=-DBNG_ARCH_$(BNG_ARCH) -DBNG_ASM_LEVEL=$(BNG_ASM_LEVEL) -CAMLOBJS=int_misc.cmo nat.cmo big_int.cmo arith_flags.cmo \ - ratio.cmo num.cmo arith_status.cmo -CMIFILES=big_int.cmi nat.cmi num.cmi ratio.cmi arith_status.cmi -COBJS=bng.$(O) nat_stubs.$(O) - -include ../Makefile - -clean:: - rm -f *~ - -bng.$(O): bng.h bng_digit.c \ - bng_amd64.c bng_ia32.c bng_ppc.c bng_sparc.c - -.PHONY: depend -depend: -ifeq "$(TOOLCHAIN)" "msvc" - $(error Dependencies cannot be regenerated using the MSVC ports) -else - $(CPP) -MM $(CPPFLAGS) *.c | sed -e 's/\.o/.$$(O)/g' > .depend - $(CAMLRUN) $(ROOTDIR)/tools/ocamldep -slash *.mli *.ml >> .depend -endif - -include .depend diff --git a/otherlibs/num/Makefile.nt b/otherlibs/num/Makefile.nt deleted file mode 100644 index ed9900bb9a..0000000000 --- a/otherlibs/num/Makefile.nt +++ /dev/null @@ -1,16 +0,0 @@ -#************************************************************************** -#* * -#* OCaml * -#* * -#* Xavier Leroy, projet Cristal, INRIA Rocquencourt * -#* * -#* Copyright 1999 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. * -#* * -#************************************************************************** - -include Makefile diff --git a/otherlibs/num/README b/otherlibs/num/README deleted file mode 100644 index a979356d2b..0000000000 --- a/otherlibs/num/README +++ /dev/null @@ -1,55 +0,0 @@ -The "libnum" library implements exact-precision arithmetic on -big integers and on rationals. - -This library is derived from Valerie Menissie-Morain's implementation -of rational arithmetic for Caml V3.1 (INRIA). Xavier Leroy (INRIA) -did the Caml Light port. Victor Manuel Gulias Fernandez did the -initial Caml Special Light port. Pierre Weis did most of the -maintenance and bug fixing. - -Initially, the low-level big integer operations were provided by the -BigNum package developed by Bernard Serpette, Jean Vuillemin and -Jean-Claude Herve (INRIA and Digital PRL). License issues forced us to -replace the BigNum package. The current implementation of low-level -big integer operations is due to Xavier Leroy. - -This library is documented in "The CAML Numbers Reference Manual" by -Valerie Menissier-Morain, technical report 141, INRIA, july 1992, -available at ftp://ftp.inria.fr/INRIA/publication/RT/RT-0141.ps.gz - - -USAGE: - -To use the bignum library from your programs, just do - - ocamlc <options> nums.cma <.cmo and .ml files> -or - ocamlopt <options> nums.cmxa <.cmx and .ml files> - -for the linking phase. - -If you'd like to have the bignum functions available at toplevel, do - - ocamlmktop -o ocamltopnum <options> nums.cma <.cmo and .ml files> - ./ocamltopnum - -As an example, try: - - open Num;; - let rec fact n = - if n = 0 then Int 1 else mult_num (num_of_int n) (fact(n-1));; - string_of_num(fact 1000);; - - -PROCESSOR-SPECIFIC OPTIMIZATIONS: - -When compiled with GCC, the low-level primitives use "inline extended asm" -to exploit useful features of the target processor (additions and -subtractions with carry; double-width multiplication, division). -Here are the processors for which such optimizations are available: - IA32 (x86) (carry, dwmult, dwdiv, 64-bit ops with SSE2 if available) - AMD64 (Opteron) (carry, dwmult, dwdiv) - PowerPC (carry, dwmult) - Alpha (dwmult) - SPARC (carry, dwmult, dwdiv) - MIPS (dwmult) diff --git a/otherlibs/num/arith_flags.ml b/otherlibs/num/arith_flags.ml deleted file mode 100644 index a1ca0b028e..0000000000 --- a/otherlibs/num/arith_flags.ml +++ /dev/null @@ -1,24 +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. *) -(* *) -(**************************************************************************) - -let error_when_null_denominator_flag = ref true;; - -let normalize_ratio_flag = ref false;; - -let normalize_ratio_when_printing_flag = ref true;; - -let floating_precision = ref 12;; - -let approx_printing_flag = ref false;; diff --git a/otherlibs/num/arith_flags.mli b/otherlibs/num/arith_flags.mli deleted file mode 100644 index 7dd6bc79d0..0000000000 --- a/otherlibs/num/arith_flags.mli +++ /dev/null @@ -1,20 +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. *) -(* *) -(**************************************************************************) - -val error_when_null_denominator_flag : bool ref -val normalize_ratio_flag : bool ref -val normalize_ratio_when_printing_flag : bool ref -val floating_precision : int ref -val approx_printing_flag : bool ref diff --git a/otherlibs/num/arith_status.ml b/otherlibs/num/arith_status.ml deleted file mode 100644 index 2fbdd4a6a4..0000000000 --- a/otherlibs/num/arith_status.ml +++ /dev/null @@ -1,100 +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 Arith_flags;; - -let get_error_when_null_denominator () = - !error_when_null_denominator_flag -and set_error_when_null_denominator choice = - error_when_null_denominator_flag := choice;; - -let get_normalize_ratio () = !normalize_ratio_flag -and set_normalize_ratio choice = normalize_ratio_flag := choice;; - -let get_normalize_ratio_when_printing () = - !normalize_ratio_when_printing_flag -and set_normalize_ratio_when_printing choice = - normalize_ratio_when_printing_flag := choice;; - -let get_floating_precision () = !floating_precision -and set_floating_precision i = floating_precision := i;; - -let get_approx_printing () = !approx_printing_flag -and set_approx_printing b = approx_printing_flag := b;; - -let arith_print_string s = print_string s; print_string " --> ";; - -let arith_print_bool = function - true -> print_string "ON" -| _ -> print_string "OFF" -;; - -let arith_status () = - print_newline (); - - arith_print_string - "Normalization during computation"; - arith_print_bool (get_normalize_ratio ()); - print_newline (); - print_string " (returned by get_normalize_ratio ())"; - print_newline (); - print_string " (modifiable with set_normalize_ratio <your choice>)"; - print_newline (); - print_newline (); - - arith_print_string - "Normalization when printing"; - arith_print_bool (get_normalize_ratio_when_printing ()); - print_newline (); - print_string - " (returned by get_normalize_ratio_when_printing ())"; - print_newline (); - print_string - " (modifiable with set_normalize_ratio_when_printing <your choice>)"; - print_newline (); - print_newline (); - - arith_print_string - "Floating point approximation when printing rational numbers"; - arith_print_bool (get_approx_printing ()); - print_newline (); - print_string - " (returned by get_approx_printing ())"; - print_newline (); - print_string - " (modifiable with set_approx_printing <your choice>)"; - print_newline (); - (if (get_approx_printing ()) - then (print_string " Default precision = "; - print_int (get_floating_precision ()); - print_newline (); - print_string " (returned by get_floating_precision ())"; - print_newline (); - print_string - " (modifiable with set_floating_precision <your choice>)"; - print_newline (); - print_newline ()) - else print_newline()); - - arith_print_string - "Error when a rational denominator is null"; - arith_print_bool (get_error_when_null_denominator ()); - print_newline (); - print_string " (returned by get_error_when_null_denominator ())"; - print_newline (); - print_string - " (modifiable with set_error_when_null_denominator <your choice>)"; - print_newline () -;; diff --git a/otherlibs/num/arith_status.mli b/otherlibs/num/arith_status.mli deleted file mode 100644 index ba604347fb..0000000000 --- a/otherlibs/num/arith_status.mli +++ /dev/null @@ -1,64 +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. *) -(* *) -(**************************************************************************) - -(** Flags that control rational arithmetic. *) - -val arith_status: unit -> unit - (** Print the current status of the arithmetic flags. *) - -val get_error_when_null_denominator : unit -> bool -(** See {!Arith_status.set_error_when_null_denominator}.*) - -val set_error_when_null_denominator : bool -> unit - (** Get or set the flag [null_denominator]. When on, attempting to - create a rational with a null denominator raises an exception. - When off, rationals with null denominators are accepted. - Initially: on. *) - -val get_normalize_ratio : unit -> bool -(** See {!Arith_status.set_normalize_ratio}.*) - -val set_normalize_ratio : bool -> unit - (** Get or set the flag [normalize_ratio]. When on, rational - numbers are normalized after each operation. When off, - rational numbers are not normalized until printed. - Initially: off. *) - -val get_normalize_ratio_when_printing : unit -> bool -(** See {!Arith_status.set_normalize_ratio_when_printing}.*) - -val set_normalize_ratio_when_printing : bool -> unit - (** Get or set the flag [normalize_ratio_when_printing]. - When on, rational numbers are normalized before being printed. - When off, rational numbers are printed as is, without normalization. - Initially: on. *) - -val get_approx_printing : unit -> bool -(** See {!Arith_status.set_approx_printing}.*) - -val set_approx_printing : bool -> unit - (** Get or set the flag [approx_printing]. - When on, rational numbers are printed as a decimal approximation. - When off, rational numbers are printed as a fraction. - Initially: off. *) - -val get_floating_precision : unit -> int -(** See {!Arith_status.set_floating_precision}.*) - -val set_floating_precision : int -> unit - (** Get or set the parameter [floating_precision]. - This parameter is the number of digits displayed when - [approx_printing] is on. - Initially: 12. *) 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 diff --git a/otherlibs/num/big_int.mli b/otherlibs/num/big_int.mli deleted file mode 100644 index 07c4072955..0000000000 --- a/otherlibs/num/big_int.mli +++ /dev/null @@ -1,276 +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. *) -(* *) -(**************************************************************************) - -(** Operations on arbitrary-precision integers. - - Big integers (type [big_int]) are signed integers of arbitrary size. -*) - -open Nat - -type big_int - (** The type of big integers. *) - -val zero_big_int : big_int -(** The big integer [0]. *) - -val unit_big_int : big_int - (** The big integer [1]. *) - -(** {6 Arithmetic operations} *) - -val minus_big_int : big_int -> big_int -(** Unary negation. *) - -val abs_big_int : big_int -> big_int -(** Absolute value. *) - -val add_big_int : big_int -> big_int -> big_int -(** Addition. *) - -val succ_big_int : big_int -> big_int -(** Successor (add 1). *) - -val add_int_big_int : int -> big_int -> big_int -(** Addition of a small integer to a big integer. *) - -val sub_big_int : big_int -> big_int -> big_int -(** Subtraction. *) - -val pred_big_int : big_int -> big_int -(** Predecessor (subtract 1). *) - -val mult_big_int : big_int -> big_int -> big_int -(** Multiplication of two big integers. *) - -val mult_int_big_int : int -> big_int -> big_int -(** Multiplication of a big integer by a small integer *) - -val square_big_int: big_int -> big_int -(** Return the square of the given big integer *) - -val sqrt_big_int: big_int -> big_int - (** [sqrt_big_int a] returns the integer square root of [a], - that is, the largest big integer [r] such that [r * r <= a]. - Raise [Invalid_argument] if [a] is negative. *) - -val quomod_big_int : big_int -> big_int -> big_int * big_int - (** Euclidean division of two big integers. - The first part of the result is the quotient, - the second part is the remainder. - Writing [(q,r) = quomod_big_int a b], we have - [a = q * b + r] and [0 <= r < |b|]. - Raise [Division_by_zero] if the divisor is zero. *) - -val div_big_int : big_int -> big_int -> big_int - (** Euclidean quotient of two big integers. - This is the first result [q] of [quomod_big_int] (see above). *) - -val mod_big_int : big_int -> big_int -> big_int - (** Euclidean modulus of two big integers. - This is the second result [r] of [quomod_big_int] (see above). *) - -val gcd_big_int : big_int -> big_int -> big_int -(** Greatest common divisor of two big integers. *) - -val power_int_positive_int: int -> int -> big_int -val power_big_int_positive_int: big_int -> int -> big_int -val power_int_positive_big_int: int -> big_int -> big_int -val power_big_int_positive_big_int: big_int -> big_int -> big_int - (** Exponentiation functions. Return the big integer - representing the first argument [a] raised to the power [b] - (the second argument). Depending - on the function, [a] and [b] can be either small integers - or big integers. Raise [Invalid_argument] if [b] is negative. *) - -(** {6 Comparisons and tests} *) - -val sign_big_int : big_int -> int - (** Return [0] if the given big integer is zero, - [1] if it is positive, and [-1] if it is negative. *) - -val compare_big_int : big_int -> big_int -> int - (** [compare_big_int a b] returns [0] if [a] and [b] are equal, - [1] if [a] is greater than [b], and [-1] if [a] is smaller - than [b]. *) - -val eq_big_int : big_int -> big_int -> bool -val le_big_int : big_int -> big_int -> bool -val ge_big_int : big_int -> big_int -> bool -val lt_big_int : big_int -> big_int -> bool -val gt_big_int : big_int -> big_int -> bool -(** Usual boolean comparisons between two big integers. *) - -val max_big_int : big_int -> big_int -> big_int -(** Return the greater of its two arguments. *) - -val min_big_int : big_int -> big_int -> big_int -(** Return the smaller of its two arguments. *) - -val num_digits_big_int : big_int -> int - (** Return the number of machine words used to store the - given big integer. *) - -val num_bits_big_int : big_int -> int - (** Return the number of significant bits in the absolute - value of the given big integer. [num_bits_big_int a] - returns 0 if [a] is 0; otherwise it returns a positive - integer [n] such that [2^(n-1) <= |a| < 2^n]. - - @since 4.03.0 *) - -(** {6 Conversions to and from strings} *) - -val string_of_big_int : big_int -> string - (** Return the string representation of the given big integer, - in decimal (base 10). *) - -val big_int_of_string : string -> big_int - (** Convert a string to a big integer, in decimal. - The string consists of an optional [-] or [+] sign, - followed by one or several decimal digits. *) -(* TODO: document error condition. *) - -val big_int_of_string_opt: string -> big_int option -(** Convert a string to a big integer, in decimal. - The string consists of an optional [-] or [+] sign, - followed by one or several decimal digits. Other the function - returns [None]. - @since 4.05 -*) - - -(** {6 Conversions to and from other numerical types} *) - -val big_int_of_int : int -> big_int -(** Convert a small integer to a big integer. *) - -val is_int_big_int : big_int -> bool - (** Test whether the given big integer is small enough to - be representable as a small integer (type [int]) - without loss of precision. On a 32-bit platform, - [is_int_big_int a] returns [true] if and only if - [a] is between 2{^30} and 2{^30}-1. On a 64-bit platform, - [is_int_big_int a] returns [true] if and only if - [a] is between -2{^62} and 2{^62}-1. *) - -val int_of_big_int : big_int -> int - (** Convert a big integer to a small integer (type [int]). - Raises [Failure "int_of_big_int"] if the big integer - is not representable as a small integer. *) - -val int_of_big_int_opt: big_int -> int option -(** Convert a big integer to a small integer (type [int]). Return - [None] if the big integer is not representable as a small - integer. - @since 4.05 -*) - -val big_int_of_int32 : int32 -> big_int -(** Convert a 32-bit integer to a big integer. *) - -val big_int_of_nativeint : nativeint -> big_int -(** Convert a native integer to a big integer. *) - -val big_int_of_int64 : int64 -> big_int -(** Convert a 64-bit integer to a big integer. *) - -val int32_of_big_int : big_int -> int32 - (** Convert a big integer to a 32-bit integer. - Raises [Failure] if the big integer is outside the - range \[-2{^31}, 2{^31}-1\]. *) - -val int32_of_big_int_opt: big_int -> int32 option -(** Convert a big integer to a 32-bit integer. Return [None] if the - big integer is outside the range \[-2{^31}, 2{^31}-1\]. - @since 4.05 -*) - -val nativeint_of_big_int : big_int -> nativeint - (** Convert a big integer to a native integer. - Raises [Failure] if the big integer is outside the - range [[Nativeint.min_int, Nativeint.max_int]]. *) - -val nativeint_of_big_int_opt: big_int -> nativeint option -(** Convert a big integer to a native integer. Return [None] if the - big integer is outside the range [[Nativeint.min_int, - Nativeint.max_int]]; - @since 4.05 -*) - -val int64_of_big_int : big_int -> int64 - (** Convert a big integer to a 64-bit integer. - Raises [Failure] if the big integer is outside the - range \[-2{^63}, 2{^63}-1\]. *) - -val int64_of_big_int_opt: big_int -> int64 option -(** Convert a big integer to a 64-bit integer. Return [None] if the - big integer is outside the range \[-2{^63}, 2{^63}-1\]. - @since 4.05 -*) - -val float_of_big_int : big_int -> float - (** Returns a floating-point number approximating the - given big integer. *) - -(** {6 Bit-oriented operations} *) - -val and_big_int : big_int -> big_int -> big_int - (** Bitwise logical 'and'. - The arguments must be positive or zero. *) - -val or_big_int : big_int -> big_int -> big_int - (** Bitwise logical 'or'. - The arguments must be positive or zero. *) - -val xor_big_int : big_int -> big_int -> big_int - (** Bitwise logical 'exclusive or'. - The arguments must be positive or zero. *) - -val shift_left_big_int : big_int -> int -> big_int - (** [shift_left_big_int b n] returns [b] shifted left by [n] bits. - Equivalent to multiplication by 2^n. *) - -val shift_right_big_int : big_int -> int -> big_int - (** [shift_right_big_int b n] returns [b] shifted right by [n] bits. - Equivalent to division by 2^n with the result being - rounded towards minus infinity. *) - -val shift_right_towards_zero_big_int : big_int -> int -> big_int - (** [shift_right_towards_zero_big_int b n] returns [b] shifted - right by [n] bits. The shift is performed on the absolute - value of [b], and the result has the same sign as [b]. - Equivalent to division by 2^n with the result being - rounded towards zero. *) - -val extract_big_int : big_int -> int -> int -> big_int - (** [extract_big_int bi ofs n] returns a nonnegative number - corresponding to bits [ofs] to [ofs + n - 1] of the - binary representation of [bi]. If [bi] is negative, - a two's complement representation is used. *) - -(**/**) - -(** {6 For internal use} *) - -val nat_of_big_int : big_int -> nat -val big_int_of_nat : nat -> big_int -val base_power_big_int: int -> int -> big_int -> big_int -val sys_big_int_of_string: string -> int -> int -> big_int -val round_futur_last_digit : bytes -> int -> int -> bool -val approx_big_int: int -> big_int -> string - -val round_big_int_to_float: big_int -> bool -> float -(** @since 4.03.0 *) diff --git a/otherlibs/num/bng.c b/otherlibs/num/bng.c deleted file mode 100644 index c4d0ea1a51..0000000000 --- a/otherlibs/num/bng.c +++ /dev/null @@ -1,433 +0,0 @@ -/**************************************************************************/ -/* */ -/* OCaml */ -/* */ -/* Xavier Leroy, projet Cristal, INRIA Rocquencourt */ -/* */ -/* Copyright 2003 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. */ -/* */ -/**************************************************************************/ - -#include "bng.h" -#include "caml/config.h" - -#if defined(__GNUC__) && BNG_ASM_LEVEL > 0 -#if defined(BNG_ARCH_ia32) -#include "bng_ia32.c" -#elif defined(BNG_ARCH_amd64) -#include "bng_amd64.c" -#elif defined(BNG_ARCH_ppc) -#include "bng_ppc.c" -#elif defined (BNG_ARCH_sparc) -#include "bng_sparc.c" -#elif defined (BNG_ARCH_arm64) -#include "bng_arm64.c" -#endif -#endif - -#include "bng_digit.c" - -/**** Operations that cannot be overridden ****/ - -/* Return number of leading zero bits in d */ -int bng_leading_zero_bits(bngdigit d) -{ - int n = BNG_BITS_PER_DIGIT; -#ifdef ARCH_SIXTYFOUR - if ((d & 0xFFFFFFFF00000000L) != 0) { n -= 32; d = d >> 32; } -#endif - if ((d & 0xFFFF0000) != 0) { n -= 16; d = d >> 16; } - if ((d & 0xFF00) != 0) { n -= 8; d = d >> 8; } - if ((d & 0xF0) != 0) { n -= 4; d = d >> 4; } - if ((d & 0xC) != 0) { n -= 2; d = d >> 2; } - if ((d & 2) != 0) { n -= 1; d = d >> 1; } - return n - d; -} - -/* Complement the digits of {a,len} */ -void bng_complement(bng a/*[alen]*/, bngsize alen) -{ - for (/**/; alen > 0; alen--, a++) *a = ~*a; -} - -/* Return number of significant digits in {a,alen}. */ -bngsize bng_num_digits(bng a/*[alen]*/, bngsize alen) -{ - while (1) { - if (alen == 0) return 1; - if (a[alen - 1] != 0) return alen; - alen--; - } -} - -/* Return 0 if {a,alen} = {b,blen} - -1 if {a,alen} < {b,blen} - 1 if {a,alen} > {b,blen}. */ -int bng_compare(bng a/*[alen]*/, bngsize alen, - bng b/*[blen]*/, bngsize blen) -{ - bngdigit da, db; - - while (alen > 0 && a[alen-1] == 0) alen--; - while (blen > 0 && b[blen-1] == 0) blen--; - if (alen > blen) return 1; - if (alen < blen) return -1; - while (alen > 0) { - alen--; - da = a[alen]; - db = b[alen]; - if (da > db) return 1; - if (da < db) return -1; - } - return 0; -} - -/**** Generic definitions of the overridable operations ****/ - -/* {a,alen} := {a, alen} + carry. Return carry out. */ -static bngcarry bng_generic_add_carry - (bng a/*[alen]*/, bngsize alen, bngcarry carry) -{ - if (carry == 0 || alen == 0) return carry; - do { - if (++(*a) != 0) return 0; - a++; - } while (--alen); - return 1; -} - -/* {a,alen} := {a,alen} + {b,blen} + carry. Return carry out. - Require alen >= blen. */ -static bngcarry bng_generic_add - (bng a/*[alen]*/, bngsize alen, - bng b/*[blen]*/, bngsize blen, - bngcarry carry) -{ - alen -= blen; - for (/**/; blen > 0; blen--, a++, b++) { - BngAdd2Carry(*a, carry, *a, *b, carry); - } - if (carry == 0 || alen == 0) return carry; - do { - if (++(*a) != 0) return 0; - a++; - } while (--alen); - return 1; -} - -/* {a,alen} := {a, alen} - carry. Return carry out. */ -static bngcarry bng_generic_sub_carry - (bng a/*[alen]*/, bngsize alen, bngcarry carry) -{ - if (carry == 0 || alen == 0) return carry; - do { - if ((*a)-- != 0) return 0; - a++; - } while (--alen); - return 1; -} - -/* {a,alen} := {a,alen} - {b,blen} - carry. Return carry out. - Require alen >= blen. */ -static bngcarry bng_generic_sub - (bng a/*[alen]*/, bngsize alen, - bng b/*[blen]*/, bngsize blen, - bngcarry carry) -{ - alen -= blen; - for (/**/; blen > 0; blen--, a++, b++) { - BngSub2Carry(*a, carry, *a, *b, carry); - } - if (carry == 0 || alen == 0) return carry; - do { - if ((*a)-- != 0) return 0; - a++; - } while (--alen); - return 1; -} - -/* {a,alen} := {a,alen} << shift. - Return the bits shifted out of the most significant digit of a. - Require 0 <= shift < BITS_PER_BNGDIGIT. */ -static bngdigit bng_generic_shift_left - (bng a/*[alen]*/, bngsize alen, - int shift) -{ - int shift2 = BNG_BITS_PER_DIGIT - shift; - bngdigit carry = 0; - if (shift > 0) { - for (/**/; alen > 0; alen--, a++) { - bngdigit d = *a; - *a = (d << shift) | carry; - carry = d >> shift2; - } - } - return carry; -} - -/* {a,alen} := {a,alen} >> shift. - Return the bits shifted out of the least significant digit of a. - Require 0 <= shift < BITS_PER_BNGDIGIT. */ -static bngdigit bng_generic_shift_right - (bng a/*[alen]*/, bngsize alen, - int shift) -{ - int shift2 = BNG_BITS_PER_DIGIT - shift; - bngdigit carry = 0; - if (shift > 0) { - for (a = a + alen - 1; alen > 0; alen--, a--) { - bngdigit d = *a; - *a = (d >> shift) | carry; - carry = d << shift2; - } - } - return carry; -} - -/* {a,alen} := {a,alen} + d * {b,blen}. Return carry out. - Require alen >= blen. */ -static bngdigit bng_generic_mult_add_digit - (bng a/*[alen]*/, bngsize alen, - bng b/*[blen]*/, bngsize blen, - bngdigit d) -{ - bngdigit out, ph, pl; - bngcarry carry; - - alen -= blen; - for (out = 0; blen > 0; blen--, a++, b++) { - bngdigit bd = *b; - /* ph:pl = double-digit product of b's current digit and d */ - BngMult(ph, pl, bd, d); - /* current digit of a += pl + out. Accumulate carries in ph. */ - BngAdd3(*a, ph, *a, pl, out); - /* prepare out for next iteration */ - out = ph; - } - if (alen == 0) return out; - /* current digit of a += out */ - BngAdd2(*a, carry, *a, out); - a++; - alen--; - /* Propagate carry */ - if (carry == 0 || alen == 0) return carry; - do { - if (++(*a) != 0) return 0; - a++; - } while (--alen); - return 1; -} - -/* {a,alen} := {a,alen} - d * {b,blen}. Return carry out. - Require alen >= blen. */ -static bngdigit bng_generic_mult_sub_digit - (bng a/*[alen]*/, bngsize alen, - bng b/*[blen]*/, bngsize blen, - bngdigit d) -{ - bngdigit out, ph, pl; - bngcarry carry; - - alen -= blen; - for (out = 0; blen > 0; blen--, a++, b++) { - bngdigit bd = *b; - /* ph:pl = double-digit product of b's current digit and d */ - BngMult(ph, pl, bd, d); - /* current digit of a -= pl + out. Accumulate carrys in ph. */ - BngSub3(*a, ph, *a, pl, out); - /* prepare out for next iteration */ - out = ph; - } - if (alen == 0) return out; - /* current digit of a -= out */ - BngSub2(*a, carry, *a, out); - a++; - alen--; - /* Propagate carry */ - if (carry == 0 || alen == 0) return carry; - do { - if ((*a)-- != 0) return 0; - a++; - } while (--alen); - return 1; -} - -/* {a,alen} := {a,alen} + {b,blen} * {c,clen}. Return carry out. - Require alen >= blen + clen. */ -static bngcarry bng_generic_mult_add - (bng a/*[alen]*/, bngsize alen, - bng b/*[blen]*/, bngsize blen, - bng c/*[clen]*/, bngsize clen) -{ - bngcarry carry; - for (carry = 0; clen > 0; clen--, c++, alen--, a++) - carry += bng_mult_add_digit(a, alen, b, blen, *c); - return carry; -} - -/* {a,alen} := 2 * {a,alen} + {b,blen}^2. Return carry out. - Require alen >= 2 * blen. */ -static bngcarry bng_generic_square_add - (bng a/*[alen]*/, bngsize alen, - bng b/*[blen]*/, bngsize blen) -{ - bngcarry carry1, carry2; - bngsize i, aofs; - bngdigit ph, pl, d; - - /* Double products */ - for (carry1 = 0, i = 1; i < blen; i++) { - aofs = 2 * i - 1; - carry1 += bng_mult_add_digit(a + aofs, alen - aofs, - b + i, blen - i, b[i - 1]); - } - /* Multiply by two */ - carry1 = (carry1 << 1) | bng_shift_left(a, alen, 1); - /* Add square of digits */ - carry2 = 0; - for (i = 0; i < blen; i++) { - d = b[i]; - BngMult(ph, pl, d, d); - BngAdd2Carry(*a, carry2, *a, pl, carry2); - a++; - BngAdd2Carry(*a, carry2, *a, ph, carry2); - a++; - } - alen -= 2 * blen; - if (alen > 0 && carry2 != 0) { - do { - if (++(*a) != 0) { carry2 = 0; break; } - a++; - } while (--alen); - } - return carry1 + carry2; -} - -/* {a,len-1} := {b,len} / d. Return {b,len} modulo d. - Require MSD of b < d. - If BngDivNeedsNormalization is defined, require d normalized. */ -static bngdigit bng_generic_div_rem_norm_digit - (bng a/*[len-1]*/, bng b/*[len]*/, bngsize len, bngdigit d) -{ - bngdigit topdigit, quo, rem; - intnat i; - - topdigit = b[len - 1]; - for (i = len - 2; i >= 0; i--) { - /* Divide topdigit:current digit of numerator by d */ - BngDiv(quo, rem, topdigit, b[i], d); - /* Quotient is current digit of result */ - a[i] = quo; - /* Iterate with topdigit = remainder */ - topdigit = rem; - } - return topdigit; -} - -#ifdef BngDivNeedsNormalization -/* {a,len-1} := {b,len} / d. Return {b,len} modulo d. - Require MSD of b < d. */ -static bngdigit bng_generic_div_rem_digit - (bng a/*[len-1]*/, bng b/*[len]*/, bngsize len, bngdigit d) -{ - bngdigit rem; - int shift; - - /* Normalize d and b */ - shift = bng_leading_zero_bits(d); - d <<= shift; - bng_shift_left(b, len, shift); - /* Do the division */ - rem = bng_div_rem_norm_digit(a, b, len, d); - /* Undo normalization on b and remainder */ - bng_shift_right(b, len, shift); - return rem >> shift; -} -#endif - -/* {n+dlen, nlen-dlen} := {n,nlen} / {d, dlen}. - {n, dlen} := {n,nlen} modulo {d, dlen}. - Require nlen > dlen and MSD of n < MSD of d. - (This implies MSD of d > 0). */ -static void bng_generic_div_rem - (bng n/*[nlen]*/, bngsize nlen, - bng d/*[dlen]*/, bngsize dlen) -{ - bngdigit topden, quo, rem; - int shift; - bngsize i, j; - - /* Normalize d */ - shift = bng_leading_zero_bits(d[dlen - 1]); - /* Note that no bits of n are lost by the following shift, - since n[nlen-1] < d[dlen-1] */ - bng_shift_left(n, nlen, shift); - bng_shift_left(d, dlen, shift); - /* Special case if d is just one digit */ - if (dlen == 1) { - *n = bng_div_rem_norm_digit(n + 1, n, nlen, *d); - } else { - topden = d[dlen - 1]; - /* Long division */ - for (j = nlen - 1; j >= dlen; j--) { - i = j - dlen; - /* At this point: - - the current numerator is n[j] : ...................... : n[0] - - to be subtracted quo times: d[dlen-1] : ... : d[0] : 0... : 0 - (there are i zeroes at the end) */ - /* Under-estimate the next digit of the quotient (quo) */ - if (topden + 1 == 0) - quo = n[j]; - else - BngDiv(quo, rem, n[j], n[j - 1], topden + 1); - /* Subtract d * quo (shifted i places) from numerator */ - n[j] -= bng_mult_sub_digit(n + i, dlen, d, dlen, quo); - /* Adjust if necessary */ - while (n[j] != 0 || bng_compare(n + i, dlen, d, dlen) >= 0) { - /* Numerator is still bigger than shifted divisor. - Increment quotient and subtract shifted divisor. */ - quo++; - n[j] -= bng_sub(n + i, dlen, d, dlen, 0); - } - /* Store quotient digit */ - n[j] = quo; - } - } - /* Undo normalization on remainder and divisor */ - bng_shift_right(n, dlen, shift); - bng_shift_right(d, dlen, shift); -} - -/**** Construction of the table of operations ****/ - -struct bng_operations bng_ops = { - bng_generic_add_carry, - bng_generic_add, - bng_generic_sub_carry, - bng_generic_sub, - bng_generic_shift_left, - bng_generic_shift_right, - bng_generic_mult_add_digit, - bng_generic_mult_sub_digit, - bng_generic_mult_add, - bng_generic_square_add, - bng_generic_div_rem_norm_digit, -#ifdef BngDivNeedsNormalization - bng_generic_div_rem_digit, -#else - bng_generic_div_rem_norm_digit, -#endif - bng_generic_div_rem -}; - -void bng_init(void) -{ -#ifdef BNG_SETUP_OPS - BNG_SETUP_OPS; -#endif -} diff --git a/otherlibs/num/bng.h b/otherlibs/num/bng.h deleted file mode 100644 index 406117dd75..0000000000 --- a/otherlibs/num/bng.h +++ /dev/null @@ -1,156 +0,0 @@ -/**************************************************************************/ -/* */ -/* OCaml */ -/* */ -/* Xavier Leroy, projet Cristal, INRIA Rocquencourt */ -/* */ -/* Copyright 2003 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. */ -/* */ -/**************************************************************************/ - -#include <string.h> -#include "caml/config.h" - -typedef uintnat bngdigit; -typedef bngdigit * bng; -typedef unsigned int bngcarry; -typedef uintnat bngsize; - -#define BNG_BITS_PER_DIGIT (sizeof(bngdigit) * 8) -#define BNG_BITS_PER_HALF_DIGIT (sizeof(bngdigit) * 4) - -struct bng_operations { - - /* {a,alen} := {a, alen} + carry. Return carry out. */ - bngcarry (*add_carry) - (bng a/*[alen]*/, bngsize alen, bngcarry carry); -#define bng_add_carry bng_ops.add_carry - - /* {a,alen} := {a,alen} + {b,blen} + carry. Return carry out. - Require alen >= blen. */ - bngcarry (*add) - (bng a/*[alen]*/, bngsize alen, - bng b/*[blen]*/, bngsize blen, - bngcarry carry); -#define bng_add bng_ops.add - - /* {a,alen} := {a, alen} - carry. Return carry out. */ - bngcarry (*sub_carry) - (bng a/*[alen]*/, bngsize alen, bngcarry carry); -#define bng_sub_carry bng_ops.sub_carry - - /* {a,alen} := {a,alen} - {b,blen} - carry. Return carry out. - Require alen >= blen. */ - bngcarry (*sub) - (bng a/*[alen]*/, bngsize alen, - bng b/*[blen]*/, bngsize blen, - bngcarry carry); -#define bng_sub bng_ops.sub - - /* {a,alen} := {a,alen} << shift. - Return the bits shifted out of the most significant digit of a. - Require 0 <= shift < BITS_PER_BNGDIGIT. */ - bngdigit (*shift_left) - (bng a/*[alen]*/, bngsize alen, - int shift); -#define bng_shift_left bng_ops.shift_left - - /* {a,alen} := {a,alen} >> shift. - Return the bits shifted out of the least significant digit of a. - Require 0 <= shift < BITS_PER_BNGDIGIT. */ - bngdigit (*shift_right) - (bng a/*[alen]*/, bngsize alen, - int shift); -#define bng_shift_right bng_ops.shift_right - - /* {a,alen} := {a,alen} + d * {b,blen}. Return carry out. - Require alen >= blen. - If alen > blen, the carry out returned is 0 or 1. - If alen == blen, the carry out returned is a full digit. */ - bngdigit (*mult_add_digit) - (bng a/*[alen]*/, bngsize alen, - bng b/*[blen]*/, bngsize blen, - bngdigit d); -#define bng_mult_add_digit bng_ops.mult_add_digit - - /* {a,alen} := {a,alen} - d * {b,blen}. Return carry out. - Require alen >= blen. - If alen > blen, the carry out returned is 0 or 1. - If alen == blen, the carry out returned is a full digit. */ - bngdigit (*mult_sub_digit) - (bng a/*[alen]*/, bngsize alen, - bng b/*[blen]*/, bngsize blen, - bngdigit d); -#define bng_mult_sub_digit bng_ops.mult_sub_digit - - /* {a,alen} := {a,alen} + {b,blen} * {c,clen}. Return carry out. - Require alen >= blen + clen. */ - bngcarry (*mult_add) - (bng a/*[alen]*/, bngsize alen, - bng b/*[blen]*/, bngsize blen, - bng c/*[clen]*/, bngsize clen); -#define bng_mult_add bng_ops.mult_add - - /* {a,alen} := 2 * {a,alen} + {b,blen}^2. Return carry out. - Require alen >= 2 * blen. */ - bngcarry (*square_add) - (bng a/*[alen]*/, bngsize alen, - bng b/*[blen]*/, bngsize blen); -#define bng_square_add bng_ops.square_add - - /* {a,len-1} := {b,len} / d. Return {b,len} modulo d. - Require d is normalized and MSD of b < d. - See div_rem_digit for a function that does not require d - to be normalized */ - bngdigit (*div_rem_norm_digit) - (bng a/*[len-1]*/, bng b/*[len]*/, bngsize len, bngdigit d); -#define bng_div_rem_norm_digit bng_ops.div_rem_norm_digit - - /* {a,len-1} := {b,len} / d. Return {b,len} modulo d. - Require MSD of b < d. */ - bngdigit (*div_rem_digit) - (bng a/*[len-1]*/, bng b/*[len]*/, bngsize len, bngdigit d); -#define bng_div_rem_digit bng_ops.div_rem_digit - - /* {n+dlen, nlen-dlen} := {n,nlen} / {d, dlen}. - {n, dlen} := {n,nlen} modulo {d, dlen}. - Require nlen > dlen and MSD of n < MSD of d (which implies d != 0). */ - void (*div_rem) - (bng n/*[nlen]*/, bngsize nlen, - bng d/*[nlen]*/, bngsize dlen); -#define bng_div_rem bng_ops.div_rem -}; - -extern struct bng_operations bng_ops; - -/* Initialize the BNG library */ -extern void bng_init(void); - -/* {a,alen} := 0 */ -#define bng_zero(a,alen) memset((a), 0, (alen) * sizeof(bngdigit)) - -/* {a,len} := {b,len} */ -#define bng_assign(a,b,len) memmove((a), (b), (len) * sizeof(bngdigit)) - -/* Complement the digits of {a,len} */ -extern void bng_complement(bng a/*[alen]*/, bngsize alen); - -/* Return number of significant digits in {a,alen}. */ -extern bngsize bng_num_digits(bng a/*[alen]*/, bngsize alen); - -/* Return 1 if {a,alen} is 0, 0 otherwise. */ -#define bng_is_zero(a,alen) (bng_num_digits(a,alen) == 0) - -/* Return 0 if {a,alen} = {b,blen} - <0 if {a,alen} < {b,blen} - >0 if {a,alen} > {b,blen}. */ -extern int bng_compare(bng a/*[alen]*/, bngsize alen, - bng b/*[blen]*/, bngsize blen); - -/* Return the number of leading zero bits in digit d. */ -extern int bng_leading_zero_bits(bngdigit d); diff --git a/otherlibs/num/bng_amd64.c b/otherlibs/num/bng_amd64.c deleted file mode 100644 index 585900e9c9..0000000000 --- a/otherlibs/num/bng_amd64.c +++ /dev/null @@ -1,195 +0,0 @@ -/**************************************************************************/ -/* */ -/* OCaml */ -/* */ -/* Xavier Leroy, projet Cristal, INRIA Rocquencourt */ -/* */ -/* Copyright 2003 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. */ -/* */ -/**************************************************************************/ - -/* Code specific to the AMD x86_64 architecture. */ - -#define BngAdd2(res,carryout,arg1,arg2) \ - asm("xorl %1, %1 \n\t" \ - "addq %3, %0 \n\t" \ - "setc %b1" \ - : "=r" (res), "=&q" (carryout) \ - : "0" (arg1), "rm" (arg2)) - -#define BngSub2(res,carryout,arg1,arg2) \ - asm("xorl %1, %1 \n\t" \ - "subq %3, %0 \n\t" \ - "setc %b1" \ - : "=r" (res), "=&q" (carryout) \ - : "0" (arg1), "rm" (arg2)) - -#define BngMult(resh,resl,arg1,arg2) \ - asm("mulq %3" \ - : "=a" (resl), "=d" (resh) \ - : "a" (arg1), "r" (arg2)) - -#define BngDiv(quo,rem,nh,nl,d) \ - asm("divq %4" \ - : "=a" (quo), "=d" (rem) \ - : "a" (nl), "d" (nh), "r" (d)) - -/* Reimplementation in asm of some of the bng operations. */ - -static bngcarry bng_amd64_add - (bng a/*[alen]*/, bngsize alen, - bng b/*[blen]*/, bngsize blen, - bngcarry carry) -{ - bngdigit tmp; - alen -= blen; - if (blen > 0) { - asm("negb %b3 \n\t" - "1: \n\t" - "movq (%0), %4 \n\t" - "adcq (%1), %4 \n\t" - "movq %4, (%0) \n\t" - "leaq 8(%0), %0 \n\t" - "leaq 8(%1), %1 \n\t" - "decq %2 \n\t" - "jnz 1b \n\t" - "setc %b3" - : "=r" (a), "=r" (b), "=r" (blen), "=q" (carry), "=r" (tmp) - : "0" (a), "1" (b), "2" (blen), "3" (carry)); - } - if (carry == 0 || alen == 0) return carry; - do { - if (++(*a) != 0) return 0; - a++; - } while (--alen); - return 1; -} - -static bngcarry bng_amd64_sub - (bng a/*[alen]*/, bngsize alen, - bng b/*[blen]*/, bngsize blen, - bngcarry carry) -{ - bngdigit tmp; - alen -= blen; - if (blen > 0) { - asm("negb %b3 \n\t" - "1: \n\t" - "movq (%0), %4 \n\t" - "sbbq (%1), %4 \n\t" - "movq %4, (%0) \n\t" - "leaq 8(%0), %0 \n\t" - "leaq 8(%1), %1 \n\t" - "decq %2 \n\t" - "jnz 1b \n\t" - "setc %b3" - : "=r" (a), "=r" (b), "=r" (blen), "=q" (carry), "=r" (tmp) - : "0" (a), "1" (b), "2" (blen), "3" (carry)); - } - if (carry == 0 || alen == 0) return carry; - do { - if ((*a)-- != 0) return 0; - a++; - } while (--alen); - return 1; -} - -static bngdigit bng_amd64_mult_add_digit - (bng a/*[alen]*/, bngsize alen, - bng b/*[blen]*/, bngsize blen, - bngdigit d) -{ - bngdigit out; - bngcarry carry; - - alen -= blen; - out = 0; - if (blen > 0) { - asm("1: \n\t" - "movq (%1), %%rax \n\t" - "mulq %7\n\t" /* rdx:rax = d * next digit of b */ - "addq (%0), %%rax \n\t" /* add next digit of a to rax */ - "adcq $0, %%rdx \n\t" /* accumulate carry in rdx */ - "addq %3, %%rax \n\t" /* add out to rax */ - "adcq $0, %%rdx \n\t" /* accumulate carry in rdx */ - "movq %%rax, (%0) \n\t" /* rax is next digit of result */ - "movq %%rdx, %3 \n\t" /* rdx is next out */ - "leaq 8(%0), %0 \n\t" - "leaq 8(%1), %1 \n\t" - "decq %2 \n\t" - "jnz 1b" - : "=&r" (a), "=&r" (b), "=&r" (blen), "=&r" (out) - : "0" (a), "1" (b), "2" (blen), "rm" (d), "3" (out) - : "rax", "rdx"); - } - if (alen == 0) return out; - /* current digit of a += out */ - BngAdd2(*a, carry, *a, out); - a++; - alen--; - /* Propagate carry */ - if (carry == 0 || alen == 0) return carry; - do { - if (++(*a) != 0) return 0; - a++; - } while (--alen); - return 1; -} - -static bngdigit bng_amd64_mult_sub_digit - (bng a/*[alen]*/, bngsize alen, - bng b/*[blen]*/, bngsize blen, - bngdigit d) -{ - bngdigit out, tmp; - bngcarry carry; - - alen -= blen; - out = 0; - if (blen > 0) { - asm("1: \n\t" - "movq (%1), %%rax \n\t" - "movq (%0), %4 \n\t" - "mulq %8\n\t" /* rdx:rax = d * next digit of b */ - "subq %%rax, %4 \n\t" /* subtract rax from next digit of a */ - "adcq $0, %%rdx \n\t" /* accumulate carry in rdx */ - "subq %3, %4 \n\t" /* subtract out */ - "adcq $0, %%rdx \n\t" /* accumulate carry in rdx */ - "movq %4, (%0) \n\t" /* store next digit of result */ - "movq %%rdx, %3 \n\t" /* rdx is next out */ - "leaq 8(%0), %0 \n\t" - "leaq 8(%1), %1 \n\t" - "decq %2 \n\t" - "jnz 1b" - : "=&r" (a), "=&r" (b), "=&rm" (blen), "=&r" (out), "=&r" (tmp) - : "0" (a), "1" (b), "2" (blen), "rm" (d), "3" (out) - : "rax", "rdx"); - } - if (alen == 0) return out; - /* current digit of a -= out */ - BngSub2(*a, carry, *a, out); - a++; - alen--; - /* Propagate carry */ - if (carry == 0 || alen == 0) return carry; - do { - if ((*a)-- != 0) return 0; - a++; - } while (--alen); - return 1; -} - -static void bng_amd64_setup_ops(void) -{ - bng_ops.add = bng_amd64_add; - bng_ops.sub = bng_amd64_sub; - bng_ops.mult_add_digit = bng_amd64_mult_add_digit; - bng_ops.mult_sub_digit = bng_amd64_mult_sub_digit; -} - -#define BNG_SETUP_OPS bng_amd64_setup_ops() diff --git a/otherlibs/num/bng_arm64.c b/otherlibs/num/bng_arm64.c deleted file mode 100644 index b900b803b1..0000000000 --- a/otherlibs/num/bng_arm64.c +++ /dev/null @@ -1,22 +0,0 @@ -/**************************************************************************/ -/* */ -/* OCaml */ -/* */ -/* Xavier Leroy, projet Gallium, INRIA Rocquencourt */ -/* */ -/* Copyright 2013 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. */ -/* */ -/**************************************************************************/ - -/* Code specific for the ARM 64 (AArch64) architecture */ - -#define BngMult(resh,resl,arg1,arg2) \ - asm("mul %0, %2, %3 \n\t" \ - "umulh %1, %2, %3" \ - : "=&r" (resl), "=&r" (resh) \ - : "r" (arg1), "r" (arg2)) diff --git a/otherlibs/num/bng_digit.c b/otherlibs/num/bng_digit.c deleted file mode 100644 index 6983af65b0..0000000000 --- a/otherlibs/num/bng_digit.c +++ /dev/null @@ -1,178 +0,0 @@ -/**************************************************************************/ -/* */ -/* OCaml */ -/* */ -/* Xavier Leroy, projet Cristal, INRIA Rocquencourt */ -/* */ -/* Copyright 2003 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. */ -/* */ -/**************************************************************************/ - -/**** Generic operations on digits ****/ - -/* These macros can be defined in the machine-specific include file. - Below are the default definitions (in plain C). - Except for BngMult, all macros are guaranteed to evaluate their - arguments exactly once. */ - -#ifndef BngAdd2 -/* res = arg1 + arg2. carryout = carry out. */ -#define BngAdd2(res,carryout,arg1,arg2) { \ - bngdigit tmp1, tmp2; \ - tmp1 = arg1; \ - tmp2 = tmp1 + (arg2); \ - carryout = (tmp2 < tmp1); \ - res = tmp2; \ -} -#endif - -#ifndef BngAdd2Carry -/* res = arg1 + arg2 + carryin. carryout = carry out. */ -#define BngAdd2Carry(res,carryout,arg1,arg2,carryin) { \ - bngdigit tmp1, tmp2, tmp3; \ - tmp1 = arg1; \ - tmp2 = tmp1 + (arg2); \ - tmp3 = tmp2 + (carryin); \ - carryout = (tmp2 < tmp1) + (tmp3 < tmp2); \ - res = tmp3; \ -} -#endif - -#ifndef BngAdd3 -/* res = arg1 + arg2 + arg3. Each carry increments carryaccu. */ -#define BngAdd3(res,carryaccu,arg1,arg2,arg3) { \ - bngdigit tmp1, tmp2, tmp3; \ - tmp1 = arg1; \ - tmp2 = tmp1 + (arg2); \ - carryaccu += (tmp2 < tmp1); \ - tmp3 = tmp2 + (arg3); \ - carryaccu += (tmp3 < tmp2); \ - res = tmp3; \ -} -#endif - -#ifndef BngSub2 -/* res = arg1 - arg2. carryout = carry out. */ -#define BngSub2(res,carryout,arg1,arg2) { \ - bngdigit tmp1, tmp2; \ - tmp1 = arg1; \ - tmp2 = arg2; \ - res = tmp1 - tmp2; \ - carryout = (tmp1 < tmp2); \ -} -#endif - -#ifndef BngSub2Carry -/* res = arg1 - arg2 - carryin. carryout = carry out. */ -#define BngSub2Carry(res,carryout,arg1,arg2,carryin) { \ - bngdigit tmp1, tmp2, tmp3; \ - tmp1 = arg1; \ - tmp2 = arg2; \ - tmp3 = tmp1 - tmp2; \ - res = tmp3 - (carryin); \ - carryout = (tmp1 < tmp2) + (tmp3 < carryin); \ -} -#endif - -#ifndef BngSub3 -/* res = arg1 - arg2 - arg3. Each carry increments carryaccu. */ -#define BngSub3(res,carryaccu,arg1,arg2,arg3) { \ - bngdigit tmp1, tmp2, tmp3, tmp4; \ - tmp1 = arg1; \ - tmp2 = arg2; \ - tmp3 = arg3; \ - tmp4 = tmp1 - tmp2; \ - res = tmp4 - tmp3; \ - carryaccu += (tmp1 < tmp2) + (tmp4 < tmp3); \ -} -#endif - -#define BngLowHalf(d) ((d) & (((bngdigit)1 << BNG_BITS_PER_HALF_DIGIT) - 1)) -#define BngHighHalf(d) ((d) >> BNG_BITS_PER_HALF_DIGIT) - -#ifndef BngMult -/* resl = low digit of product arg1 * arg2 - resh = high digit of product arg1 * arg2. */ -#if SIZEOF_PTR == 4 && defined(ARCH_UINT64_TYPE) -#define BngMult(resh,resl,arg1,arg2) { \ - ARCH_UINT64_TYPE p = (ARCH_UINT64_TYPE)(arg1) * (ARCH_UINT64_TYPE)(arg2); \ - resh = p >> 32; \ - resl = p; \ -} -#else -#define BngMult(resh,resl,arg1,arg2) { \ - bngdigit p11 = BngLowHalf(arg1) * BngLowHalf(arg2); \ - bngdigit p12 = BngLowHalf(arg1) * BngHighHalf(arg2); \ - bngdigit p21 = BngHighHalf(arg1) * BngLowHalf(arg2); \ - bngdigit p22 = BngHighHalf(arg1) * BngHighHalf(arg2); \ - resh = p22 + (p12 >> BNG_BITS_PER_HALF_DIGIT) \ - + (p21 >> BNG_BITS_PER_HALF_DIGIT); \ - BngAdd3(resl, resh, \ - p11, p12 << BNG_BITS_PER_HALF_DIGIT, p21 << BNG_BITS_PER_HALF_DIGIT); \ -} -#endif -#endif - -#ifndef BngDiv -/* Divide the double-width number nh:nl by d. - Require d != 0 and nh < d. - Store quotient in quo, remainder in rem. - Can be slow if d is not normalized. */ -#define BngDiv(quo,rem,nh,nl,d) bng_div_aux(&(quo),&(rem),nh,nl,d) -#define BngDivNeedsNormalization - -static void bng_div_aux(bngdigit * quo, bngdigit * rem, - bngdigit nh, bngdigit nl, bngdigit d) -{ - bngdigit dl, dh, ql, qh, pl, ph, nsaved; - - dl = BngLowHalf(d); - dh = BngHighHalf(d); - /* Under-estimate the top half of the quotient (qh) */ - qh = nh / (dh + 1); - /* Shift nh:nl right by BNG_BITS_PER_HALF_DIGIT bits, - so that we focus on the top 1.5 digits of the numerator. - Then, subtract (qh * d) from nh:nl. */ - nsaved = BngLowHalf(nl); - ph = qh * dh; - pl = qh * dl; - nh -= ph; /* Subtract before shifting so that carry propagates for free */ - nl = (nl >> BNG_BITS_PER_HALF_DIGIT) | (nh << BNG_BITS_PER_HALF_DIGIT); - nh = (nh >> BNG_BITS_PER_HALF_DIGIT); - nh -= (nl < pl); /* Borrow */ - nl -= pl; - /* Adjust estimate qh until nh:nl < 0:d */ - while (nh != 0 || nl >= d) { - nh -= (nl < d); /* Borrow */ - nl -= d; - qh++; - } - /* Under-estimate the bottom half of the quotient (ql) */ - ql = nl / (dh + 1); - /* Shift nh:nl left by BNG_BITS_PER_HALF_DIGIT bits, restoring the - low bits we saved earlier, so that we focus on the bottom 1.5 digit - of the numerator. Then, subtract (ql * d) from nh:nl. */ - ph = ql * dh; - pl = ql * dl; - nl -= ph; /* Subtract before shifting so that carry propagates for free */ - nh = (nl >> BNG_BITS_PER_HALF_DIGIT); - nl = (nl << BNG_BITS_PER_HALF_DIGIT) | nsaved; - nh -= (nl < pl); /* Borrow */ - nl -= pl; - /* Adjust estimate ql until nh:nl < 0:d */ - while (nh != 0 || nl >= d) { - nh -= (nl < d); /* Borrow */ - nl -= d; - ql++; - } - /* We're done */ - *quo = (qh << BNG_BITS_PER_HALF_DIGIT) | ql; - *rem = nl; -} - -#endif diff --git a/otherlibs/num/bng_ia32.c b/otherlibs/num/bng_ia32.c deleted file mode 100644 index 6b6cabd2c7..0000000000 --- a/otherlibs/num/bng_ia32.c +++ /dev/null @@ -1,411 +0,0 @@ -/**************************************************************************/ -/* */ -/* OCaml */ -/* */ -/* Xavier Leroy, projet Cristal, INRIA Rocquencourt */ -/* */ -/* Copyright 2003 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. */ -/* */ -/**************************************************************************/ - -/* Code specific to the Intel IA32 (x86) architecture. */ - -#define BngAdd2(res,carryout,arg1,arg2) \ - asm("xorl %1, %1 \n\t" \ - "addl %3, %0 \n\t" \ - "setc %b1" \ - : "=r" (res), "=&q" (carryout) \ - : "0" (arg1), "rm" (arg2)) - -#define BngSub2(res,carryout,arg1,arg2) \ - asm("xorl %1, %1 \n\t" \ - "subl %3, %0 \n\t" \ - "setc %b1" \ - : "=r" (res), "=&q" (carryout) \ - : "0" (arg1), "rm" (arg2)) - -#define BngMult(resh,resl,arg1,arg2) \ - asm("mull %3" \ - : "=a" (resl), "=d" (resh) \ - : "a" (arg1), "r" (arg2)) - -#define BngDiv(quo,rem,nh,nl,d) \ - asm("divl %4" \ - : "=a" (quo), "=d" (rem) \ - : "a" (nl), "d" (nh), "r" (d)) - -/* Reimplementation in asm of some of the bng operations. */ - -static bngcarry bng_ia32_add - (bng a/*[alen]*/, bngsize alen, - bng b/*[blen]*/, bngsize blen, - bngcarry carry) -{ - bngdigit tmp; - alen -= blen; - if (blen > 0) { - asm("negb %b3 \n\t" - "1: \n\t" - "movl (%0), %4 \n\t" - "adcl (%1), %4 \n\t" - "movl %4, (%0) \n\t" - "leal 4(%0), %0 \n\t" - "leal 4(%1), %1 \n\t" - "decl %2 \n\t" - "jnz 1b \n\t" - "setc %b3" - : "+&r" (a), "+&r" (b), "+&r" (blen), "+&q" (carry), "=&r" (tmp)); - } - if (carry == 0 || alen == 0) return carry; - do { - if (++(*a) != 0) return 0; - a++; - } while (--alen); - return 1; -} - -static bngcarry bng_ia32_sub - (bng a/*[alen]*/, bngsize alen, - bng b/*[blen]*/, bngsize blen, - bngcarry carry) -{ - bngdigit tmp; - alen -= blen; - if (blen > 0) { - asm("negb %b3 \n\t" - "1: \n\t" - "movl (%0), %4 \n\t" - "sbbl (%1), %4 \n\t" - "movl %4, (%0) \n\t" - "leal 4(%0), %0 \n\t" - "leal 4(%1), %1 \n\t" - "decl %2 \n\t" - "jnz 1b \n\t" - "setc %b3" - : "+&r" (a), "+&r" (b), "+&r" (blen), "+&q" (carry), "=&r" (tmp)); - } - if (carry == 0 || alen == 0) return carry; - do { - if ((*a)-- != 0) return 0; - a++; - } while (--alen); - return 1; -} - -static bngdigit bng_ia32_mult_add_digit - (bng a/*[alen]*/, bngsize alen, - bng b/*[blen]*/, bngsize blen, - bngdigit d) -{ - bngdigit out; - bngcarry carry; - - alen -= blen; - out = 0; - if (blen > 0) { - asm("1: \n\t" - "movl (%1), %%eax \n\t" - "mull %4\n\t" /* edx:eax = d * next digit of b */ - "addl (%0), %%eax \n\t" /* add next digit of a to eax */ - "adcl $0, %%edx \n\t" /* accumulate carry in edx */ - "addl %3, %%eax \n\t" /* add out to eax */ - "adcl $0, %%edx \n\t" /* accumulate carry in edx */ - "movl %%eax, (%0) \n\t" /* eax is next digit of result */ - "movl %%edx, %3 \n\t" /* edx is next out */ - "leal 4(%0), %0 \n\t" - "leal 4(%1), %1 \n\t" - "decl %2 \n\t" - "jnz 1b" - : "+&r" (a), "+&r" (b), "+&r" (blen), "=m" (out) - : "m" (d) - : "eax", "edx"); - } - if (alen == 0) return out; - /* current digit of a += out */ - BngAdd2(*a, carry, *a, out); - a++; - alen--; - /* Propagate carry */ - if (carry == 0 || alen == 0) return carry; - do { - if (++(*a) != 0) return 0; - a++; - } while (--alen); - return 1; -} - -static bngdigit bng_ia32_mult_sub_digit - (bng a/*[alen]*/, bngsize alen, - bng b/*[blen]*/, bngsize blen, - bngdigit d) -{ - bngdigit out, tmp; - bngcarry carry; - - alen -= blen; - out = 0; - if (blen > 0) { - asm("1: \n\t" - "movl (%1), %%eax \n\t" - "movl (%0), %4 \n\t" - "mull %5\n\t" /* edx:eax = d * next digit of b */ - "subl %%eax, %4 \n\t" /* subtract eax from next digit of a */ - "adcl $0, %%edx \n\t" /* accumulate carry in edx */ - "subl %3, %4 \n\t" /* subtract out */ - "adcl $0, %%edx \n\t" /* accumulate carry in edx */ - "movl %4, (%0) \n\t" /* store next digit of result */ - "movl %%edx, %3 \n\t" /* edx is next out */ - "leal 4(%0), %0 \n\t" - "leal 4(%1), %1 \n\t" - "decl %2 \n\t" - "jnz 1b" - : "+&r" (a), "+&r" (b), "=m" (blen), "=m" (out), "=&r" (tmp) - : "m" (d) - : "eax", "edx"); - } - if (alen == 0) return out; - /* current digit of a -= out */ - BngSub2(*a, carry, *a, out); - a++; - alen--; - /* Propagate carry */ - if (carry == 0 || alen == 0) return carry; - do { - if ((*a)-- != 0) return 0; - a++; - } while (--alen); - return 1; -} - -/* This is another asm implementation of some of the bng operations, - using SSE2 operations to provide 64-bit arithmetic. - This is faster than the plain IA32 code above on the Pentium 4. - (Arithmetic operations with carry are slow on the Pentium 4). */ - -#if BNG_ASM_LEVEL >= 2 - -static bngcarry bng_ia32sse2_add - (bng a/*[alen]*/, bngsize alen, - bng b/*[blen]*/, bngsize blen, - bngcarry carry) -{ - alen -= blen; - if (blen > 0) { - asm("movd %3, %%mm0 \n\t" /* MM0 is carry */ - "1: \n\t" - "movd (%0), %%mm1 \n\t" /* MM1 is next digit of a */ - "movd (%1), %%mm2 \n\t" /* MM2 is next digit of b */ - "paddq %%mm1, %%mm0 \n\t" /* Add carry (64 bits) */ - "paddq %%mm2, %%mm0 \n\t" /* Add digits (64 bits) */ - "movd %%mm0, (%0) \n\t" /* Store low 32 bits of result */ - "psrlq $32, %%mm0 \n\t" /* Next carry is top 32 bits of results */ - "addl $4, %0\n\t" - "addl $4, %1\n\t" - "subl $1, %2\n\t" - "jne 1b \n\t" - "movd %%mm0, %3 \n\t" - "emms" - : "+&r" (a), "+&r" (b), "+&r" (blen), "+&rm" (carry)); - } - if (carry == 0 || alen == 0) return carry; - do { - if (++(*a) != 0) return 0; - a++; - } while (--alen); - return 1; -} - -static bngcarry bng_ia32sse2_sub - (bng a/*[alen]*/, bngsize alen, - bng b/*[blen]*/, bngsize blen, - bngcarry carry) -{ - alen -= blen; - if (blen > 0) { - asm("movd %3, %%mm0 \n\t" /* MM0 is carry */ - "1: \n\t" - "movd (%0), %%mm1 \n\t" /* MM1 is next digit of a */ - "movd (%1), %%mm2 \n\t" /* MM2 is next digit of b */ - "psubq %%mm0, %%mm1 \n\t" /* Subtract carry (64 bits) */ - "psubq %%mm2, %%mm1 \n\t" /* Subtract digits (64 bits) */ - "movd %%mm1, (%0) \n\t" /* Store low 32 bits of result */ - "psrlq $63, %%mm1 \n\t" /* Next carry is sign bit of result */ - "movq %%mm1, %%mm0 \n\t" - "addl $4, %0\n\t" - "addl $4, %1\n\t" - "subl $1, %2\n\t" - "jne 1b \n\t" - "movd %%mm0, %3 \n\t" - "emms" - : "+&r" (a), "+&r" (b), "+&r" (blen), "+&rm" (carry)); - } - if (carry == 0 || alen == 0) return carry; - do { - if ((*a)-- != 0) return 0; - a++; - } while (--alen); - return 1; -} - -static bngdigit bng_ia32sse2_mult_add_digit - (bng a/*[alen]*/, bngsize alen, - bng b/*[blen]*/, bngsize blen, - bngdigit d) -{ - bngdigit out; - bngcarry carry; - - alen -= blen; - out = 0; - if (blen > 0) { - asm("pxor %%mm0, %%mm0 \n\t" /* MM0 is carry */ - "movd %4, %%mm7 \n\t" /* MM7 is digit d */ - "1: \n\t" - "movd (%0), %%mm1 \n\t" /* MM1 is next digit of a */ - "movd (%1), %%mm2 \n\t" /* MM2 is next digit of b */ - "pmuludq %%mm7, %%mm2 \n\t" /* MM2 = d * digit of b */ - "paddq %%mm1, %%mm0 \n\t" /* Add product and carry ... */ - "paddq %%mm2, %%mm0 \n\t" /* ... and digit of a */ - "movd %%mm0, (%0) \n\t" /* Store low 32 bits of result */ - "psrlq $32, %%mm0 \n\t" /* Next carry is high 32 bits result */ - "addl $4, %0\n\t" - "addl $4, %1\n\t" - "subl $1, %2\n\t" - "jne 1b \n\t" - "movd %%mm0, %3 \n\t" - "emms" - : "+&r" (a), "+&r" (b), "+&r" (blen), "=&rm" (out) - : "m" (d)); - } - if (alen == 0) return out; - /* current digit of a += out */ - BngAdd2(*a, carry, *a, out); - a++; - alen--; - /* Propagate carry */ - if (carry == 0 || alen == 0) return carry; - do { - if (++(*a) != 0) return 0; - a++; - } while (--alen); - return 1; -} - -static bngdigit bng_ia32sse2_mult_sub_digit - (bng a/*[alen]*/, bngsize alen, - bng b/*[blen]*/, bngsize blen, - bngdigit d) -{ - static unsigned long long bias1 = 0xFFFFFFFF00000000ULL - 0xFFFFFFFFULL; - static unsigned long bias2 = 0xFFFFFFFFUL; - bngdigit out; - bngcarry carry; - - alen -= blen; - out = 0; - if (blen > 0) { - /* Carry C is represented by ENC(C) = 0xFFFFFFFF - C (one's complement) */ - asm("movd %6, %%mm0 \n\t" /* MM0 is carry (initially 0xFFFFFFFF) */ - "movq %5, %%mm6 \n\t" /* MM6 is magic constant bias1 */ - "movd %4, %%mm7 \n\t" /* MM7 is digit d */ - "1: \n\t" - "movd (%0), %%mm1 \n\t" /* MM1 is next digit of a */ - "movd (%1), %%mm2 \n\t" /* MM2 is next digit of b */ - "paddq %%mm6, %%mm1 \n\t" /* bias digit of a */ - "pmuludq %%mm7, %%mm2 \n\t" /* MM2 = d * digit of b */ - /* Compute - digit of a + ENC(carry) + 0xFFFFFFFF00000000 - 0xFFFFFFFF - product - = digit of a - carry + 0xFFFFFFFF00000000 - product - = digit of a - carry - productlow + (ENC(nextcarry) << 32) */ - "psubq %%mm2, %%mm1 \n\t" - "paddq %%mm1, %%mm0 \n\t" - "movd %%mm0, (%0) \n\t" /* Store low 32 bits of result */ - "psrlq $32, %%mm0 \n\t" /* Next carry is 32 high bits of result */ - "addl $4, %0\n\t" - "addl $4, %1\n\t" - "subl $1, %2\n\t" - "jne 1b \n\t" - "movd %%mm0, %3 \n\t" - "emms" - : "+&r" (a), "+&r" (b), "+&r" (blen), "=&rm" (out) - : "m" (d), "m" (bias1), "m" (bias2)); - out = ~out; /* Undo encoding on out digit */ - } - if (alen == 0) return out; - /* current digit of a -= out */ - BngSub2(*a, carry, *a, out); - a++; - alen--; - /* Propagate carry */ - if (carry == 0 || alen == 0) return carry; - do { - if ((*a)-- != 0) return 0; - a++; - } while (--alen); - return 1; -} - -/* Detect whether SSE2 instructions are supported */ - -static int bng_ia32_sse2_supported(void) -{ - unsigned int flags, newflags, max_id, capabilities; - -#define EFLAG_CPUID 0x00200000 -#define CPUID_IDENTIFY 0 -#define CPUID_CAPABILITIES 1 -#define SSE2_CAPABILITY 26 - - /* Check if processor has CPUID instruction */ - asm("pushfl \n\t" - "popl %0" - : "=r" (flags) : ); - newflags = flags ^ EFLAG_CPUID; /* CPUID detection flag */ - asm("pushfl \n\t" - "pushl %1 \n\t" - "popfl \n\t" - "pushfl \n\t" - "popl %0 \n\t" - "popfl" - : "=r" (flags) : "r" (newflags)); - /* If CPUID detection flag cannot be changed, CPUID instruction is not - available */ - if ((flags & EFLAG_CPUID) != (newflags & EFLAG_CPUID)) return 0; - /* See if SSE2 extensions are supported */ - asm("pushl %%ebx \n\t" /* need to preserve %ebx for PIC */ - "cpuid \n\t" - "popl %%ebx" - : "=a" (max_id) : "a" (CPUID_IDENTIFY): "ecx", "edx"); - if (max_id < 1) return 0; - asm("pushl %%ebx \n\t" - "cpuid \n\t" - "popl %%ebx" - : "=d" (capabilities) : "a" (CPUID_CAPABILITIES) : "ecx"); - return capabilities & (1 << SSE2_CAPABILITY); -} - -#endif - -static void bng_ia32_setup_ops(void) -{ -#if BNG_ASM_LEVEL >= 2 - if (bng_ia32_sse2_supported()) { - bng_ops.add = bng_ia32sse2_add; - bng_ops.sub = bng_ia32sse2_sub; - bng_ops.mult_add_digit = bng_ia32sse2_mult_add_digit; - bng_ops.mult_sub_digit = bng_ia32sse2_mult_sub_digit; - return; - } -#endif - bng_ops.add = bng_ia32_add; - bng_ops.sub = bng_ia32_sub; - bng_ops.mult_add_digit = bng_ia32_mult_add_digit; - bng_ops.mult_sub_digit = bng_ia32_mult_sub_digit; -} - -#define BNG_SETUP_OPS bng_ia32_setup_ops() diff --git a/otherlibs/num/bng_ppc.c b/otherlibs/num/bng_ppc.c deleted file mode 100644 index f4c098cf08..0000000000 --- a/otherlibs/num/bng_ppc.c +++ /dev/null @@ -1,94 +0,0 @@ -/**************************************************************************/ -/* */ -/* OCaml */ -/* */ -/* Xavier Leroy, projet Cristal, INRIA Rocquencourt */ -/* */ -/* Copyright 2003 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. */ -/* */ -/**************************************************************************/ - -/* Code specific to the PowerPC architecture. */ - -#define BngAdd2(res,carryout,arg1,arg2) \ - asm("addc %0, %2, %3 \n\t" \ - "li %1, 0 \n\t" \ - "addze %1, %1" \ - : "=r" (res), "=r" (carryout) \ - : "r" (arg1), "r" (arg2)) - -#define BngAdd2Carry(res,carryout,arg1,arg2,carryin) \ - asm("addic %1, %4, -1 \n\t" \ - "adde %0, %2, %3 \n\t" \ - "li %1, 0 \n\t" \ - "addze %1, %1" \ - : "=r" (res), "=&r" (carryout) \ - : "r" (arg1), "r" (arg2), "1" (carryin)) - -#define BngAdd3(res,carryaccu,arg1,arg2,arg3) \ - asm("addc %0, %2, %3 \n\t" \ - "addze %1, %1 \n\t" \ - "addc %0, %0, %4 \n\t" \ - "addze %1, %1" \ - : "=&r" (res), "=&r" (carryaccu) \ - : "r" (arg1), "r" (arg2), "r" (arg3), "1" (carryaccu)) - -/* The "subtract" instructions interpret carry differently than what we - need: the processor carry bit CA is 1 if no carry occured, - 0 if a carry occured. In other terms, CA = !carry. - Thus, subfe rd,ra,rb computes rd = ra - rb - !CA - subfe rd,rd,rd sets rd = - !CA - subfe rd,rd,rd; neg rd, rd sets rd = !CA and recovers "our" carry. */ - -#define BngSub2(res,carryout,arg1,arg2) \ - asm("subfc %0, %3, %2 \n\t" \ - "subfe %1, %1, %1\n\t" \ - "neg %1, %1" \ - : "=r" (res), "=r" (carryout) \ - : "r" (arg1), "r" (arg2)) - -#define BngSub2Carry(res,carryout,arg1,arg2,carryin) \ - asm("subfic %1, %4, 0 \n\t" \ - "subfe %0, %3, %2 \n\t" \ - "subfe %1, %1, %1 \n\t" \ - "neg %1, %1" \ - : "=r" (res), "=&r" (carryout) \ - : "r" (arg1), "r" (arg2), "1" (carryin)) - -/* Here is what happens with carryaccu: - neg %1, %1 carryaccu = -carryaccu - addze %1, %1 carryaccu += !carry1 - addze %1, %1 carryaccu += !carry2 - subifc %1, %1, 2 carryaccu = 2 - carryaccu - Thus, carryaccu_final = carryaccu_initial + 2 - (1 - carry1) - (1 - carry2) - = carryaccu_initial + carry1 + carry2 -*/ - -#define BngSub3(res,carryaccu,arg1,arg2,arg3) \ - asm("neg %1, %1 \n\t" \ - "subfc %0, %3, %2 \n\t" \ - "addze %1, %1 \n\t" \ - "subfc %0, %4, %0 \n\t" \ - "addze %1, %1 \n\t" \ - "subfic %1, %1, 2 \n\t" \ - : "=&r" (res), "=&r" (carryaccu) \ - : "r" (arg1), "r" (arg2), "r" (arg3), "1" (carryaccu)) - -#if defined(__ppc64__) || defined(__PPC64__) -#define BngMult(resh,resl,arg1,arg2) \ - asm("mulld %0, %2, %3 \n\t" \ - "mulhdu %1, %2, %3" \ - : "=&r" (resl), "=r" (resh) \ - : "r" (arg1), "r" (arg2)) -#else -#define BngMult(resh,resl,arg1,arg2) \ - asm("mullw %0, %2, %3 \n\t" \ - "mulhwu %1, %2, %3" \ - : "=&r" (resl), "=r" (resh) \ - : "r" (arg1), "r" (arg2)) -#endif diff --git a/otherlibs/num/bng_sparc.c b/otherlibs/num/bng_sparc.c deleted file mode 100644 index c007cb7742..0000000000 --- a/otherlibs/num/bng_sparc.c +++ /dev/null @@ -1,77 +0,0 @@ -/**************************************************************************/ -/* */ -/* OCaml */ -/* */ -/* Xavier Leroy, projet Cristal, INRIA Rocquencourt */ -/* */ -/* Copyright 2003 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. */ -/* */ -/**************************************************************************/ - -/* Code specific to the SPARC (V8 and above) architecture. */ - -#define BngAdd2(res,carryout,arg1,arg2) \ - asm("addcc %2, %3, %0 \n\t" \ - "addx %%g0, 0, %1" \ - : "=r" (res), "=r" (carryout) \ - : "r" (arg1), "r" (arg2) \ - : "cc") - -#define BngAdd2Carry(res,carryout,arg1,arg2,carryin) \ - asm("subcc %%g0, %4, %%g0 \n\t" \ - "addxcc %2, %3, %0 \n\t" \ - "addx %%g0, 0, %1" \ - : "=r" (res), "=r" (carryout) \ - : "r" (arg1), "r" (arg2), "r" (carryin) \ - : "cc") - -#define BngAdd3(res,carryaccu,arg1,arg2,arg3) \ - asm("addcc %2, %3, %0 \n\t" \ - "addx %1, 0, %1 \n\t" \ - "addcc %0, %4, %0 \n\t" \ - "addx %1, 0, %1" \ - : "=r" (res), "=r" (carryaccu) \ - : "r" (arg1), "r" (arg2), "r" (arg3), "1" (carryaccu) \ - : "cc") - -#define BngSub2(res,carryout,arg1,arg2) \ - asm("subcc %2, %3, %0 \n\t" \ - "addx %%g0, 0, %1" \ - : "=r" (res), "=r" (carryout) \ - : "r" (arg1), "r" (arg2) \ - : "cc") - -#define BngSub2Carry(res,carryout,arg1,arg2,carryin) \ - asm("subcc %%g0, %4, %%g0 \n\t" \ - "subxcc %2, %3, %0 \n\t" \ - "addx %%g0, 0, %1" \ - : "=r" (res), "=r" (carryout) \ - : "r" (arg1), "r" (arg2), "r" (carryin) \ - : "cc") - -#define BngSub3(res,carryaccu,arg1,arg2,arg3) \ - asm("subcc %2, %3, %0 \n\t" \ - "addx %1, 0, %1 \n\t" \ - "subcc %0, %4, %0 \n\t" \ - "addx %1, 0, %1" \ - : "=r" (res), "=r" (carryaccu) \ - : "r" (arg1), "r" (arg2), "r" (arg3), "1" (carryaccu) \ - : "cc") - -#define BngMult(resh,resl,arg1,arg2) \ - asm("umul %2, %3, %0 \n\t" \ - "rd %%y, %1" \ - : "=r" (resl), "=r" (resh) \ - : "r" (arg1), "r" (arg2)) - -#define BngDiv(quo,rem,nh,nl,d) \ - asm("wr %1, %%y \n\t" \ - "udiv %2, %3, %0" \ - : "=r" (quo) \ - : "r" (nh), "r" (nl), "r" (d)); \ - rem = nl - d * quo diff --git a/otherlibs/num/int_misc.ml b/otherlibs/num/int_misc.ml deleted file mode 100644 index d7d7190ea7..0000000000 --- a/otherlibs/num/int_misc.ml +++ /dev/null @@ -1,36 +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. *) -(* *) -(**************************************************************************) - -(* Some extra operations on integers *) - -let rec gcd_int i1 i2 = - if i2 = 0 then abs i1 else gcd_int i2 (i1 mod i2) -;; - -let rec num_bits_int_aux n = - if n = 0 then 0 else succ(num_bits_int_aux (n lsr 1));; - -let num_bits_int n = num_bits_int_aux (abs n);; - -let sign_int i = if i = 0 then 0 else if i > 0 then 1 else -1;; - -let length_of_int = Sys.word_size - 2;; - -let monster_int = 1 lsl length_of_int;; -let biggest_int = monster_int - 1;; -let least_int = - biggest_int;; - -let compare_int n1 n2 = - if n1 == n2 then 0 else if n1 > n2 then 1 else -1;; diff --git a/otherlibs/num/int_misc.mli b/otherlibs/num/int_misc.mli deleted file mode 100644 index 1ee11ba5f0..0000000000 --- a/otherlibs/num/int_misc.mli +++ /dev/null @@ -1,25 +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. *) -(* *) -(**************************************************************************) - -(* Some extra operations on integers *) - -val gcd_int: int -> int -> int -val num_bits_int: int -> int -val compare_int: int -> int -> int -val sign_int: int -> int -val length_of_int: int -val biggest_int: int -val least_int: int -val monster_int: int diff --git a/otherlibs/num/nat.h b/otherlibs/num/nat.h deleted file mode 100644 index 45e7b9577c..0000000000 --- a/otherlibs/num/nat.h +++ /dev/null @@ -1,18 +0,0 @@ -/**************************************************************************/ -/* */ -/* OCaml */ -/* */ -/* Xavier Leroy, projet Cristal, INRIA Rocquencourt */ -/* */ -/* Copyright 1999 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. */ -/* */ -/**************************************************************************/ - -/* Nats are represented as unstructured blocks with tag Custom_tag. */ - -#define Digit_val(nat,pos) (((bng) Data_custom_val(nat))[pos]) diff --git a/otherlibs/num/nat.ml b/otherlibs/num/nat.ml deleted file mode 100644 index c7a2669840..0000000000 --- a/otherlibs/num/nat.ml +++ /dev/null @@ -1,594 +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 - -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 - -*******) - -(* - 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 digits = "0123456789ABCDEF" - -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) diff --git a/otherlibs/num/nat.mli b/otherlibs/num/nat.mli deleted file mode 100644 index 803a6537d2..0000000000 --- a/otherlibs/num/nat.mli +++ /dev/null @@ -1,89 +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. *) -(* *) -(**************************************************************************) - -(* Module [Nat]: operations on natural numbers *) - -type nat - -(* Natural numbers (type [nat]) are positive integers of arbitrary size. - All operations on [nat] are performed in-place. *) - -external create_nat: int -> nat = "create_nat" -val make_nat: int -> 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" -val copy_nat: nat -> int -> int -> 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" -val length_nat : nat -> int -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" -val is_zero_nat: nat -> int -> int -> bool -val is_nat_int: nat -> int -> int -> bool -val int_of_nat: nat -> int -val nat_of_int: int -> nat -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" -val eq_nat : nat -> int -> int -> nat -> int -> int -> bool -val le_nat : nat -> int -> int -> nat -> int -> int -> bool -val lt_nat : nat -> int -> int -> nat -> int -> int -> bool -val ge_nat : nat -> int -> int -> nat -> int -> int -> bool -val gt_nat : nat -> int -> int -> nat -> int -> int -> bool -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" -val gcd_nat : nat -> int -> int -> nat -> int -> int -> int -val sqrt_nat : nat -> int -> int -> nat -val string_of_nat : nat -> string -val nat_of_string : string -> nat -val sys_nat_of_string : int -> string -> int -> int -> nat -val float_of_nat : nat -> float -val make_power_base : int -> nat -> int * int -val power_base_int : int -> int -> nat -val length_of_digit: int diff --git a/otherlibs/num/nat_stubs.c b/otherlibs/num/nat_stubs.c deleted file mode 100644 index 5a07a801df..0000000000 --- a/otherlibs/num/nat_stubs.c +++ /dev/null @@ -1,421 +0,0 @@ -/**************************************************************************/ -/* */ -/* OCaml */ -/* */ -/* Xavier Leroy, 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. */ -/* */ -/**************************************************************************/ - -#define CAML_INTERNALS - -#include "caml/alloc.h" -#include "caml/config.h" -#include "caml/custom.h" -#include "caml/intext.h" -#include "caml/fail.h" -#include "caml/hash.h" -#include "caml/memory.h" -#include "caml/mlvalues.h" - -#include "bng.h" -#include "nat.h" - -/* Stub code for the Nat module. */ - -static intnat hash_nat(value); -static void serialize_nat(value, uintnat *, uintnat *); -static uintnat deserialize_nat(void * dst); - -static struct custom_operations nat_operations = { - "_nat", - custom_finalize_default, - custom_compare_default, - hash_nat, - serialize_nat, - deserialize_nat, - custom_compare_ext_default -}; - -CAMLprim value initialize_nat(value unit) -{ - bng_init(); - caml_register_custom_operations(&nat_operations); - return Val_unit; -} - -CAMLprim value create_nat(value size) -{ - mlsize_t sz = Long_val(size); - - return caml_alloc_custom(&nat_operations, sz * sizeof(value), 0, 1); -} - -CAMLprim value length_nat(value nat) -{ - return Val_long(Wosize_val(nat) - 1); -} - -CAMLprim value set_to_zero_nat(value nat, value ofs, value len) -{ - bng_zero(&Digit_val(nat, Long_val(ofs)), Long_val(len)); - return Val_unit; -} - -CAMLprim value blit_nat(value nat1, value ofs1, - value nat2, value ofs2, - value len) -{ - bng_assign(&Digit_val(nat1, Long_val(ofs1)), - &Digit_val(nat2, Long_val(ofs2)), - Long_val(len)); - return Val_unit; -} - -CAMLprim value set_digit_nat(value nat, value ofs, value digit) -{ - Digit_val(nat, Long_val(ofs)) = Long_val(digit); - return Val_unit; -} - -CAMLprim value nth_digit_nat(value nat, value ofs) -{ - return Val_long(Digit_val(nat, Long_val(ofs))); -} - -CAMLprim value set_digit_nat_native(value nat, value ofs, value digit) -{ - Digit_val(nat, Long_val(ofs)) = Nativeint_val(digit); - return Val_unit; -} - -CAMLprim value nth_digit_nat_native(value nat, value ofs) -{ - return caml_copy_nativeint(Digit_val(nat, Long_val(ofs))); -} - -CAMLprim value num_digits_nat(value nat, value ofs, value len) -{ - return Val_long(bng_num_digits(&Digit_val(nat, Long_val(ofs)), - Long_val(len))); -} - -CAMLprim value num_leading_zero_bits_in_digit(value nat, value ofs) -{ - return - Val_long(bng_leading_zero_bits(Digit_val(nat, Long_val(ofs)))); -} - -CAMLprim value is_digit_int(value nat, value ofs) -{ - return Val_bool(Digit_val(nat, Long_val(ofs)) <= Max_long); -} - -CAMLprim value is_digit_zero(value nat, value ofs) -{ - return Val_bool(Digit_val(nat, Long_val(ofs)) == 0); -} - -CAMLprim value is_digit_normalized(value nat, value ofs) -{ - return - Val_bool(Digit_val(nat, Long_val(ofs)) - & ((bngdigit)1 << (BNG_BITS_PER_DIGIT-1))); -} - -CAMLprim value is_digit_odd(value nat, value ofs) -{ - return Val_bool(Digit_val(nat, Long_val(ofs)) & 1); -} - -CAMLprim value incr_nat(value nat, value ofs, value len, value carry_in) -{ - return Val_long(bng_add_carry(&Digit_val(nat, Long_val(ofs)), - Long_val(len), Long_val(carry_in))); -} - -value add_nat_native(value nat1, value ofs1, value len1, - value nat2, value ofs2, value len2, value carry_in) -{ - return Val_long(bng_add(&Digit_val(nat1, Long_val(ofs1)), Long_val(len1), - &Digit_val(nat2, Long_val(ofs2)), Long_val(len2), - Long_val(carry_in))); -} - -CAMLprim value add_nat(value *argv, int argn) -{ - return add_nat_native(argv[0], argv[1], argv[2], argv[3], - argv[4], argv[5], argv[6]); -} - -CAMLprim value complement_nat(value nat, value ofs, value len) -{ - bng_complement(&Digit_val(nat, Long_val(ofs)), Long_val(len)); - return Val_unit; -} - -CAMLprim value decr_nat(value nat, value ofs, value len, value carry_in) -{ - return Val_long(1 ^ bng_sub_carry(&Digit_val(nat, Long_val(ofs)), - Long_val(len), 1 ^ Long_val(carry_in))); -} - -value sub_nat_native(value nat1, value ofs1, value len1, - value nat2, value ofs2, value len2, value carry_in) -{ - return Val_long(1 ^ bng_sub(&Digit_val(nat1, Long_val(ofs1)), Long_val(len1), - &Digit_val(nat2, Long_val(ofs2)), Long_val(len2), - 1 ^ Long_val(carry_in))); -} - -CAMLprim value sub_nat(value *argv, int argn) -{ - return sub_nat_native(argv[0], argv[1], argv[2], argv[3], - argv[4], argv[5], argv[6]); -} - -value mult_digit_nat_native(value nat1, value ofs1, value len1, - value nat2, value ofs2, value len2, - value nat3, value ofs3) -{ - return - Val_long(bng_mult_add_digit( - &Digit_val(nat1, Long_val(ofs1)), Long_val(len1), - &Digit_val(nat2, Long_val(ofs2)), Long_val(len2), - Digit_val(nat3, Long_val(ofs3)))); -} - -CAMLprim value mult_digit_nat(value *argv, int argn) -{ - return mult_digit_nat_native(argv[0], argv[1], argv[2], argv[3], - argv[4], argv[5], argv[6], argv[7]); -} - -value mult_nat_native(value nat1, value ofs1, value len1, - value nat2, value ofs2, value len2, - value nat3, value ofs3, value len3) -{ - return - Val_long(bng_mult_add(&Digit_val(nat1, Long_val(ofs1)), Long_val(len1), - &Digit_val(nat2, Long_val(ofs2)), Long_val(len2), - &Digit_val(nat3, Long_val(ofs3)), Long_val(len3))); -} - -CAMLprim value mult_nat(value *argv, int argn) -{ - return mult_nat_native(argv[0], argv[1], argv[2], argv[3], - argv[4], argv[5], argv[6], argv[7], argv[8]); -} - -value square_nat_native(value nat1, value ofs1, value len1, - value nat2, value ofs2, value len2) -{ - return - Val_long(bng_square_add(&Digit_val(nat1, Long_val(ofs1)), Long_val(len1), - &Digit_val(nat2, Long_val(ofs2)), Long_val(len2))); -} - -CAMLprim value square_nat(value *argv, int argn) -{ - return square_nat_native(argv[0], argv[1], argv[2], - argv[3], argv[4], argv[5]); -} - -value shift_left_nat_native(value nat1, value ofs1, value len1, - value nat2, value ofs2, value nbits) -{ - Digit_val(nat2, Long_val(ofs2)) = - bng_shift_left(&Digit_val(nat1, Long_val(ofs1)), Long_val(len1), - Long_val(nbits)); - return Val_unit; -} - -CAMLprim value shift_left_nat(value *argv, int argn) -{ - return shift_left_nat_native(argv[0], argv[1], argv[2], - argv[3], argv[4], argv[5]); -} - -value div_digit_nat_native(value natq, value ofsq, - value natr, value ofsr, - value nat1, value ofs1, value len1, - value nat2, value ofs2) -{ - Digit_val(natr, Long_val(ofsr)) = - bng_div_rem_digit(&Digit_val(natq, Long_val(ofsq)), - &Digit_val(nat1, Long_val(ofs1)), Long_val(len1), - Digit_val(nat2, Long_val(ofs2))); - return Val_unit; -} - -CAMLprim value div_digit_nat(value *argv, int argn) -{ - return div_digit_nat_native(argv[0], argv[1], argv[2], argv[3], - argv[4], argv[5], argv[6], argv[7], argv[8]); -} - -value div_nat_native(value nat1, value ofs1, value len1, - value nat2, value ofs2, value len2) -{ - bng_div_rem(&Digit_val(nat1, Long_val(ofs1)), Long_val(len1), - &Digit_val(nat2, Long_val(ofs2)), Long_val(len2)); - return Val_unit; -} - -CAMLprim value div_nat(value *argv, int argn) -{ - return div_nat_native(argv[0], argv[1], argv[2], - argv[3], argv[4], argv[5]); -} - -value shift_right_nat_native(value nat1, value ofs1, value len1, - value nat2, value ofs2, value nbits) -{ - Digit_val(nat2, Long_val(ofs2)) = - bng_shift_right(&Digit_val(nat1, Long_val(ofs1)), Long_val(len1), - Long_val(nbits)); - return Val_unit; -} - -CAMLprim value shift_right_nat(value *argv, int argn) -{ - return shift_right_nat_native(argv[0], argv[1], argv[2], - argv[3], argv[4], argv[5]); -} - -CAMLprim value compare_digits_nat(value nat1, value ofs1, - value nat2, value ofs2) -{ - bngdigit d1 = Digit_val(nat1, Long_val(ofs1)); - bngdigit d2 = Digit_val(nat2, Long_val(ofs2)); - if (d1 > d2) return Val_int(1); - if (d1 < d2) return Val_int(-1); - return Val_int(0); -} - -value compare_nat_native(value nat1, value ofs1, value len1, - value nat2, value ofs2, value len2) -{ - return - Val_int(bng_compare(&Digit_val(nat1, Long_val(ofs1)), Long_val(len1), - &Digit_val(nat2, Long_val(ofs2)), Long_val(len2))); -} - -CAMLprim value compare_nat(value *argv, int argn) -{ - return compare_nat_native(argv[0], argv[1], argv[2], - argv[3], argv[4], argv[5]); -} - -CAMLprim value land_digit_nat(value nat1, value ofs1, value nat2, value ofs2) -{ - Digit_val(nat1, Long_val(ofs1)) &= Digit_val(nat2, Long_val(ofs2)); - return Val_unit; -} - -CAMLprim value lor_digit_nat(value nat1, value ofs1, value nat2, value ofs2) -{ - Digit_val(nat1, Long_val(ofs1)) |= Digit_val(nat2, Long_val(ofs2)); - return Val_unit; -} - -CAMLprim value lxor_digit_nat(value nat1, value ofs1, value nat2, value ofs2) -{ - Digit_val(nat1, Long_val(ofs1)) ^= Digit_val(nat2, Long_val(ofs2)); - return Val_unit; -} - -/* The wire format for a nat is: - - 32-bit word: number of 32-bit words in nat - - N 32-bit words (big-endian format) - For little-endian platforms, the memory layout between 32-bit and 64-bit - machines is identical, so we can write the nat using caml_serialize_block_4. - For big-endian 64-bit platforms, we need to swap the two 32-bit halves - of 64-bit words to obtain the correct behavior. */ - -static void serialize_nat(value nat, - uintnat * wsize_32, - uintnat * wsize_64) -{ - mlsize_t len = Wosize_val(nat) - 1; - -#ifdef ARCH_SIXTYFOUR - len = len * 2; /* two 32-bit words per 64-bit digit */ - if (len >= ((mlsize_t)1 << 32)) - caml_failwith("output_value: nat too big"); -#endif - caml_serialize_int_4((int32_t) len); -#if defined(ARCH_SIXTYFOUR) && defined(ARCH_BIG_ENDIAN) - { int32_t * p; - mlsize_t i; - for (i = len, p = Data_custom_val(nat); i > 0; i -= 2, p += 2) { - caml_serialize_int_4(p[1]); /* low 32 bits of 64-bit digit */ - caml_serialize_int_4(p[0]); /* high 32 bits of 64-bit digit */ - } - } -#else - caml_serialize_block_4(Data_custom_val(nat), len); -#endif - *wsize_32 = len * 4; - *wsize_64 = len * 4; -} - -static uintnat deserialize_nat(void * dst) -{ - mlsize_t len; - - len = caml_deserialize_uint_4(); -#if defined(ARCH_SIXTYFOUR) && defined(ARCH_BIG_ENDIAN) - { uint32_t * p; - mlsize_t i; - for (i = len, p = dst; i > 1; i -= 2, p += 2) { - p[1] = caml_deserialize_uint_4(); /* low 32 bits of 64-bit digit */ - p[0] = caml_deserialize_uint_4(); /* high 32 bits of 64-bit digit */ - } - if (i > 0){ - p[1] = caml_deserialize_uint_4(); /* low 32 bits of 64-bit digit */ - p[0] = 0; /* high 32 bits of 64-bit digit */ - ++ len; - } - } -#else - caml_deserialize_block_4(dst, len); -#if defined(ARCH_SIXTYFOUR) - if (len & 1){ - ((uint32_t *) dst)[len] = 0; - ++ len; - } -#endif -#endif - return len * 4; -} - -static intnat hash_nat(value v) -{ - bngsize len, i; - uint32_t h; - - len = bng_num_digits(&Digit_val(v,0), Wosize_val(v) - 1); - h = 0; - for (i = 0; i < len; i++) { - bngdigit d = Digit_val(v, i); -#ifdef ARCH_SIXTYFOUR - /* Mix the two 32-bit halves as if we were on a 32-bit platform, - namely low 32 bits first, then high 32 bits. - Also, ignore final 32 bits if they are zero. */ - h = caml_hash_mix_uint32(h, (uint32_t) d); - d = d >> 32; - if (d == 0 && i + 1 == len) break; - h = caml_hash_mix_uint32(h, (uint32_t) d); -#else - h = caml_hash_mix_uint32(h, d); -#endif - } - return h; -} diff --git a/otherlibs/num/num.ml b/otherlibs/num/num.ml deleted file mode 100644 index 46b70a137f..0000000000 --- a/otherlibs/num/num.ml +++ /dev/null @@ -1,450 +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 -open Big_int -open Arith_flags -open Ratio - -type num = Int of int | Big_int of big_int | Ratio of ratio - (* The type of numbers. *) - -let biggest_INT = big_int_of_int biggest_int -and least_INT = big_int_of_int least_int - -(* Coercion big_int -> num *) -let num_of_big_int bi = - if le_big_int bi biggest_INT && ge_big_int bi least_INT - then Int (int_of_big_int bi) - else Big_int bi - -let normalize_num = function - Int i -> Int i -| Big_int bi -> num_of_big_int bi -| Ratio r -> if is_integer_ratio r - then num_of_big_int (numerator_ratio r) - else Ratio r - -let cautious_normalize_num_when_printing n = - if (!normalize_ratio_when_printing_flag) then (normalize_num n) else n - -let num_of_ratio r = - ignore (normalize_ratio r); - if not (is_integer_ratio r) then Ratio r - else if is_int_big_int (numerator_ratio r) then - Int (int_of_big_int (numerator_ratio r)) - else Big_int (numerator_ratio r) - -(* Operations on num *) - -let add_num a b = match (a,b) with - ((Int int1), (Int int2)) -> - let r = int1 + int2 in - if (int1 lxor int2) lor (int1 lxor (r lxor (-1))) < 0 - then Int r (* No overflow *) - else Big_int(add_big_int (big_int_of_int int1) (big_int_of_int int2)) - | ((Int i), (Big_int bi)) -> - num_of_big_int (add_int_big_int i bi) - | ((Big_int bi), (Int i)) -> - num_of_big_int (add_int_big_int i bi) - - | ((Int i), (Ratio r)) -> - Ratio (add_int_ratio i r) - | ((Ratio r), (Int i)) -> - Ratio (add_int_ratio i r) - - | ((Big_int bi1), (Big_int bi2)) -> num_of_big_int (add_big_int bi1 bi2) - - | ((Big_int bi), (Ratio r)) -> - Ratio (add_big_int_ratio bi r) - | ((Ratio r), (Big_int bi)) -> - Ratio (add_big_int_ratio bi r) - - | ((Ratio r1), (Ratio r2)) -> num_of_ratio (add_ratio r1 r2) - -let ( +/ ) = add_num - -let minus_num = function - Int i -> if i = monster_int - then Big_int (minus_big_int (big_int_of_int i)) - else Int (-i) -| Big_int bi -> Big_int (minus_big_int bi) -| Ratio r -> Ratio (minus_ratio r) - -let sub_num n1 n2 = add_num n1 (minus_num n2) - -let ( -/ ) = sub_num - -let mult_num a b = match (a,b) with - ((Int int1), (Int int2)) -> - if num_bits_int int1 + num_bits_int int2 < length_of_int - then Int (int1 * int2) - else num_of_big_int (mult_big_int (big_int_of_int int1) - (big_int_of_int int2)) - - | ((Int i), (Big_int bi)) -> - num_of_big_int (mult_int_big_int i bi) - | ((Big_int bi), (Int i)) -> - num_of_big_int (mult_int_big_int i bi) - - | ((Int i), (Ratio r)) -> - num_of_ratio (mult_int_ratio i r) - | ((Ratio r), (Int i)) -> - num_of_ratio (mult_int_ratio i r) - - | ((Big_int bi1), (Big_int bi2)) -> - num_of_big_int (mult_big_int bi1 bi2) - - | ((Big_int bi), (Ratio r)) -> - num_of_ratio (mult_big_int_ratio bi r) - | ((Ratio r), (Big_int bi)) -> - num_of_ratio (mult_big_int_ratio bi r) - - | ((Ratio r1), (Ratio r2)) -> - num_of_ratio (mult_ratio r1 r2) - -let ( */ ) = mult_num - -let square_num = function - Int i -> if 2 * num_bits_int i < length_of_int - then Int (i * i) - else num_of_big_int (square_big_int (big_int_of_int i)) - | Big_int bi -> Big_int (square_big_int bi) - | Ratio r -> Ratio (square_ratio r) - -let div_num n1 n2 = - match n1 with - | Int i1 -> - begin match n2 with - | Int i2 -> - num_of_ratio (create_ratio (big_int_of_int i1) (big_int_of_int i2)) - | Big_int bi2 -> num_of_ratio (create_ratio (big_int_of_int i1) bi2) - | Ratio r2 -> num_of_ratio (div_int_ratio i1 r2) end - - | Big_int bi1 -> - begin match n2 with - | Int i2 -> num_of_ratio (create_ratio bi1 (big_int_of_int i2)) - | Big_int bi2 -> num_of_ratio (create_ratio bi1 bi2) - | Ratio r2 -> num_of_ratio (div_big_int_ratio bi1 r2) end - - | Ratio r1 -> - begin match n2 with - | Int i2 -> num_of_ratio (div_ratio_int r1 i2) - | Big_int bi2 -> num_of_ratio (div_ratio_big_int r1 bi2) - | Ratio r2 -> num_of_ratio (div_ratio r1 r2) end -;; - -let ( // ) = div_num - -let floor_num = function - Int _ as n -> n -| Big_int _ as n -> n -| Ratio r -> num_of_big_int (floor_ratio r) - -(* Coercion with ratio type *) -let ratio_of_num = function - Int i -> ratio_of_int i -| Big_int bi -> ratio_of_big_int bi -| Ratio r -> r -;; - -(* Euclidean division and remainder. The specification is: - - a = b * quo_num a b + mod_num a b - quo_num a b is an integer (Z) - 0 <= mod_num a b < |b| - -A correct but slow implementation is: - - quo_num a b = - if b >= 0 then floor_num (div_num a b) - else minus_num (floor_num (div_num a (minus_num b))) - - mod_num a b = - sub_num a (mult_num b (quo_num a b)) - - However, this definition is vastly inefficient (cf PR #3473): - we define here a better way of computing the same thing. - - PR#6753: the previous implementation was based on - quo_num a b = floor_num (div_num a b) - which is incorrect for negative b. -*) - -let quo_num n1 n2 = - match n1, n2 with - | Int i1, Int i2 -> - let q = i1 / i2 and r = i1 mod i2 in - Int (if r >= 0 then q else if i2 > 0 then q - 1 else q + 1) - | Int i1, Big_int bi2 -> - num_of_big_int (div_big_int (big_int_of_int i1) bi2) - | Int i1, Ratio r2 -> - num_of_big_int (report_sign_ratio r2 - (floor_ratio (div_int_ratio i1 (abs_ratio r2)))) - | Big_int bi1, Int i2 -> - num_of_big_int (div_big_int bi1 (big_int_of_int i2)) - | Big_int bi1, Big_int bi2 -> - num_of_big_int (div_big_int bi1 bi2) - | Big_int bi1, Ratio r2 -> - num_of_big_int (report_sign_ratio r2 - (floor_ratio (div_big_int_ratio bi1 (abs_ratio r2)))) - | Ratio r1, _ -> - let r2 = ratio_of_num n2 in - num_of_big_int (report_sign_ratio r2 - (floor_ratio (div_ratio r1 (abs_ratio r2)))) - -let mod_num n1 n2 = - match n1, n2 with - | Int i1, Int i2 -> - let r = i1 mod i2 in - Int (if r >= 0 then r else if i2 > 0 then r + i2 else r - i2) - | Int i1, Big_int bi2 -> - num_of_big_int (mod_big_int (big_int_of_int i1) bi2) - | Big_int bi1, Int i2 -> - num_of_big_int (mod_big_int bi1 (big_int_of_int i2)) - | Big_int bi1, Big_int bi2 -> - num_of_big_int (mod_big_int bi1 bi2) - | _, _ -> - sub_num n1 (mult_num n2 (quo_num n1 n2)) - -let power_num_int a b = match (a,b) with - ((Int i), n) -> - (match sign_int n with - 0 -> Int 1 - | 1 -> num_of_big_int (power_int_positive_int i n) - | _ -> Ratio (create_normalized_ratio - unit_big_int (power_int_positive_int i (-n)))) -| ((Big_int bi), n) -> - (match sign_int n with - 0 -> Int 1 - | 1 -> num_of_big_int (power_big_int_positive_int bi n) - | _ -> Ratio (create_normalized_ratio - unit_big_int (power_big_int_positive_int bi (-n)))) -| ((Ratio r), n) -> - (match sign_int n with - 0 -> Int 1 - | 1 -> Ratio (power_ratio_positive_int r n) - | _ -> Ratio (power_ratio_positive_int - (inverse_ratio r) (-n))) - -let power_num_big_int a b = match (a,b) with - ((Int i), n) -> - (match sign_big_int n with - 0 -> Int 1 - | 1 -> num_of_big_int (power_int_positive_big_int i n) - | _ -> Ratio (create_normalized_ratio - unit_big_int - (power_int_positive_big_int i (minus_big_int n)))) -| ((Big_int bi), n) -> - (match sign_big_int n with - 0 -> Int 1 - | 1 -> num_of_big_int (power_big_int_positive_big_int bi n) - | _ -> Ratio (create_normalized_ratio - unit_big_int - (power_big_int_positive_big_int bi (minus_big_int n)))) -| ((Ratio r), n) -> - (match sign_big_int n with - 0 -> Int 1 - | 1 -> Ratio (power_ratio_positive_big_int r n) - | _ -> Ratio (power_ratio_positive_big_int - (inverse_ratio r) (minus_big_int n))) - -let power_num a b = match (a,b) with - (n, (Int i)) -> power_num_int n i -| (n, (Big_int bi)) -> power_num_big_int n bi -| _ -> invalid_arg "power_num" - -let ( **/ ) = power_num - -let is_integer_num = function - Int _ -> true -| Big_int _ -> true -| Ratio r -> is_integer_ratio r - -(* integer_num, floor_num, round_num, ceiling_num rendent des nums *) -let integer_num = function - Int _ as n -> n -| Big_int _ as n -> n -| Ratio r -> num_of_big_int (integer_ratio r) - -and round_num = function - Int _ as n -> n -| Big_int _ as n -> n -| Ratio r -> num_of_big_int (round_ratio r) - -and ceiling_num = function - Int _ as n -> n -| Big_int _ as n -> n -| Ratio r -> num_of_big_int (ceiling_ratio r) - -(* Comparisons on nums *) - -let sign_num = function - Int i -> sign_int i -| Big_int bi -> sign_big_int bi -| Ratio r -> sign_ratio r - -let eq_num a b = match (a,b) with - ((Int int1), (Int int2)) -> int1 = int2 - -| ((Int i), (Big_int bi)) -> eq_big_int (big_int_of_int i) bi -| ((Big_int bi), (Int i)) -> eq_big_int (big_int_of_int i) bi - -| ((Int i), (Ratio r)) -> eq_big_int_ratio (big_int_of_int i) r -| ((Ratio r), (Int i)) -> eq_big_int_ratio (big_int_of_int i) r - -| ((Big_int bi1), (Big_int bi2)) -> eq_big_int bi1 bi2 - -| ((Big_int bi), (Ratio r)) -> eq_big_int_ratio bi r -| ((Ratio r), (Big_int bi)) -> eq_big_int_ratio bi r - -| ((Ratio r1), (Ratio r2)) -> eq_ratio r1 r2 - -let ( =/ ) = eq_num - -let ( <>/ ) a b = not(eq_num a b) - -let compare_num a b = match (a,b) with - ((Int int1), (Int int2)) -> compare_int int1 int2 - -| ((Int i), (Big_int bi)) -> compare_big_int (big_int_of_int i) bi -| ((Big_int bi), (Int i)) -> compare_big_int bi (big_int_of_int i) - -| ((Int i), (Ratio r)) -> compare_big_int_ratio (big_int_of_int i) r -| ((Ratio r), (Int i)) -> -(compare_big_int_ratio (big_int_of_int i) r) - -| ((Big_int bi1), (Big_int bi2)) -> compare_big_int bi1 bi2 - -| ((Big_int bi), (Ratio r)) -> compare_big_int_ratio bi r -| ((Ratio r), (Big_int bi)) -> -(compare_big_int_ratio bi r) - -| ((Ratio r1), (Ratio r2)) -> compare_ratio r1 r2 - -let lt_num num1 num2 = compare_num num1 num2 < 0 -and le_num num1 num2 = compare_num num1 num2 <= 0 -and gt_num num1 num2 = compare_num num1 num2 > 0 -and ge_num num1 num2 = compare_num num1 num2 >= 0 - -let ( </ ) = lt_num -and ( <=/ ) = le_num -and ( >/ ) = gt_num -and ( >=/ ) = ge_num - -let max_num num1 num2 = if lt_num num1 num2 then num2 else num1 -and min_num num1 num2 = if gt_num num1 num2 then num2 else num1 - -(* Coercions with basic types *) - -(* Coercion with int type *) -let int_of_num = function - Int i -> i -| Big_int bi -> int_of_big_int bi -| Ratio r -> int_of_ratio r - -let int_of_num_opt = function - Int i -> Some i -| Big_int bi -> int_of_big_int_opt bi -| Ratio r -> (try Some (int_of_ratio r) with Failure _ -> None) - -and num_of_int i = - if i = monster_int - then Big_int (big_int_of_int i) - else Int i - -(* Coercion with nat type *) -let nat_of_num = function - Int i -> nat_of_int i -| Big_int bi -> nat_of_big_int bi -| Ratio r -> nat_of_ratio r - -and num_of_nat nat = - if (is_nat_int nat 0 (length_nat nat)) - then Int (nth_digit_nat nat 0) - else Big_int (big_int_of_nat nat) - -let nat_of_num_opt x = - try Some (nat_of_num x) with Failure _ -> None - -(* Coercion with big_int type *) -let big_int_of_num = function - Int i -> big_int_of_int i -| Big_int bi -> bi -| Ratio r -> big_int_of_ratio r - -let big_int_of_num_opt x = - try Some (big_int_of_num x) with Failure _ -> None - -let string_of_big_int_for_num bi = - if !approx_printing_flag - then approx_big_int !floating_precision bi - else string_of_big_int bi - -(* Coercion with string type *) - -let string_of_normalized_num = function - Int i -> string_of_int i -| Big_int bi -> string_of_big_int_for_num bi -| Ratio r -> string_of_ratio r -let string_of_num n = - string_of_normalized_num (cautious_normalize_num_when_printing n) - -let num_of_string s = - try - let flag = !normalize_ratio_flag in - normalize_ratio_flag := true; - let r = ratio_of_string s in - normalize_ratio_flag := flag; - if eq_big_int (denominator_ratio r) unit_big_int - then num_of_big_int (numerator_ratio r) - else Ratio r - with Failure _ -> - failwith "num_of_string" - -let num_of_string_opt s = - try Some (num_of_string s) with Failure _ -> None - -(* Coercion with float type *) -let float_of_num = function - Int i -> float i -| Big_int bi -> float_of_big_int bi -| Ratio r -> float_of_ratio r - -let succ_num = function - Int i -> if i = biggest_int - then Big_int (succ_big_int (big_int_of_int i)) - else Int (succ i) -| Big_int bi -> num_of_big_int (succ_big_int bi) -| Ratio r -> Ratio (add_int_ratio 1 r) - -and pred_num = function - Int i -> if i = monster_int - then Big_int (pred_big_int (big_int_of_int i)) - else Int (pred i) -| Big_int bi -> num_of_big_int (pred_big_int bi) -| Ratio r -> Ratio (add_int_ratio (-1) r) - -let abs_num = function - Int i -> if i = monster_int - then Big_int (minus_big_int (big_int_of_int i)) - else Int (abs i) - | Big_int bi -> Big_int (abs_big_int bi) - | Ratio r -> Ratio (abs_ratio r) - -let approx_num_fix n num = approx_ratio_fix n (ratio_of_num num) -and approx_num_exp n num = approx_ratio_exp n (ratio_of_num num) - -let incr_num r = r := succ_num !r -and decr_num r = r := pred_num !r diff --git a/otherlibs/num/num.mli b/otherlibs/num/num.mli deleted file mode 100644 index 4d3793b985..0000000000 --- a/otherlibs/num/num.mli +++ /dev/null @@ -1,191 +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. *) -(* *) -(**************************************************************************) - -(** Operation on arbitrary-precision numbers. - - Numbers (type [num]) are arbitrary-precision rational numbers, - plus the special elements [1/0] (infinity) and [0/0] (undefined). -*) - -open Nat -open Big_int -open Ratio - -(** The type of numbers. *) -type num = - Int of int - | Big_int of big_int - | Ratio of ratio - - -(** {6 Arithmetic operations} *) - - -val ( +/ ) : num -> num -> num -(** Same as {!Num.add_num}.*) - -val add_num : num -> num -> num -(** Addition *) - -val minus_num : num -> num -(** Unary negation. *) - -val ( -/ ) : num -> num -> num -(** Same as {!Num.sub_num}.*) - -val sub_num : num -> num -> num -(** Subtraction *) - -val ( */ ) : num -> num -> num -(** Same as {!Num.mult_num}.*) - -val mult_num : num -> num -> num -(** Multiplication *) - -val square_num : num -> num -(** Squaring *) - -val ( // ) : num -> num -> num -(** Same as {!Num.div_num}.*) - -val div_num : num -> num -> num -(** Division *) - -val quo_num : num -> num -> num -(** Euclidean division: quotient. *) - -val mod_num : num -> num -> num -(** Euclidean division: remainder. *) - -val ( **/ ) : num -> num -> num -(** Same as {!Num.power_num}. *) - -val power_num : num -> num -> num -(** Exponentiation *) - -val abs_num : num -> num -(** Absolute value. *) - -val succ_num : num -> num -(** [succ n] is [n+1] *) - -val pred_num : num -> num -(** [pred n] is [n-1] *) - -val incr_num : num ref -> unit -(** [incr r] is [r:=!r+1], where [r] is a reference to a number. *) - -val decr_num : num ref -> unit -(** [decr r] is [r:=!r-1], where [r] is a reference to a number. *) - -val is_integer_num : num -> bool -(** Test if a number is an integer *) - -(** The four following functions approximate a number by an integer : *) - -val integer_num : num -> num -(** [integer_num n] returns the integer closest to [n]. In case of ties, - rounds towards zero. *) - -val floor_num : num -> num -(** [floor_num n] returns the largest integer smaller or equal to [n]. *) - -val round_num : num -> num -(** [round_num n] returns the integer closest to [n]. In case of ties, - rounds off zero. *) - -val ceiling_num : num -> num -(** [ceiling_num n] returns the smallest integer bigger or equal to [n]. *) - - -val sign_num : num -> int -(** Return [-1], [0] or [1] according to the sign of the argument. *) - -(** {7 Comparisons between numbers} *) - -val ( =/ ) : num -> num -> bool -val ( </ ) : num -> num -> bool -val ( >/ ) : num -> num -> bool -val ( <=/ ) : num -> num -> bool -val ( >=/ ) : num -> num -> bool -val ( <>/ ) : num -> num -> bool -val eq_num : num -> num -> bool -val lt_num : num -> num -> bool -val le_num : num -> num -> bool -val gt_num : num -> num -> bool -val ge_num : num -> num -> bool - -val compare_num : num -> num -> int -(** Return [-1], [0] or [1] if the first argument is less than, - equal to, or greater than the second argument. *) - -val max_num : num -> num -> num -(** Return the greater of the two arguments. *) - -val min_num : num -> num -> num -(** Return the smaller of the two arguments. *) - - -(** {6 Coercions with strings} *) - -val string_of_num : num -> string -(** Convert a number to a string, using fractional notation. *) - -val approx_num_fix : int -> num -> string -(** See {!Num.approx_num_exp}.*) - -val approx_num_exp : int -> num -> string -(** Approximate a number by a decimal. The first argument is the - required precision. The second argument is the number to - approximate. {!Num.approx_num_fix} uses decimal notation; the first - argument is the number of digits after the decimal point. - [approx_num_exp] uses scientific (exponential) notation; the - first argument is the number of digits in the mantissa. *) - -val num_of_string : string -> num -(** Convert a string to a number. - Raise [Failure "num_of_string"] if the given string is not - a valid representation of an integer *) - -val num_of_string_opt: string -> num option -(** Convert a string to a number. - Return [None] if the given string is not - a valid representation of an integer. - - @since 4.05 -*) - -(** {6 Coercions between numerical types} *) - -(* TODO: document the functions below (truncating behavior and error conditions). *) - -val int_of_num : num -> int -val int_of_num_opt: num -> int option -(** @since 4.05.0 *) - -val num_of_int : int -> num -val nat_of_num : num -> nat -val nat_of_num_opt: num -> nat option -(** @since 4.05.0 *) - -val num_of_nat : nat -> num -val num_of_big_int : big_int -> num -val big_int_of_num : num -> big_int -val big_int_of_num_opt: num -> big_int option -(** @since 4.05.0 *) - -val ratio_of_num : num -> ratio -val num_of_ratio : ratio -> num -val float_of_num : num -> float diff --git a/otherlibs/num/ratio.ml b/otherlibs/num/ratio.ml deleted file mode 100644 index 04f9c5e8ff..0000000000 --- a/otherlibs/num/ratio.ml +++ /dev/null @@ -1,619 +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 -open Big_int -open Arith_flags - -(* Definition of the type ratio : - Conventions : - - the denominator is always a positive number - - the sign of n/0 is the sign of n -These convention is automatically respected when a ratio is created with -the create_ratio primitive -*) - -type ratio = { mutable numerator : big_int; - mutable denominator : big_int; - mutable normalized : bool} - -let failwith_zero name = - let s = "infinite or undefined rational number" in - failwith (if String.length name = 0 then s else name ^ " " ^ s) - -let numerator_ratio r = r.numerator -and denominator_ratio r = r.denominator - -let null_denominator r = sign_big_int r.denominator = 0 - -let verify_null_denominator r = - if sign_big_int r.denominator = 0 - then (if !error_when_null_denominator_flag - then (failwith_zero "") - else true) - else false - -let sign_ratio r = sign_big_int r.numerator - -(* Physical normalization of rational numbers *) -(* 1/0, 0/0 and -1/0 are the normalized forms for n/0 numbers *) -let normalize_ratio r = - if r.normalized then r - else if verify_null_denominator r then begin - r.numerator <- big_int_of_int (sign_big_int r.numerator); - r.normalized <- true; - r - end else begin - let p = gcd_big_int r.numerator r.denominator in - if eq_big_int p unit_big_int - then begin - r.normalized <- true; r - end else begin - r.numerator <- div_big_int (r.numerator) p; - r.denominator <- div_big_int (r.denominator) p; - r.normalized <- true; r - end - end - -let cautious_normalize_ratio r = - if (!normalize_ratio_flag) then (normalize_ratio r) else r - -let cautious_normalize_ratio_when_printing r = - if (!normalize_ratio_when_printing_flag) then (normalize_ratio r) else r - -let create_ratio bi1 bi2 = - match sign_big_int bi2 with - -1 -> cautious_normalize_ratio - { numerator = minus_big_int bi1; - denominator = minus_big_int bi2; - normalized = false } - | 0 -> if !error_when_null_denominator_flag - then (failwith_zero "create_ratio") - else cautious_normalize_ratio - { numerator = bi1; denominator = bi2; normalized = false } - | _ -> cautious_normalize_ratio - { numerator = bi1; denominator = bi2; normalized = false } - -let create_normalized_ratio bi1 bi2 = - match sign_big_int bi2 with - -1 -> { numerator = minus_big_int bi1; - denominator = minus_big_int bi2; - normalized = true } -| 0 -> if !error_when_null_denominator_flag - then failwith_zero "create_normalized_ratio" - else { numerator = bi1; denominator = bi2; normalized = true } -| _ -> { numerator = bi1; denominator = bi2; normalized = true } - -let is_normalized_ratio r = r.normalized - -let report_sign_ratio r bi = - if sign_ratio r = -1 - then minus_big_int bi - else bi - -let abs_ratio r = - { numerator = abs_big_int r.numerator; - denominator = r.denominator; - normalized = r.normalized } - -let is_integer_ratio r = - eq_big_int ((normalize_ratio r).denominator) unit_big_int - -(* Operations on rational numbers *) - -let add_ratio r1 r2 = - if !normalize_ratio_flag then begin - let p = gcd_big_int ((normalize_ratio r1).denominator) - ((normalize_ratio r2).denominator) in - if eq_big_int p unit_big_int then - {numerator = add_big_int (mult_big_int (r1.numerator) r2.denominator) - (mult_big_int (r2.numerator) r1.denominator); - denominator = mult_big_int (r1.denominator) r2.denominator; - normalized = true} - else begin - let d1 = div_big_int (r1.denominator) p - and d2 = div_big_int (r2.denominator) p in - let n = add_big_int (mult_big_int (r1.numerator) d2) - (mult_big_int d1 r2.numerator) in - let p' = gcd_big_int n p in - { numerator = div_big_int n p'; - denominator = mult_big_int d1 (div_big_int (r2.denominator) p'); - normalized = true } - end - end else - { numerator = add_big_int (mult_big_int (r1.numerator) r2.denominator) - (mult_big_int (r1.denominator) r2.numerator); - denominator = mult_big_int (r1.denominator) r2.denominator; - normalized = false } - -let minus_ratio r = - { numerator = minus_big_int (r.numerator); - denominator = r.denominator; - normalized = r.normalized } - -let add_int_ratio i r = - ignore (cautious_normalize_ratio r); - { numerator = add_big_int (mult_int_big_int i r.denominator) r.numerator; - denominator = r.denominator; - normalized = r.normalized } - -let add_big_int_ratio bi r = - ignore (cautious_normalize_ratio r); - { numerator = add_big_int (mult_big_int bi r.denominator) r.numerator ; - denominator = r.denominator; - normalized = r.normalized } - -let sub_ratio r1 r2 = add_ratio r1 (minus_ratio r2) - -let mult_ratio r1 r2 = - if !normalize_ratio_flag then begin - let p1 = gcd_big_int ((normalize_ratio r1).numerator) - ((normalize_ratio r2).denominator) - and p2 = gcd_big_int (r2.numerator) r1.denominator in - let (n1, d2) = - if eq_big_int p1 unit_big_int - then (r1.numerator, r2.denominator) - else (div_big_int (r1.numerator) p1, div_big_int (r2.denominator) p1) - and (n2, d1) = - if eq_big_int p2 unit_big_int - then (r2.numerator, r1.denominator) - else (div_big_int r2.numerator p2, div_big_int r1.denominator p2) in - { numerator = mult_big_int n1 n2; - denominator = mult_big_int d1 d2; - normalized = true } - end else - { numerator = mult_big_int (r1.numerator) r2.numerator; - denominator = mult_big_int (r1.denominator) r2.denominator; - normalized = false } - -let mult_int_ratio i r = - if !normalize_ratio_flag then - begin - let p = gcd_big_int ((normalize_ratio r).denominator) (big_int_of_int i) in - if eq_big_int p unit_big_int - then { numerator = mult_big_int (big_int_of_int i) r.numerator; - denominator = r.denominator; - normalized = true } - else { numerator = mult_big_int (div_big_int (big_int_of_int i) p) - r.numerator; - denominator = div_big_int (r.denominator) p; - normalized = true } - end - else - { numerator = mult_int_big_int i r.numerator; - denominator = r.denominator; - normalized = false } - -let mult_big_int_ratio bi r = - if !normalize_ratio_flag then - begin - let p = gcd_big_int ((normalize_ratio r).denominator) bi in - if eq_big_int p unit_big_int - then { numerator = mult_big_int bi r.numerator; - denominator = r.denominator; - normalized = true } - else { numerator = mult_big_int (div_big_int bi p) r.numerator; - denominator = div_big_int (r.denominator) p; - normalized = true } - end - else - { numerator = mult_big_int bi r.numerator; - denominator = r.denominator; - normalized = false } - -let square_ratio r = - ignore (cautious_normalize_ratio r); - { numerator = square_big_int r.numerator; - denominator = square_big_int r.denominator; - normalized = r.normalized } - -let inverse_ratio r = - if !error_when_null_denominator_flag && (sign_big_int r.numerator) = 0 - then failwith_zero "inverse_ratio" - else {numerator = report_sign_ratio r r.denominator; - denominator = abs_big_int r.numerator; - normalized = r.normalized} - -let div_ratio r1 r2 = - mult_ratio r1 (inverse_ratio r2) - -(* Integer part of a rational number *) -(* Odd function *) -let integer_ratio r = - if null_denominator r then failwith_zero "integer_ratio" - else if sign_ratio r = 0 then zero_big_int - else report_sign_ratio r (div_big_int (abs_big_int r.numerator) - (abs_big_int r.denominator)) - -(* Floor of a rational number *) -(* Always less or equal to r *) -let floor_ratio r = - ignore (verify_null_denominator r); - div_big_int (r.numerator) r.denominator - -(* Round of a rational number *) -(* Odd function, 1/2 -> 1 *) -let round_ratio r = - ignore (verify_null_denominator r); - let abs_num = abs_big_int r.numerator in - let bi = div_big_int abs_num r.denominator in - report_sign_ratio r - (if sign_big_int - (sub_big_int - (mult_int_big_int - 2 - (sub_big_int abs_num (mult_big_int (r.denominator) bi))) - r.denominator) = -1 - then bi - else succ_big_int bi) - -let ceiling_ratio r = - if (is_integer_ratio r) - then r.numerator - else succ_big_int (floor_ratio r) - - -(* Comparison operators on rational numbers *) -let eq_ratio r1 r2 = - ignore (normalize_ratio r1); - ignore (normalize_ratio r2); - eq_big_int (r1.numerator) r2.numerator && - eq_big_int (r1.denominator) r2.denominator - -let compare_ratio r1 r2 = - if verify_null_denominator r1 then - let sign_num_r1 = sign_big_int r1.numerator in - if (verify_null_denominator r2) - then - let sign_num_r2 = sign_big_int r2.numerator in - if sign_num_r1 = 1 && sign_num_r2 = -1 then 1 - else if sign_num_r1 = -1 && sign_num_r2 = 1 then -1 - else 0 - else sign_num_r1 - else if verify_null_denominator r2 then - -(sign_big_int r2.numerator) - else match compare_int (sign_big_int r1.numerator) - (sign_big_int r2.numerator) with - 1 -> 1 - | -1 -> -1 - | _ -> if eq_big_int (r1.denominator) r2.denominator - then compare_big_int (r1.numerator) r2.numerator - else compare_big_int - (mult_big_int (r1.numerator) r2.denominator) - (mult_big_int (r1.denominator) r2.numerator) - - -let lt_ratio r1 r2 = compare_ratio r1 r2 < 0 -and le_ratio r1 r2 = compare_ratio r1 r2 <= 0 -and gt_ratio r1 r2 = compare_ratio r1 r2 > 0 -and ge_ratio r1 r2 = compare_ratio r1 r2 >= 0 - -let max_ratio r1 r2 = if lt_ratio r1 r2 then r2 else r1 -and min_ratio r1 r2 = if gt_ratio r1 r2 then r2 else r1 - -let eq_big_int_ratio bi r = - (is_integer_ratio r) && eq_big_int bi r.numerator - -let compare_big_int_ratio bi r = - ignore (normalize_ratio r); - if (verify_null_denominator r) - then -(sign_big_int r.numerator) - else compare_big_int (mult_big_int bi r.denominator) r.numerator - -let lt_big_int_ratio bi r = compare_big_int_ratio bi r < 0 -and le_big_int_ratio bi r = compare_big_int_ratio bi r <= 0 -and gt_big_int_ratio bi r = compare_big_int_ratio bi r > 0 -and ge_big_int_ratio bi r = compare_big_int_ratio bi r >= 0 - -(* Coercions *) - -(* Coercions with type int *) -let int_of_ratio r = - if ((is_integer_ratio r) && (is_int_big_int r.numerator)) - then (int_of_big_int r.numerator) - else failwith "integer argument required" - -and ratio_of_int i = - { numerator = big_int_of_int i; - denominator = unit_big_int; - normalized = true } - -(* Coercions with type nat *) -let ratio_of_nat nat = - { numerator = big_int_of_nat nat; - denominator = unit_big_int; - normalized = true } - -and nat_of_ratio r = - ignore (normalize_ratio r); - if not (is_integer_ratio r) then - failwith "nat_of_ratio" - else if sign_big_int r.numerator > -1 then - nat_of_big_int (r.numerator) - else failwith "nat_of_ratio" - -(* Coercions with type big_int *) -let ratio_of_big_int bi = - { numerator = bi; denominator = unit_big_int; normalized = true } - -and big_int_of_ratio r = - ignore (normalize_ratio r); - if is_integer_ratio r - then r.numerator - else failwith "big_int_of_ratio" - -let div_int_ratio i r = - ignore (verify_null_denominator r); - mult_int_ratio i (inverse_ratio r) - -let div_ratio_int r i = - div_ratio r (ratio_of_int i) - -let div_big_int_ratio bi r = - ignore (verify_null_denominator r); - mult_big_int_ratio bi (inverse_ratio r) - -let div_ratio_big_int r bi = - div_ratio r (ratio_of_big_int bi) - -(* Functions on type string *) -(* giving floating point approximations of rational numbers *) - -(* Compares strings that contains only digits, have the same length, - from index i to index i + l *) -let rec compare_num_string s1 s2 i len = - if i >= len then 0 else - let c1 = int_of_char s1.[i] - and c2 = int_of_char s2.[i] in - match compare_int c1 c2 with - | 0 -> compare_num_string s1 s2 (succ i) len - | c -> c;; - -(* Position of the leading digit of the decimal expansion *) -(* of a strictly positive rational number *) -(* if the decimal expansion of a non null rational r is equal to *) -(* sigma for k=-P to N of r_k*10^k then msd_ratio r = N *) -(* Nota : for a big_int we have msd_ratio = nums_digits_big_int -1 *) - -(* Tests if s has only zeros characters from index i to index lim *) -let rec only_zeros s i lim = - i >= lim || s.[i] == '0' && only_zeros s (succ i) lim;; - -(* Nota : for a big_int we have msd_ratio = nums_digits_big_int -1 *) -let msd_ratio r = - ignore (cautious_normalize_ratio r); - if null_denominator r then failwith_zero "msd_ratio" - else if sign_big_int r.numerator == 0 then 0 - else begin - let str_num = string_of_big_int r.numerator - and str_den = string_of_big_int r.denominator in - let size_num = String.length str_num - and size_den = String.length str_den in - let size_min = min size_num size_den in - let m = size_num - size_den in - let cmp = compare_num_string str_num str_den 0 size_min in - match cmp with - | 1 -> m - | -1 -> pred m - | _ -> - if m >= 0 then m else - if only_zeros str_den size_min size_den then m - else pred m - end -;; - -(* Decimal approximations of rational numbers *) - -(* Approximation with fix decimal point *) -(* This is an odd function and the last digit is round off *) -(* Format integer_part . decimal_part_with_n_digits *) -let approx_ratio_fix n r = - (* Don't need to normalize *) - if (null_denominator r) then failwith_zero "approx_ratio_fix" - else - let sign_r = sign_ratio r in - if sign_r = 0 - then "+0" (* r = 0 *) - else - (* r.numerator and r.denominator are not null numbers - s1 contains one more digit than desired for the round off operation *) - if n >= 0 then begin - let s1 = - string_of_nat - (nat_of_big_int - (div_big_int - (base_power_big_int - 10 (succ n) (abs_big_int r.numerator)) - r.denominator)) in - (* Round up and add 1 in front if needed *) - let s2 = - if round_futur_last_digit (Bytes.unsafe_of_string s1) 0 - (String.length s1) - then "1" ^ s1 - else s1 in - let l2 = String.length s2 - 1 in - (* if s2 without last digit is xxxxyyy with n 'yyy' digits: - <sign> xxxx . yyy - if s2 without last digit is yy with <= n digits: - <sign> 0 . 0yy *) - if l2 > n then begin - let s = Bytes.make (l2 + 2) '0' in - Bytes.set s 0 (if sign_r = -1 then '-' else '+'); - String.blit s2 0 s 1 (l2 - n); - Bytes.set s (l2 - n + 1) '.'; - String.blit s2 (l2 - n) s (l2 - n + 2) n; - Bytes.unsafe_to_string s - end else begin - let s = Bytes.make (n + 3) '0' in - Bytes.set s 0 (if sign_r = -1 then '-' else '+'); - Bytes.set s 2 '.'; - String.blit s2 0 s (n + 3 - l2) l2; - Bytes.unsafe_to_string s - end - end else begin - (* Dubious; what is this code supposed to do? *) - let s = string_of_big_int - (div_big_int - (abs_big_int r.numerator) - (base_power_big_int - 10 (-n) r.denominator)) in - let len = succ (String.length s) in - let s' = Bytes.make len '0' in - Bytes.set s' 0 (if sign_r = -1 then '-' else '+'); - String.blit s 0 s' 1 (pred len); - Bytes.unsafe_to_string s' - end - -(* Number of digits of the decimal representation of an int *) -let num_decimal_digits_int n = - String.length (string_of_int n) - -(* Approximation with floating decimal point *) -(* This is an odd function and the last digit is round off *) -(* Format (+/-)(0. n_first_digits e msd)/(1. n_zeros e (msd+1) *) -let approx_ratio_exp n r = - (* Don't need to normalize *) - if (null_denominator r) then failwith_zero "approx_ratio_exp" - else if n <= 0 then invalid_arg "approx_ratio_exp" - else - let sign_r = sign_ratio r - and i = ref (n + 3) in - if sign_r = 0 then - String.concat "" ["+0."; String.make n '0'; "e0"] - else - let msd = msd_ratio (abs_ratio r) in - let k = n - msd in - let s = - (let nat = nat_of_big_int - (if k < 0 - then - div_big_int (abs_big_int r.numerator) - (base_power_big_int 10 (- k) - r.denominator) - else - div_big_int (base_power_big_int - 10 k (abs_big_int r.numerator)) - r.denominator) in - string_of_nat nat) in - if round_futur_last_digit (Bytes.unsafe_of_string s) 0 - (String.length s) - then - let m = num_decimal_digits_int (succ msd) in - let str = Bytes.make (n + m + 4) '0' in - (String.blit (if sign_r = -1 then "-1." else "+1.") 0 str 0 3); - Bytes.set str !i ('e'); - incr i; - (if m = 0 - then Bytes.set str !i '0' - else String.blit (string_of_int (succ msd)) 0 str !i m); - Bytes.unsafe_to_string str - else - let m = num_decimal_digits_int (succ msd) - and p = n + 3 in - let str = Bytes.make (succ (m + p)) '0' in - (String.blit (if sign_r = -1 then "-0." else "+0.") 0 str 0 3); - (String.blit s 0 str 3 n); - Bytes.set str p 'e'; - (if m = 0 - then Bytes.set str (succ p) '0' - else (String.blit (string_of_int (succ msd)) 0 str (succ p) m)); - Bytes.unsafe_to_string str - -(* String approximation of a rational with a fixed number of significant *) -(* digits printed *) -let float_of_rational_string r = - let s = approx_ratio_exp !floating_precision r in - if String.get s 0 = '+' - then (String.sub s 1 (pred (String.length s))) - else s - -(* Coercions with type string *) -let string_of_ratio r = - ignore (cautious_normalize_ratio_when_printing r); - if !approx_printing_flag - then float_of_rational_string r - else string_of_big_int r.numerator ^ "/" ^ string_of_big_int r.denominator - -(* XL: j'ai puissamment simplifie "ratio_of_string" en virant la notation - scientifique. - | I have strongly simplified "ratio_of_string" by deleting scientific notation -*) - -let ratio_of_string s = - try - let n = String.index s '/' in - create_ratio (sys_big_int_of_string s 0 n) - (sys_big_int_of_string s (n+1) (String.length s - n - 1)) - with Not_found -> - { numerator = big_int_of_string s; - denominator = unit_big_int; - normalized = true } - -(* Coercion with type float *) - -let float_of_ratio r = - let p = r.numerator and q = r.denominator in - (* Special cases 0/0, 0/q and p/0 *) - if sign_big_int q = 0 then begin - match sign_big_int p with - | 0 -> nan - | 1 -> infinity - | -1 -> neg_infinity - | _ -> assert false - end - else if sign_big_int p = 0 then 0.0 - else begin - let np = num_bits_big_int p and nq = num_bits_big_int q in - if np <= 53 && nq <= 53 then - (* p and q convert to floats exactly; use FP division to get the - correctly-rounded result. *) - Int64.to_float (int64_of_big_int p) - /. Int64.to_float (int64_of_big_int q) - else begin - let ap = abs_big_int p in - (* |p| is in [2^(np-1), 2^np) - q is in [2^(nq-1), 2^nq) - hence |p|/q is in (2^(np-nq-1), 2^(np-nq+1)). - We define n such that |p|/q*2^n is in [2^54, 2^56). - >= 2^54 so that the round to odd technique applies. - < 2^56 so that the integral part is representable as an int64. *) - let n = 55 - (np - nq) in - (* Scaling |p|/q by 2^n *) - let (p', q') = - if n >= 0 - then (shift_left_big_int ap n, q) - else (ap, shift_left_big_int q (-n)) in - (* Euclidean division of p' by q' *) - let (quo, rem) = quomod_big_int p' q' in - (* quo is the integral part of |p|/q*2^n - rem/q' is the fractional part. *) - (* Round quo to float *) - let f = round_big_int_to_float quo (sign_big_int rem = 0) in - (* Apply exponent *) - let f = ldexp f (-n) in - (* Apply sign *) - if sign_big_int p < 0 then -. f else f - end - end - - -let power_ratio_positive_int r n = - create_ratio (power_big_int_positive_int (r.numerator) n) - (power_big_int_positive_int (r.denominator) n) - -let power_ratio_positive_big_int r bi = - create_ratio (power_big_int_positive_big_int (r.numerator) bi) - (power_big_int_positive_big_int (r.denominator) bi) diff --git a/otherlibs/num/ratio.mli b/otherlibs/num/ratio.mli deleted file mode 100644 index 4a765051d7..0000000000 --- a/otherlibs/num/ratio.mli +++ /dev/null @@ -1,93 +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. *) -(* *) -(**************************************************************************) - -(** Operation on rational numbers. - - This module is used to support the implementation of {!Num} and - should not be called directly. *) - -open Nat -open Big_int - -(* Rationals (type [ratio]) are arbitrary-precision rational numbers, - plus the special elements [1/0] (infinity) and [0/0] (undefined). - In constrast with numbers (type [num]), the special cases of - small integers and big integers are not optimized specially. *) - -type ratio - -(**/**) - -val null_denominator : ratio -> bool -val numerator_ratio : ratio -> big_int -val denominator_ratio : ratio -> big_int -val sign_ratio : ratio -> int -val normalize_ratio : ratio -> ratio -val cautious_normalize_ratio : ratio -> ratio -val cautious_normalize_ratio_when_printing : ratio -> ratio -val create_ratio : big_int -> big_int -> ratio (* assumes nothing *) -val create_normalized_ratio : big_int -> big_int -> ratio - (* assumes normalized argument *) -val is_normalized_ratio : ratio -> bool -val report_sign_ratio : ratio -> big_int -> big_int -val abs_ratio : ratio -> ratio -val is_integer_ratio : ratio -> bool -val add_ratio : ratio -> ratio -> ratio -val minus_ratio : ratio -> ratio -val add_int_ratio : int -> ratio -> ratio -val add_big_int_ratio : big_int -> ratio -> ratio -val sub_ratio : ratio -> ratio -> ratio -val mult_ratio : ratio -> ratio -> ratio -val mult_int_ratio : int -> ratio -> ratio -val mult_big_int_ratio : big_int -> ratio -> ratio -val square_ratio : ratio -> ratio -val inverse_ratio : ratio -> ratio -val div_ratio : ratio -> ratio -> ratio -val integer_ratio : ratio -> big_int -val floor_ratio : ratio -> big_int -val round_ratio : ratio -> big_int -val ceiling_ratio : ratio -> big_int -val eq_ratio : ratio -> ratio -> bool -val compare_ratio : ratio -> ratio -> int -val lt_ratio : ratio -> ratio -> bool -val le_ratio : ratio -> ratio -> bool -val gt_ratio : ratio -> ratio -> bool -val ge_ratio : ratio -> ratio -> bool -val max_ratio : ratio -> ratio -> ratio -val min_ratio : ratio -> ratio -> ratio -val eq_big_int_ratio : big_int -> ratio -> bool -val compare_big_int_ratio : big_int -> ratio -> int -val lt_big_int_ratio : big_int -> ratio -> bool -val le_big_int_ratio : big_int -> ratio -> bool -val gt_big_int_ratio : big_int -> ratio -> bool -val ge_big_int_ratio : big_int -> ratio -> bool -val int_of_ratio : ratio -> int -val ratio_of_int : int -> ratio -val ratio_of_nat : nat -> ratio -val nat_of_ratio : ratio -> nat -val ratio_of_big_int : big_int -> ratio -val big_int_of_ratio : ratio -> big_int -val div_int_ratio : int -> ratio -> ratio -val div_ratio_int : ratio -> int -> ratio -val div_big_int_ratio : big_int -> ratio -> ratio -val div_ratio_big_int : ratio -> big_int -> ratio -val approx_ratio_fix : int -> ratio -> string -val approx_ratio_exp : int -> ratio -> string -val float_of_rational_string : ratio -> string -val string_of_ratio : ratio -> string -val ratio_of_string : string -> ratio -val float_of_ratio : ratio -> float -val power_ratio_positive_int : ratio -> int -> ratio -val power_ratio_positive_big_int : ratio -> big_int -> ratio |