[Toss-devel-svn] SF.net SVN: toss:[1677] trunk/Toss
Status: Beta
Brought to you by:
lukaszkaiser
Menu
▾
▴
[Toss-devel-svn] SF.net SVN: toss:[1677] trunk/Toss
|
From: <luk...@us...> - 2012-02-18 02:07:15
|
Revision: 1677
http://toss.svn.sourceforge.net/toss/?rev=1677&view=rev
Author: lukaszkaiser
Date: 2012-02-18 02:07:02 +0000 (Sat, 18 Feb 2012)
Log Message:
-----------
Replacing Num.
Modified Paths:
--------------
trunk/Toss/Makefile
trunk/Toss/Server/Tests.ml
trunk/Toss/Solver/Num/Integers.ml
trunk/Toss/Solver/Num/Integers.mli
trunk/Toss/Solver/Num/IntegersTest.ml
trunk/Toss/Solver/Num/Naturals.ml
trunk/Toss/Solver/Num/NaturalsTest.ml
trunk/Toss/Solver/Num/NumbersTest.ml
trunk/Toss/Solver/Num/Rationals.ml
Modified: trunk/Toss/Makefile
===================================================================
--- trunk/Toss/Makefile 2012-02-17 02:36:56 UTC (rev 1676)
+++ trunk/Toss/Makefile 2012-02-18 02:07:02 UTC (rev 1677)
@@ -47,8 +47,8 @@
# TODO: Hard-coded path to js_of_ocaml.
OCB_LFLAG=-lflags -I,/usr/local/lib/ocaml/3.12.0/js_of_ocaml,-I,+oUnit,-I,+js_of_ocaml,-I,+site-lib/oUnit,-I,+site-lib/js_of_ocaml
OCB_CFLAG=-cflags -I,/usr/local/lib/ocaml/3.12.0/js_of_ocaml,-I,+oUnit,-I,+js_of_ocaml,-I,+site-lib/oUnit,-I,+site-lib/js_of_ocaml,-g
-OCB_LIB=-libs str,nums,unix,oUnit
-OCB_LIBJS=-libs str,js_of_ocaml
+OCB_LIB=-libs str,unix,oUnit
+OCB_LIBJS=-libs js_of_ocaml
OCB_PP=-pp "camlp4o -I /usr/local/lib/ocaml/3.12.0 -I /opt/local/lib/ocaml/site-lib ../caml_extensions/pa_let_try.cmo ../caml_extensions/pa_log.cmo pa_macro.cmo js_of_ocaml/pa_js.cmo"
OCB_PPJS=-pp "camlp4o -unsafe -I /usr/local/lib/ocaml/3.12.0 -I /opt/local/lib/ocaml/site-lib ../caml_extensions/pa_let_try.cmo ../caml_extensions/pa_log.cmo pa_macro.cmo -DJAVASCRIPT js_of_ocaml/pa_js.cmo"
OCAMLBUILD=ocamlbuild -log build.log -j 8 -menhir ../menhir_conf \
Modified: trunk/Toss/Server/Tests.ml
===================================================================
--- trunk/Toss/Server/Tests.ml 2012-02-17 02:36:56 UTC (rev 1676)
+++ trunk/Toss/Server/Tests.ml 2012-02-18 02:07:02 UTC (rev 1677)
@@ -14,9 +14,9 @@
let solver_tests = "Solver", [
"NaturalsTest", [NaturalsTest.tests];
- "IntegersTest", [IntegersTest.tests];
+ "IntegersTest", [IntegersTest.tests; IntegersTest.bigtests];
"RationalsTest", [RationalsTest.tests];
- "NumbersTest", [NumbersTest.tests];
+ "NumbersTest", [NumbersTest.tests; NumbersTest.bigtests];
"StructureTest", [StructureTest.tests];
"AssignmentsTest", [AssignmentsTest.tests];
"SolverTest", [SolverTest.tests; SolverTest.bigtests];
Modified: trunk/Toss/Solver/Num/Integers.ml
===================================================================
--- trunk/Toss/Solver/Num/Integers.ml 2012-02-17 02:36:56 UTC (rev 1676)
+++ trunk/Toss/Solver/Num/Integers.ml 2012-02-18 02:07:02 UTC (rev 1677)
@@ -13,7 +13,7 @@
open MiscNum
-open Naturals.N
+open Naturals
type big_int =
{ sign : int;
@@ -76,27 +76,31 @@
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
+let string_of_big_int bi =
+ if bi.sign = -1
+ then "-" ^ string_of_nat bi.abs_value
+ else string_of_nat bi.abs_value
+
+
(* Operations on big_int *)
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 *)
+ 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) (bi2.abs_value) 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;
- add_nat res (bi1.abs_value);
- 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;
- add_nat res (bi2.abs_value);
- res)}
-
+ | -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;
+ add_nat res (bi1.abs_value);
+ 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;
+ add_nat res (bi2.abs_value);
+ res)}
else (* Subtract absolute values if signs are different *)
match compare_nat (bi1.abs_value) (bi2.abs_value) with
| 0 -> zero_big_int
@@ -113,12 +117,12 @@
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;
- incr_nat res)
- else set_digit_nat res 0 (abs i));
- res }
+ if i = monster_int then
+ let res = (create_nat 2) in
+ (set_digit_nat res 0 biggest_int;
+ incr_nat res; res)
+ else let res = (create_nat 1) in (set_digit_nat res 0 (abs i); res)
+ }
let big_int_of_nat nat =
let length = num_digits_nat nat in
@@ -146,34 +150,34 @@
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
+ 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 (mult_nat res (bi2.abs_value) (bi1.abs_value); res)
- else (mult_nat res (bi1.abs_value) (bi2.abs_value); res)
+ if size_bi2 > size_bi1
+ then (mult_nat res (bi2.abs_value) (bi1.abs_value); res)
+ else (mult_nat res (bi1.abs_value) (bi2.abs_value); res)
}
(* (quotient, rest) 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) (bi2.abs_value) with
- -1 -> (* 1/2 -> 0, reste 1, -1/2 -> -1, reste 1 *)
- (* 1/-2 -> 0, reste 1, -1/-2 -> 1, reste 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)
+ let size_bi1 = num_digits_big_int bi1
+ and size_bi2 = num_digits_big_int bi2 in
+ match compare_nat (bi1.abs_value) (bi2.abs_value) with
+ | -1 -> (* 1/2 -> 0, reste 1, -1/2 -> -1, reste 1 *)
+ (* 1/-2 -> 0, reste 1, -1/-2 -> 1, reste 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 bi1_negatif = (bi1.sign = -1) in
let size_q =
if bi1_negatif
then succ (max (succ (size_bi1 - size_bi2)) 1)
@@ -236,26 +240,25 @@
let 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
+ 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 }
+ { sign = 1; abs_value = bi1.abs_value }
else
- { sign = 1;
- abs_value =
- match compare_nat (bi1.abs_value) (bi2.abs_value) with
- | 0 -> bi1.abs_value
- | 1 ->
- let res = copy_nat (bi1.abs_value) 0 size_bi1 in
- let len = gcd_nat res (bi2.abs_value) in
- copy_nat res 0 len
- | _ ->
- let res = copy_nat (bi2.abs_value) 0 size_bi2 in
- let len = gcd_nat res (bi1.abs_value) in
- copy_nat res 0 len
- }
+ { sign = 1;
+ abs_value =
+ match compare_nat (bi1.abs_value) (bi2.abs_value) with
+ | 0 -> bi1.abs_value
+ | 1 ->
+ let res = copy_nat (bi1.abs_value) 0 size_bi1 in
+ let len = gcd_nat res (bi2.abs_value) in
+ copy_nat res 0 len
+ | _ ->
+ let res = copy_nat (bi2.abs_value) 0 size_bi2 in
+ let len = gcd_nat res (bi1.abs_value) in
+ copy_nat res 0 len
+ }
(* Coercion operators *)
@@ -264,11 +267,9 @@
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 monster_nat with
- | 0 -> bi.sign == -1
- | -1 -> true
- | _ -> false;;
+ num_digits_big_int bi == 1 || (
+ num_digits_big_int bi == 2 && bi.sign == -1 &&
+ compare_nat bi.abs_value monster_nat == 0)
let int_of_big_int bi =
try let n = int_of_nat bi.abs_value in
@@ -286,12 +287,6 @@
(* 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 =
if len < 1 then failwith "sys_big_int_of_string";
let n = nat_of_string s ofs len in
Modified: trunk/Toss/Solver/Num/Integers.mli
===================================================================
--- trunk/Toss/Solver/Num/Integers.mli 2012-02-17 02:36:56 UTC (rev 1676)
+++ trunk/Toss/Solver/Num/Integers.mli 2012-02-18 02:07:02 UTC (rev 1677)
@@ -98,7 +98,7 @@
val big_int_of_int : int -> big_int
(** Convert a small integer to a big integer. *)
-val big_int_of_nat : Naturals.N.nat -> big_int
+val big_int_of_nat : Naturals.nat -> big_int
(** Convert a natural to a big integer. *)
val is_int_big_int : big_int -> bool
@@ -121,7 +121,7 @@
(** {2 For internal use} *)
-val nat_of_big_int : big_int -> Naturals.N.nat
+val nat_of_big_int : big_int -> Naturals.nat
val approx_big_int: int -> big_int -> string
val base_power_big_int: int -> int -> big_int -> big_int
val round_futur_last_digit : string -> int -> int -> bool
Modified: trunk/Toss/Solver/Num/IntegersTest.ml
===================================================================
--- trunk/Toss/Solver/Num/IntegersTest.ml 2012-02-17 02:36:56 UTC (rev 1676)
+++ trunk/Toss/Solver/Num/IntegersTest.ml 2012-02-18 02:07:02 UTC (rev 1677)
@@ -1,14 +1,22 @@
open OUnit
open Integers
-let eq_bool (b1, b2) = assert_equal ~printer:string_of_bool b1 b2
+let eq_bool ?i (b1, b2) = match i with
+ | None -> assert_equal ~printer:string_of_bool b1 b2
+ | Some j -> assert_equal ~msg:(Printf.sprintf "test %i" j)
+ ~printer:string_of_bool b1 b2
let eq_int (i1, i2) = assert_equal ~printer:string_of_int i1 i2
let eq_string (s1, s2) = assert_equal ~printer:(fun x -> x) s1 s2
-let eq_big_int (bi1, bi2) = eq_bool (Integers.eq_big_int bi1 bi2, true)
+let eq_big_int ?i (bi1, bi2) =
+ let eq = Integers.eq_big_int bi1 bi2 in
+ if eq then () else
+ let is = match i with None -> "" | Some j -> Printf.sprintf "test %i: " j in
+ let msg = is ^ (string_of_big_int bi1) ^ " <> " ^ (string_of_big_int bi2) in
+ assert_equal ~msg ~printer:string_of_bool eq true
let failwith_test f x except =
try let _ = ignore (f x) in eq_string ("worked", "failed") with
- e -> eq_bool (e = except, true)
+ e -> if e = except then () else raise e
let rec gcd_int i1 i2 = if i2 = 0 then abs i1 else gcd_int i2 (i1 mod i2)
@@ -17,6 +25,19 @@
let biggest_int = monster_int - 1
let least_int = - biggest_int
+let pi_100_digits =
+"3141592653 :10
+5897932384 :20
+6264338327 :30
+9502884197 :40
+1693993751 :50
+0582097494 :60
+4592307816 :70
+4062862089 :80
+9862803482 :90
+5342117067 :100
+"
+
let pi_1000_digits =
"3141592653 :10
5897932384 :20
@@ -120,6 +141,55 @@
9216420198 :1000
"
+let pi_digits n_digits =
+ (* Pi digits computed with the streaming algorithm given on pages 4, 6
+ & 7 of "Unbounded Spigot Algorithms for the Digits of Pi", Jeremy
+ Gibbons, August 2004. *)
+ let ( !$ ) = big_int_of_int
+ and ( +$ ) = add_big_int
+ and ( *$ ) = mult_big_int in
+
+ let zero = zero_big_int
+ and one = unit_big_int
+ and three = !$ 3
+ and four = !$ 4
+ and ten = !$ 10
+ and neg_ten = !$(-10) in
+
+ (* Linear Fractional (aka Moebius) Transformations *)
+ let floor_ev (q, r, s, t) x = div_big_int (q *$ x +$ r) (s *$ x +$ t) in
+ let unit = (one, zero, zero, one) in
+ let comp (q, r, s, t) (q', r', s', t') =
+ (q *$ q' +$ r *$ s', q *$ r' +$ r *$ t',
+ s *$ q' +$ t *$ s', s *$ r' +$ t *$ t') in
+
+ let next z = floor_ev z three in
+ let safe z n = Integers.eq_big_int n (floor_ev z four) in
+ let prod z n = comp (ten, neg_ten *$ n, zero, one) z in
+ let cons z k = let den = 2 * k + 1 in
+ comp z (!$ k, !$(2 * den), zero, !$ den) in
+
+ let rec digit k z n row col acc =
+ if n > 0 then
+ let y = next z in
+ if safe z y then
+ if col = 10 then (
+ let row = row + 10 in
+ digit k (prod z y) (n - 1) row 1
+ ((Printf.sprintf "\t:%i\n%s" row (string_of_big_int y)) :: acc)
+ ) else (
+ digit k (prod z y) (n - 1) row (col + 1)
+ ((string_of_big_int y) :: acc)
+ )
+ else digit (k + 1) (cons z k) n row col acc
+ else
+ (Printf.sprintf "%*s\t:%i\n" (10 - col) "" (row + col)) :: acc in
+
+ let digits n = digit 1 unit n 0 0 [] in
+ String.concat "" (List.rev (digits n_digits))
+
+
+
let tests = "Integers" >::: [
"compare_big_int" >::
(fun () ->
@@ -140,76 +210,76 @@
"add_big_int" >::
(fun () ->
- eq_big_int (add_big_int zero_big_int zero_big_int, zero_big_int);
- eq_big_int (add_big_int zero_big_int (big_int_of_int 1),
- big_int_of_int 1);
- eq_big_int (add_big_int (big_int_of_int 1) zero_big_int,
- big_int_of_int 1);
- eq_big_int (add_big_int zero_big_int (big_int_of_int (-1)),
- big_int_of_int (-1));
- eq_big_int (add_big_int (big_int_of_int (-1)) zero_big_int,
- big_int_of_int (-1));
- eq_big_int (add_big_int (big_int_of_int 1) (big_int_of_int 1),
- big_int_of_int 2);
- eq_big_int (add_big_int (big_int_of_int 1) (big_int_of_int 2),
- big_int_of_int 3);
- eq_big_int (add_big_int (big_int_of_int 2) (big_int_of_int 1),
- big_int_of_int 3);
- eq_big_int (add_big_int (big_int_of_int (-1)) (big_int_of_int (-1)),
- big_int_of_int (-2));
- eq_big_int (add_big_int (big_int_of_int (-1)) (big_int_of_int (-2)),
- big_int_of_int (-3));
- eq_big_int (add_big_int (big_int_of_int (-2)) (big_int_of_int (-1)),
- big_int_of_int (-3));
- eq_big_int (add_big_int (big_int_of_int 1) (big_int_of_int (-1)),
- zero_big_int);
- eq_big_int (add_big_int (big_int_of_int (-1)) (big_int_of_int 1),
- zero_big_int);
- eq_big_int (add_big_int (big_int_of_int 1) (big_int_of_int (-2)),
- big_int_of_int (-1));
- eq_big_int (add_big_int (big_int_of_int (-2)) (big_int_of_int 1),
- big_int_of_int (-1));
- eq_big_int (add_big_int (big_int_of_int (-1)) (big_int_of_int 2),
- big_int_of_int 1);
- eq_big_int (add_big_int (big_int_of_int 2) (big_int_of_int (-1)),
- big_int_of_int 1);
+ eq_big_int ~i:1 (add_big_int zero_big_int zero_big_int, zero_big_int);
+ eq_big_int ~i:2 (add_big_int zero_big_int (big_int_of_int 1),
+ big_int_of_int 1);
+ eq_big_int ~i:3 (add_big_int (big_int_of_int 1) zero_big_int,
+ big_int_of_int 1);
+ eq_big_int ~i:4 (add_big_int zero_big_int (big_int_of_int (-1)),
+ big_int_of_int (-1));
+ eq_big_int ~i:5 (add_big_int (big_int_of_int (-1)) zero_big_int,
+ big_int_of_int (-1));
+ eq_big_int ~i:6 (add_big_int (big_int_of_int 1) (big_int_of_int 1),
+ big_int_of_int 2);
+ eq_big_int ~i:7 (add_big_int (big_int_of_int 1) (big_int_of_int 2),
+ big_int_of_int 3);
+ eq_big_int ~i:8 (add_big_int (big_int_of_int 2) (big_int_of_int 1),
+ big_int_of_int 3);
+ eq_big_int ~i:9 (add_big_int (big_int_of_int (-1)) (big_int_of_int (-1)),
+ big_int_of_int (-2));
+ eq_big_int ~i:10 (add_big_int (big_int_of_int (-1)) (big_int_of_int (-2)),
+ big_int_of_int (-3));
+ eq_big_int ~i:11 (add_big_int (big_int_of_int (-2)) (big_int_of_int (-1)),
+ big_int_of_int (-3));
+ eq_big_int ~i:12 (add_big_int (big_int_of_int 1) (big_int_of_int (-1)),
+ zero_big_int);
+ eq_big_int ~i:13 (add_big_int (big_int_of_int (-1)) (big_int_of_int 1),
+ zero_big_int);
+ eq_big_int ~i:14 (add_big_int (big_int_of_int 1) (big_int_of_int (-2)),
+ big_int_of_int (-1));
+ eq_big_int ~i:15 (add_big_int (big_int_of_int (-2)) (big_int_of_int 1),
+ big_int_of_int (-1));
+ eq_big_int ~i:16 (add_big_int (big_int_of_int (-1)) (big_int_of_int 2),
+ big_int_of_int 1);
+ eq_big_int ~i:17 (add_big_int (big_int_of_int 2) (big_int_of_int (-1)),
+ big_int_of_int 1);
);
"sub_big_int" >::
(fun () ->
- eq_big_int (sub_big_int zero_big_int zero_big_int, zero_big_int);
- eq_big_int (sub_big_int zero_big_int (big_int_of_int 1),
- big_int_of_int (-1));
- eq_big_int (sub_big_int (big_int_of_int 1) zero_big_int,
- big_int_of_int 1);
- eq_big_int (sub_big_int zero_big_int (big_int_of_int (-1)),
- big_int_of_int 1);
- eq_big_int (sub_big_int (big_int_of_int (-1)) zero_big_int,
- big_int_of_int (-1));
- eq_big_int (sub_big_int (big_int_of_int 1) (big_int_of_int 1),
- zero_big_int);
- eq_big_int (sub_big_int (big_int_of_int 1) (big_int_of_int 2),
- big_int_of_int (-1));
- eq_big_int (sub_big_int (big_int_of_int 2) (big_int_of_int 1),
- big_int_of_int 1);
- eq_big_int (sub_big_int (big_int_of_int (-1)) (big_int_of_int (-1)),
- zero_big_int);
- eq_big_int (sub_big_int (big_int_of_int (-1)) (big_int_of_int (-2)),
- big_int_of_int 1);
- eq_big_int (sub_big_int (big_int_of_int (-2)) (big_int_of_int (-1)),
- big_int_of_int (-1));
- eq_big_int (sub_big_int (big_int_of_int 1) (big_int_of_int (-1)),
- big_int_of_int 2);
- eq_big_int (sub_big_int (big_int_of_int (-1)) (big_int_of_int 1),
- big_int_of_int (-2));
- eq_big_int (sub_big_int (big_int_of_int 1) (big_int_of_int (-2)),
- big_int_of_int 3);
- eq_big_int (sub_big_int (big_int_of_int (-2)) (big_int_of_int 1),
- big_int_of_int (-3));
- eq_big_int (sub_big_int (big_int_of_int (-1)) (big_int_of_int 2),
- big_int_of_int (-3));
- eq_big_int (sub_big_int (big_int_of_int 2) (big_int_of_int (-1)),
- big_int_of_int 3);
+ eq_big_int ~i:1 (sub_big_int zero_big_int zero_big_int, zero_big_int);
+ eq_big_int ~i:2 (sub_big_int zero_big_int (big_int_of_int 1),
+ big_int_of_int (-1));
+ eq_big_int ~i:3 (sub_big_int (big_int_of_int 1) zero_big_int,
+ big_int_of_int 1);
+ eq_big_int ~i:4 (sub_big_int zero_big_int (big_int_of_int (-1)),
+ big_int_of_int 1);
+ eq_big_int ~i:5 (sub_big_int (big_int_of_int (-1)) zero_big_int,
+ big_int_of_int (-1));
+ eq_big_int ~i:6 (sub_big_int (big_int_of_int 1) (big_int_of_int 1),
+ zero_big_int);
+ eq_big_int ~i:7 (sub_big_int (big_int_of_int 1) (big_int_of_int 2),
+ big_int_of_int (-1));
+ eq_big_int ~i:8 (sub_big_int (big_int_of_int 2) (big_int_of_int 1),
+ big_int_of_int 1);
+ eq_big_int ~i:9 (sub_big_int (big_int_of_int (-1)) (big_int_of_int (-1)),
+ zero_big_int);
+ eq_big_int ~i:10 (sub_big_int (big_int_of_int (-1)) (big_int_of_int (-2)),
+ big_int_of_int 1);
+ eq_big_int ~i:11 (sub_big_int (big_int_of_int (-2)) (big_int_of_int (-1)),
+ big_int_of_int (-1));
+ eq_big_int ~i:12 (sub_big_int (big_int_of_int 1) (big_int_of_int (-1)),
+ big_int_of_int 2);
+ eq_big_int ~i:13 (sub_big_int (big_int_of_int (-1)) (big_int_of_int 1),
+ big_int_of_int (-2));
+ eq_big_int ~i:14 (sub_big_int (big_int_of_int 1) (big_int_of_int (-2)),
+ big_int_of_int 3);
+ eq_big_int ~i:15 (sub_big_int (big_int_of_int (-2)) (big_int_of_int 1),
+ big_int_of_int (-3));
+ eq_big_int ~i:16 (sub_big_int (big_int_of_int (-1)) (big_int_of_int 2),
+ big_int_of_int (-3));
+ eq_big_int ~i:17 (sub_big_int (big_int_of_int 2) (big_int_of_int (-1)),
+ big_int_of_int 3);
);
"mult_int_big_int" >::
@@ -240,43 +310,43 @@
(fun () ->
let (quotient, modulo) =
quomod_big_int (big_int_of_int 1) (big_int_of_int 1) in
- eq_big_int (quotient, big_int_of_int 1);
- eq_big_int (modulo, zero_big_int);
+ eq_big_int ~i:1 (quotient, big_int_of_int 1);
+ eq_big_int ~i:2 (modulo, zero_big_int);
let (quotient, modulo) =
quomod_big_int (big_int_of_int 1) (big_int_of_int (-1)) in
- eq_big_int (quotient, big_int_of_int (-1));
- eq_big_int (modulo, zero_big_int);
+ eq_big_int ~i:3 (quotient, big_int_of_int (-1));
+ eq_big_int ~i:4 (modulo, zero_big_int);
let (quotient, modulo) =
quomod_big_int (big_int_of_int (-1)) (big_int_of_int 1) in
- eq_big_int (quotient, big_int_of_int (-1));
- eq_big_int (modulo, zero_big_int);
+ eq_big_int ~i:5 (quotient, big_int_of_int (-1));
+ eq_big_int ~i:6 (modulo, zero_big_int);
let (quotient, modulo) =
quomod_big_int (big_int_of_int 3) (big_int_of_int 2) in
- eq_big_int (quotient, big_int_of_int 1);
- eq_big_int (modulo, big_int_of_int 1);
+ eq_big_int ~i:7 (quotient, big_int_of_int 1);
+ eq_big_int ~i:8 (modulo, big_int_of_int 1);
let (quotient, modulo) =
quomod_big_int (big_int_of_int 5) (big_int_of_int 3) in
- eq_big_int (quotient, big_int_of_int 1);
- eq_big_int (modulo, big_int_of_int 2);
+ eq_big_int ~i:9 (quotient, big_int_of_int 1);
+ eq_big_int ~i:10 (modulo, big_int_of_int 2);
let (quotient, modulo) =
quomod_big_int (big_int_of_int (-5)) (big_int_of_int 3) in
- eq_big_int (quotient, big_int_of_int (-2));
- eq_big_int (modulo, big_int_of_int 1);
+ eq_big_int ~i:11 (quotient, big_int_of_int (-2));
+ eq_big_int ~i:12 (modulo, big_int_of_int 1);
let (quotient, modulo) =
quomod_big_int (big_int_of_int 1) (big_int_of_int 2) in
- eq_big_int (quotient, zero_big_int);
- eq_big_int (modulo, big_int_of_int 1);
+ eq_big_int ~i:13 (quotient, zero_big_int);
+ eq_big_int ~i:14 (modulo, big_int_of_int 1);
let (quotient, modulo) =
quomod_big_int (big_int_of_int (-1)) (big_int_of_int 3) in
- eq_big_int (quotient, minus_big_int unit_big_int);
- eq_big_int (modulo, big_int_of_int 2);
+ eq_big_int ~i:14 (quotient, minus_big_int unit_big_int);
+ eq_big_int ~i:15 (modulo, big_int_of_int 2);
failwith_test
(quomod_big_int (big_int_of_int 1)) zero_big_int
@@ -284,23 +354,23 @@
let (quotient, modulo) =
quomod_big_int (big_int_of_int 10) (big_int_of_int 20) in
- eq_big_int (quotient, big_int_of_int 0);
- eq_big_int (modulo, big_int_of_int 10);
+ eq_big_int ~i:16 (quotient, big_int_of_int 0);
+ eq_big_int ~i:17 (modulo, big_int_of_int 10);
let (quotient, modulo) =
quomod_big_int (big_int_of_int (-10)) (big_int_of_int 20) in
- eq_big_int (quotient, big_int_of_int (-1));
- eq_big_int (modulo, big_int_of_int 10);
+ eq_big_int ~i:18 (quotient, big_int_of_int (-1));
+ eq_big_int ~i:19 (modulo, big_int_of_int 10);
let (quotient, modulo) =
quomod_big_int (big_int_of_int 10) (big_int_of_int (-20)) in
- eq_big_int (quotient, big_int_of_int 0);
- eq_big_int (modulo, big_int_of_int 10);
+ eq_big_int ~i:20 (quotient, big_int_of_int 0);
+ eq_big_int ~i:21 (modulo, big_int_of_int 10);
let (quotient, modulo) =
quomod_big_int (big_int_of_int (-10)) (big_int_of_int (-20)) in
- eq_big_int (quotient, big_int_of_int 1);
- eq_big_int (modulo, big_int_of_int 10);
+ eq_big_int ~i:22 (quotient, big_int_of_int 1);
+ eq_big_int ~i:23 (modulo, big_int_of_int 10);
);
"gcd_big_int" >::
@@ -353,27 +423,24 @@
"is_int_big_int" >::
(fun () ->
- eq_bool (is_int_big_int (big_int_of_int 1), true);
- eq_bool (is_int_big_int (big_int_of_int (-1)), true);
- eq_bool (is_int_big_int (add_big_int (big_int_of_int 1)
- (big_int_of_int biggest_int)),
- false);
+ eq_bool ~i:1 (is_int_big_int (big_int_of_int 1), true);
+ eq_bool ~i:2 (is_int_big_int (big_int_of_int (-1)), true);
+ eq_bool ~i:3 (is_int_big_int (add_big_int (big_int_of_int 1)
+ (big_int_of_int biggest_int)), false);
eq_int (int_of_big_int (big_int_of_int monster_int), monster_int);
- eq_bool (is_int_big_int (big_int_of_string (string_of_int biggest_int)),
- true);
- eq_bool (is_int_big_int (big_int_of_string (string_of_int least_int)),
- true);
- eq_bool (is_int_big_int (big_int_of_string (string_of_int monster_int)),
- true);
- eq_bool (is_int_big_int (add_big_int (big_int_of_int 1)
- (big_int_of_int (biggest_int))),
- false);
- eq_bool (is_int_big_int (add_big_int (big_int_of_int 1)
- (big_int_of_int (biggest_int))),
- false);
- eq_bool (is_int_big_int
- (minus_big_int (big_int_of_string(string_of_int monster_int))),
- false);
+ eq_bool ~i:4 (is_int_big_int (big_int_of_string
+ (string_of_int biggest_int)), true);
+ eq_bool ~i:5 (is_int_big_int (big_int_of_string
+ (string_of_int least_int)), true);
+ eq_bool ~i:6 (is_int_big_int (big_int_of_int monster_int), true);
+ eq_bool ~i:7 (is_int_big_int (big_int_of_string
+ (string_of_int monster_int)), true);
+ eq_bool ~i:8 (is_int_big_int (add_big_int (big_int_of_int 1)
+ (big_int_of_int (biggest_int))), false);
+ eq_bool ~i:9 (is_int_big_int (add_big_int (big_int_of_int 1)
+ (big_int_of_int (biggest_int))), false);
+ eq_bool ~i:10 (is_int_big_int (minus_big_int (
+ big_int_of_string(string_of_int monster_int))), false);
);
"string_of_big_int" >::
@@ -383,13 +450,13 @@
"big_int_of_string" >::
(fun () ->
- eq_big_int (big_int_of_string "1", big_int_of_int 1);
- eq_big_int (big_int_of_string "-1", big_int_of_int (-1));
- eq_big_int (big_int_of_string "0", zero_big_int);
+ eq_big_int ~i:1 (big_int_of_string "1", big_int_of_int 1);
+ eq_big_int ~i:2 (big_int_of_string "-1", big_int_of_int (-1));
+ eq_big_int ~i:3 (big_int_of_string "0", zero_big_int);
failwith_test big_int_of_string "sdjdkfighdgf" (Failure "invalid digit");
- eq_big_int (big_int_of_string "123", big_int_of_int 123);
- eq_big_int (big_int_of_string "3456", big_int_of_int 3456);
- eq_big_int (big_int_of_string "-3456", big_int_of_int (-3456));
+ eq_big_int ~i:4 (big_int_of_string "123", big_int_of_int 123);
+ eq_big_int ~i:5 (big_int_of_string "3456", big_int_of_int 3456);
+ eq_big_int ~i:6 (big_int_of_string "-3456", big_int_of_int (-3456));
let implode = List.fold_left (^) "" in
let l = List.rev [
@@ -414,18 +481,19 @@
] in
let bi1 = big_int_of_string (implode (List.rev l)) in
let bi2 = big_int_of_string (implode (List.rev ("3" :: List.tl l))) in
- eq_big_int (bi1, (add_big_int (mult_big_int bi2 (big_int_of_string "10"))
- (big_int_of_string "2")));
+ eq_big_int ~i:7 (bi1, (add_big_int (mult_big_int bi2
+ (big_int_of_string "10"))
+ (big_int_of_string "2")));
);
"power_base_int" >::
(fun () ->
- eq_big_int (big_int_of_nat (Naturals.N.power_base_int 10 0),unit_big_int);
- eq_big_int (big_int_of_nat (Naturals.N.power_base_int 10 8),
+ eq_big_int (big_int_of_nat (Naturals.power_base_int 10 0),unit_big_int);
+ eq_big_int (big_int_of_nat (Naturals.power_base_int 10 8),
big_int_of_int 100000000);
- eq_big_int (big_int_of_nat(Naturals.N.power_base_int 2 (length_of_int+2)),
- big_int_of_nat (let nat = Naturals.N.make_nat 2 in
- Naturals.N.set_digit_nat nat 1 1;
+ eq_big_int (big_int_of_nat(Naturals.power_base_int 2 (length_of_int)),
+ big_int_of_nat (let nat = Naturals.make_nat 2 in
+ Naturals.set_digit_nat nat 1 1;
nat));
);
@@ -438,53 +506,16 @@
big_int_of_int 1230);
);
- "pi digits" >::
+ "pi 100 digits" >::
(fun () ->
- (* Pi digits computed with the streaming algorithm given on pages 4, 6
- & 7 of "Unbounded Spigot Algorithms for the Digits of Pi", Jeremy
- Gibbons, August 2004. *)
- let ( !$ ) = big_int_of_int
- and ( +$ ) = add_big_int
- and ( *$ ) = mult_big_int in
+ eq_string (pi_digits 100, pi_100_digits);
+ );
+]
- let zero = zero_big_int
- and one = unit_big_int
- and three = !$ 3
- and four = !$ 4
- and ten = !$ 10
- and neg_ten = !$(-10) in
-
- (* Linear Fractional (aka Moebius) Transformations *)
- let floor_ev (q, r, s, t) x = div_big_int (q *$ x +$ r) (s *$ x +$ t) in
- let unit = (one, zero, zero, one) in
- let comp (q, r, s, t) (q', r', s', t') =
- (q *$ q' +$ r *$ s', q *$ r' +$ r *$ t',
- s *$ q' +$ t *$ s', s *$ r' +$ t *$ t') in
-
- let next z = floor_ev z three in
- let safe z n = Integers.eq_big_int n (floor_ev z four) in
- let prod z n = comp (ten, neg_ten *$ n, zero, one) z in
- let cons z k = let den = 2 * k + 1 in
- comp z (!$ k, !$(2 * den), zero, !$ den) in
-
- let rec digit k z n row col acc =
- if n > 0 then
- let y = next z in
- if safe z y then
- if col = 10 then (
- let row = row + 10 in
- digit k (prod z y) (n - 1) row 1
- ((Printf.sprintf "\t:%i\n%s" row (string_of_big_int y)) :: acc)
- ) else (
- digit k (prod z y) (n - 1) row (col + 1)
- ((string_of_big_int y) :: acc)
- )
- else digit (k + 1) (cons z k) n row col acc
- else
- (Printf.sprintf "%*s\t:%i\n" (10 - col) "" (row + col)) :: acc in
-
- let digits n = digit 1 unit n 0 0 [] in
- eq_string (String.concat "" (List.rev (digits 1000)), pi_1000_digits);
+let bigtests = "IntegersBig" >::: [
+ "pi 1000 digits" >::
+ (fun () ->
+ eq_string (pi_digits 1000, pi_1000_digits);
);
]
Modified: trunk/Toss/Solver/Num/Naturals.ml
===================================================================
--- trunk/Toss/Solver/Num/Naturals.ml 2012-02-17 02:36:56 UTC (rev 1676)
+++ trunk/Toss/Solver/Num/Naturals.ml 2012-02-18 02:07:02 UTC (rev 1677)
@@ -1,201 +1,225 @@
-IFDEF JAVASCRIPT THEN
open Array
-module N = (struct
- type nat = int array
- let max_int_length = Sys.word_size - 2 (* should be even *)
+type nat = int array
- let create_nat s = make s 0
+let max_int_length = Sys.word_size - 2 (* should be even *)
+let max_power10_int = 1000000000
+let sprint_full_length_int i = Printf.sprintf "%.9i" i
- let set_to_zero_nat s i1 i2 =
- for i = i1 to i2 do s.(i) <- 0; done
+let create_nat s = make s 0
- let make_nat len =
- if len < 0 then invalid_arg "make_nat" else create_nat len
+let set_to_zero_nat s i1 i2 =
+ for i = i1 to i2 do s.(i) <- 0; done
- let blit_nat n1 i1 n2 i2 i3 = blit n1 i1 n2 i2 i3
+let make_nat len =
+ if len < 0 then invalid_arg "make_nat" else
+ if len = 0 then create_nat 1 else create_nat len
- let copy_nat nat offset length =
- let res = create_nat length in blit nat offset res 0 length; res
+let blit_nat n1 i1 n2 i2 i3 = blit n2 i2 n1 i1 i3
- let set_digit_nat n d x = n.(d) <- x
+let copy_nat nat offset length =
+ let res = create_nat length in blit nat offset res 0 length; res
- let num_digits_nat n =
- let l = ref ((length n) - 1) in
- while (!l >= 0 && n.(!l) = 0) do l := !l - 1 done;
- !l + 1
+let set_digit_nat n d x = n.(d) <- x
- let is_zero_nat n i1 i2 = num_digits_nat (copy_nat n i1 i2) = 0
+let one_nat = make 1 1
- let shrink n =
- let m = num_digits_nat n in if m = length n then n else
- let res = make m 0 in blit n 0 res 0 m; res
+let num_digits_nat n =
+ let l = ref ((length n) - 1) in
+ while (!l >= 0 && n.(!l) = 0) do l := !l - 1 done;
+ !l + 1
- let int_of_nat n =
- if num_digits_nat n > 1 then failwith "int_of_nat" else n.(0)
- let nat_of_int i = make 1 i
+let is_zero_nat n i1 i2 = num_digits_nat (copy_nat n i1 i2) = 0
- let compare_nat n m =
- let rec compare_from i =
- let ni = if i < length n then n.(i) else 0 in
- let mi = if i < length m then m.(i) else 0 in
- if ni < mi then -1 else if ni > mi then 1 else
- if i = 0 then 0 else compare_from (i-1) in
- compare_from ((max (length n) (length m)) - 1)
+let int_of_nat n =
+ if num_digits_nat n > 1 then failwith "int_of_nat" else n.(0)
+let nat_of_int i = make 1 i
- let add_nat n x = (* n := n + x *)
- let rec add_carry i carry =
- if i >= length n then (if carry <> 0 then failwith "overflow") else
- if i >= length x then (
- let res = n.(i) + carry in
- if res >= 0 then n.(i) <- res else (
- n.(i) <- 0; add_carry (i+1) 1
- )
+let compare_nat n m =
+ let rec compare_from i =
+ let ni = if i < length n then n.(i) else 0 in
+ let mi = if i < length m then m.(i) else 0 in
+ if ni < mi then -1 else if ni > mi then 1 else
+ if i = 0 then 0 else compare_from (i-1) in
+ compare_from ((max (length n) (length m)) - 1)
+
+let add_nat_off off n x = (* n := n + (x shifted by off) *)
+ let rec add_carry i carry =
+ if i + off >= length n then (if carry <> 0 then failwith "overflow") else
+ if i >= length x then (
+ let res = n.(i+off) + carry in
+ if res >= 0 then n.(i+off) <- res else (
+ n.(i+off) <- 0; add_carry (i+1) 1
+ )
+ ) else (
+ let res = n.(i+off) + x.(i) + carry in
+ if res >= 0 then (
+ n.(i+off) <- res;
+ add_carry (i+1) 0
) else (
- let res = n.(i) + x.(i) + carry in
- if res >= 0 then (
- n.(i) <- res;
- add_carry (i+1) 0
- ) else (
- let mid = n.(i) - max_int - 1 in
- n.(i) <- mid + x.(i) + carry;
- add_carry (i+1) 1
- )
- ) in
- add_carry 0 0
+ let mid = n.(i+off) - max_int - 1 in
+ n.(i+off) <- mid + x.(i) + carry;
+ add_carry (i+1) 1
+ )
+ ) in
+ add_carry 0 0
- let incr_nat n = add_nat n (make 1 1)
+let add_nat n x = add_nat_off 0 n x
- let sub_nat n x = (* n := n - x *)
- let rec sub_carry i carry =
- if i >= length n then (if carry <> 0 then failwith "sub too big") else
- if i >= length x then (
- let res = n.(i) + carry in
- if res >= 0 then n.(i) <- res else (
- sub_carry (i+1) (-1)
- )
- ) else
- let res = n.(i) - x.(i) + carry in
- if res >= 0 then (
- n.(i) <- res;
- sub_carry (i+1) 0;
- ) else (
- n.(i) <- res + 1;
- n.(i) <- n.(i) + max_int;
- sub_carry (i+1) (-1)
- ) in
- sub_carry 0 0
+let incr_nat n = add_nat n one_nat
- let half_int = 1 lsl (max_int_length / 2)
+let sub_nat n x = (* n := n - x *)
+ let rec sub_carry i carry =
+ if i >= length n then (if carry <> 0 then failwith "sub too big") else
+ if i >= length x then (
+ let res = n.(i) + carry in
+ if res >= 0 then n.(i) <- res else (
+ sub_carry (i+1) (-1)
+ )
+ ) else
+ let res = n.(i) - x.(i) + carry in
+ if res >= 0 then (
+ n.(i) <- res;
+ sub_carry (i+1) 0;
+ ) else (
+ n.(i) <- res + 1;
+ n.(i) <- n.(i) + max_int;
+ sub_carry (i+1) (-1)
+ ) in
+ sub_carry 0 0
- let mult_digit_nat_off off n x i = (* n := x * i shift by offset *)
- let add_to_pos j x = (* n.(j) <- n.(j)+x; unsafe-add overflow to n.(j+1) *)
- if n.(j) + x >= 0 then n.(j) <- n.(j) + x else (
- let mid = x - max_int - 1 in
- n.(j) <- n.(j) + mid;
- n.(j+1) <- n.(j+1) + 1;
- ) in
- let i0, i1 = i mod half_int, i / half_int in
- let rec mult_digit j =
- if (j >= length x) then () else (
- let x0, x1 = x.(j) mod half_int, x.(j) / half_int in
- let res0, resa, resb = x0 * i0, x1 * i0, x0 * i1 in
- let res0a, res1a = (resa mod half_int) * half_int, resa / half_int in
- let res0b, res1b = (resb mod half_int) * half_int, resb / half_int in
- let next = x1 * i1 + res1a + res1b in
- if next > 0 then n.(off+j+1) <- next;
- add_to_pos (off+j) res0;
- add_to_pos (off+j) res0a;
- add_to_pos (off+j) res0b;
- mult_digit (j+1);
- ) in
- for i = 0 to off do n.(i) <- 0 done;
- mult_digit 0
+let half_int = 1 lsl (max_int_length / 2)
- let mult_digit_nat n x i = (* n := x*i *)
- fill n 0 (length n) 0;
- mult_digit_nat_off 0 n x i
+let one_arr = make 1 0
+let add_nat_off_digit off n digit =
+ one_arr.(0) <- digit;
+ add_nat_off off n one_arr
- let mult_nat n x1 x2 = (* n := x1 * x2 *)
- fill n 0 (length n) 0;
- let interim = make (length n) 0 in
- for j = 0 to (length x2) - 1 do
- mult_digit_nat_off j interim x1 x2.(j);
- add_nat n interim;
+let mult_digit_nat_off off n x i = (* n += x*i shift by offset *)
+ let i0, i1 = i mod half_int, i / half_int in
+ let rec mult_digit j =
+ if (j >= length x) then () else (
+ let x0, x1 = x.(j) mod half_int, x.(j) / half_int in
+ let res0, resa, resb = x0 * i0, x1 * i0, x0 * i1 in
+ let res0a, res1a = (resa mod half_int) * half_int, resa / half_int in
+ let res0b, res1b = (resb mod half_int) * half_int, resb / half_int in
+ let next = x1 * i1 + res1a + res1b in
+ if next > 0 then add_nat_off_digit (off+j+1) n next;
+ add_nat_off_digit (off+j) n res0;
+ add_nat_off_digit (off+j) n res0a;
+ add_nat_off_digit (off+j) n res0b;
+ mult_digit (j+1);
+ ) in
+ if i = 0 then () else if i = 1 then add_nat_off off n x else
+ mult_digit 0
+
+let mult_digit_nat n x i = (* n += x*i *)
+ mult_digit_nat_off 0 n x i
+
+let mult_nat n x1 x2 = (* n += x1 * x2 *)
+ let (nd1, nd2) = num_digits_nat x1, num_digits_nat x2 in
+ if nd1 = 0 || nd2 = 0 then () else (
+ for j = 0 to nd2 - 1 do
+ mult_digit_nat_off j n x1 x2.(j);
done
+ )
- let div_nat n xin = (* n := n / x *)
- let res, x = make_nat (length n), shrink xin in
- let lx = length x in
- let rec approx_backshift i n m =
- if m = max_int then
- let (d, b) = approx_backshift 0 n ((m / 2) + 1) in (d, b+1)
- else if m+1 > n * (1 lsl i) then
- approx_backshift (i+1) n m
- else (n * (1 lsl i)) / (m+1), i in
- let rec div_subs cur =
- match compare_nat cur x with
- | y when y < 0 -> ()
- | 0 -> add_nat res (make 1 1)
- | _->
- let l = length cur - lx in
- let (i, b) = approx_backshift 0 cur.((length cur) - 1) x.(lx - 1) in
- let xmult, resmult = make (l+1) 0, make ((length cur)+1) 0 in
- if b = 0 then (
- xmult.(l) <- i;
+let prealloc_len = 1000
+let prealloc_res, prealloc_resmult = make prealloc_len 0, make prealloc_len 0
+
+let div_nat_fn ln n x = (* put (n := n mod x) and return (n / x) *)
+ let (res, resmult) =
+ if prealloc_len > ln then (
+ fill prealloc_res 0 (ln+1) 0;
+ fill prealloc_resmult 0 (ln+1) 0;
+ prealloc_res, prealloc_resmult
+ ) else (make_nat ln, make_nat ln) in
+ let lx = num_digits_nat x in
+ let rec approx_backshift i n m add =
+ if m = max_int then
+ let (d, b) = approx_backshift 0 n ((m / 2) + 1) add in (d, b+1)
+ else
+ let shn = n * (1 lsl i) in
+ if shn > 0 && m+add > shn then
+ approx_backshift (i+1) n m add
+ else if shn > 0 then
+ shn / (m+add), i
+ else 1, i in
+ let rec div_subs cur =
+ match compare_nat cur x with
+ | y when y < 0 -> ()
+ | 0 -> add_nat res one_nat; fill cur 0 ln 0
+ | _ ->
+ let lc = num_digits_nat cur in
+ let l = lc - lx in
+ let (i, b) = approx_backshift 0 cur.(lc-1) x.(lx-1)
+ (if lx = 1 then 0 else 1) in
+ if b = 0 then (
+ add_nat_off_digit l res i;
+ mult_digit_nat_off l resmult x i;
+ ) else (
+ if l = 0 then (
+ add_nat_off_digit l res 1;
+ mult_digit_nat_off l resmult x 1;
) else (
- xmult.(l-1) <- i * (1 lsl (max_int_length - b));
- );
- add_nat res xmult;
- mult_nat resmult xmult x;
- sub_nat cur resmult;
- div_subs (shrink cur) in
- div_subs (shrink n);
- for i = 0 to (length n) - 1 do n.(i) <- res.(i) done
+ let d = i * (1 lsl (max_int_length - b)) in
+ add_nat_off_digit (l-1) res d;
+ mult_digit_nat_off (l-1) resmult x d;
+ )
+ );
+ sub_nat cur resmult;
+ fill resmult 0 ln 0;
+ div_subs cur in
+ div_subs n;
+ res
- let mod_nat a b = (* a-(b*[a/b]) *)
- let div = copy a in
- div_nat div b;
- let prod = make (length a) 0 in
- mult_nat prod div b;
- let res = copy a in
- sub_nat res prod;
- res
+(* After div, n contains
+ - in the (length x) least significant digits, the remainder
+ - in the (length n)-(length x) most significant digits, the quotient *)
+let div_nat n x =
+ let ln = num_digits_nat n in
+ let lnplus = if ln = length n then ln else ln+1 in
+ let quo, lx = div_nat_fn ln n x, length x in
+ for i = lx to lnplus - 1 do
+ n.(i) <- quo.(i - lx)
+ done
+
+let gcd_nat a b = (* set a := gcd a b, return the length *)
+ let rec gcd_n a b = if num_digits_nat b = 0 then a else (
+ ignore (div_nat_fn (num_digits_nat a) a b); gcd_n b a) in
+ let res = gcd_n a (copy b) in
+ blit res 0 a 0 (min (length res) (length a));
+ num_digits_nat a
- let rec gcd_nat_fn a b =
- if compare_nat a b < 0 then gcd_nat_fn b a else
- if num_digits_nat b = 0 then a else gcd_nat_fn b (mod_nat a b)
+let shrink ?max n =
+ let nd = num_digits_nat n in
+ let m = match max with None -> nd | Some m -> min m nd in
+ if m = length n then n else
+ if m = 0 then make 1 0 else let res = make m 0 in blit n 0 res 0 m; res
- let gcd_nat a b = (* set a := gcd a b, return length of a *)
- let res = shrink (gcd_nat_fn a b) in
- fill a 0 (length a) 0;
- blit res 0 a 0 (length res);
- length res
+let rec power_base_int n i = (* n ^ i *)
+ if i < 0 then invalid_arg "negative power" else
+ if i = 0 then make 1 1 else if i = 1 then make 1 n else
+ let r = power_base_int n (i / 2) in
+ let rsq = make (2 * (length r)) 0 in
+ mult_nat rsq r r;
+ if i mod 2 = 0 then shrink rsq else
+ let m = make (2 * (length r) + 1) 0 in
+ (mult_nat m rsq (make 1 n); shrink m)
+
+let string_of_nat n =
+ let rec string_rec m =
+ let lm = length m in
+ if lm = 1 then string_of_int m.(0) else (
+ let quo = div_nat_fn (num_digits_nat m) m (make 1 max_power10_int) in
+ let s = string_rec (shrink ~max:lm quo) in
+ s ^ (sprint_full_length_int m.(0))
+ ) in
+ if num_digits_nat n = 0 then "0" else string_rec (copy (shrink n))
- let rec power_base_int n i = (* n ^ i *)
- if i < 0 then invalid_arg "negative power" else
- if i = 0 then make 1 1 else if i = 1 then make 1 n else
- let r = power_base_int n (i / 2) in
- let rsq = make (2 * (length r)) 0 in
- mult_nat rsq r r;
- if i mod 2 = 0 then shrink rsq else
- let m = make (2 * (length r) + 1) 0 in
- (mult_nat m rsq (make 1 n); shrink m)
-
- let string_of_nat n =
- let rec string_rec m =
- if length m = 1 then string_of_int m.(0) else (
- let mcp, mdiv = copy m, make (length m) 0 in
- div_nat m (make 1 1000000);
- mult_nat mdiv m (make 1 1000000);
- sub_nat mcp mdiv;
- let s = string_rec (shrink m) in
- s ^ (Printf.sprintf "%.6i" mcp.(0))
- ) in
- if num_digits_nat n = 0 then "0" else string_rec (copy (shrink n))
-
- let max_int_str_len = String.length (string_of_int max_int)
- let rec nat_of_string s ofs len =
+let max_int_str_len = String.length (string_of_int max_int)
+let rec nat_of_string s ofs len =
+ try
if len < max_int_str_len then
make 1 (int_of_string (String.sub s ofs len))
else (
@@ -208,54 +232,4 @@
add_nat res n;
res
)
-end)
-ELSE
-module N = (struct
- type nat = Nat.nat * int (* store number size *)
-
- let create_nat s = (Nat.create_nat s, s)
- let set_to_zero_nat (s, l) i1 i2 = Nat.set_to_zero_nat s i1 i2
- 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
-
- let blit_nat (n1, l1) i1 (n2, l2) i2 i3 = Nat.blit_nat n1 i1 n2 i2 i3
- let copy_nat nat offset length =
- let res = create_nat (length) in
- blit_nat res 0 nat offset length;
- res
-
- let set_digit_nat (n, l) d x = Nat.set_digit_nat n d x
- let num_digits_nat (n, l) =
- Nat.num_digits_nat n 0 (Nat.length_nat n)
-
- let is_zero_nat (n, _) i1 i2 = Nat.is_zero_nat n i1 i2
- let int_of_nat (n, l) = Nat.int_of_nat n
- let nat_of_int i = (Nat.nat_of_int i, 1)
- let incr_nat (n, l) = ignore (Nat.incr_nat n 0 l 1)
-
- let add_nat (n, ln) (x, lx) =
- ignore (Nat.add_nat n 0 ln x 0 lx 0)
-
- let sub_nat (n, ln) (x, lx) = ignore (Nat.sub_nat n 0 ln x 0 lx 1)
-
- let mult_digit_nat (n, ln) (x, lx) i =
- ignore (Nat.mult_digit_nat n 0 ln x 0 lx (Nat.nat_of_int i) 0)
-
- let mult_nat (n, ln) (x1, l1) (x2, l2) =
- ignore (Nat.mult_nat n 0 ln x1 0 l1 x2 0 l2)
-
- let div_nat (n, ln) (x, lx) = Nat.div_nat n 0 ln x 0 lx
-
- let compare_nat (n, ln) (m, lm) = Nat.compare_nat n 0 ln m 0 lm
- let gcd_nat (n, ln) (m, lm) = Nat.gcd_nat n 0 ln m 0 lm
- let string_of_nat (n, _) = Nat.string_of_nat n
- let nat_of_string s ofs len =
- let n = Nat.sys_nat_of_string 10 s ofs len in
- (n, Nat.num_digits_nat n 0 (Nat.length_nat n))
-
- let power_base_int i j =
- let n = Nat.power_base_int i j in
- (n, Nat.num_digits_nat n 0 (Nat.length_nat n))
-end)
-ENDIF
+ with Failure s -> failwith "invalid digit"
Modified: trunk/Toss/Solver/Num/NaturalsTest.ml
===================================================================
--- trunk/Toss/Solver/Num/NaturalsTest.ml 2012-02-17 02:36:56 UTC (rev 1676)
+++ trunk/Toss/Solver/Num/NaturalsTest.ml 2012-02-18 02:07:02 UTC (rev 1677)
@@ -1,5 +1,5 @@
open OUnit
-open Naturals.N
+open Naturals
let eq_bool (b1, b2) = assert_equal ~printer:string_of_bool b2 b1
let eq_int (i1, i2) = assert_equal ~printer:string_of_int i2 i1
@@ -42,11 +42,12 @@
incr_nat n;
equal_nat (n, nat_of_int 2);
- (*let n = create_nat 2 in
+ let n = create_nat 2 in
set_digit_nat n 0 max_int;
+ set_digit_nat n 1 0;
incr_nat n;
sub_nat n (nat_of_int 2);
- equal_nat (n, nat_of_int (max_int - 1));*)
+ equal_nat (n, nat_of_int (max_int - 1));
);
"is_zero_nat" >::
Modified: trunk/Toss/Solver/Num/NumbersTest.ml
===================================================================
--- trunk/Toss/Solver/Num/NumbersTest.ml 2012-02-17 02:36:56 UTC (rev 1676)
+++ trunk/Toss/Solver/Num/NumbersTest.ml 2012-02-18 02:07:02 UTC (rev 1677)
@@ -17,7 +17,51 @@
let biggest_int = monster_int - 1
let least_int = - biggest_int
+let pi_digits n_digits =
+ (* Pi digits computed with the streaming algorithm given on pages 4, 6
+ & 7 of "Unbounded Spigot Algorithms for the Digits of Pi", Jeremy
+ Gibbons, August 2004. *)
+ let zero = num_of_int 0
+ and one = num_of_int 1
+ and three = num_of_int 3
+ and four = num_of_int 4
+ and ten = num_of_int 10
+ and neg_ten = num_of_int (-10) in
+ (* Linear Fractional (aka Moebius) Transformations *)
+ let floor_ev (q, r, s, t) x = quo_num (q */ x +/ r) (s */ x +/ t) in
+ let unit = (one, zero, zero, one) in
+ let comp (q, r, s, t) (q', r', s', t') =
+ (q */ q' +/ r */ s', q */ r' +/ r */ t',
+ s */ q' +/ t */ s', s */ r' +/ t */ t') in
+
+ let next z = floor_ev z three in
+ let safe z n = (n =/ (floor_ev z four)) in
+ let prod z n = comp (ten, neg_ten */ n, zero, one) z in
+ let cons z k =
+ let den = 2 * k + 1 in
+ comp z (num_of_int k, num_of_int (2 * den), zero, num_of_int den) in
+
+ let rec digit k z n row col acc =
+ if n > 0 then
+ let y = next z in
+ if safe z y then
+ if col = 10 then (
+ let row = row + 10 in
+ digit k (prod z y) (n - 1) row 1
+ ((Printf.sprintf "\t:%i\n%s" row (string_of_num y)) :: acc)
+ ) else (
+ digit k (prod z y) (n - 1) row (col + 1)
+ ((string_of_num y) :: acc)
+ )
+ else digit (k + 1) (cons z k) n row col acc
+ else
+ (Printf.sprintf "%*s\t:%i\n" (10 - col) "" (row + col)) :: acc in
+
+ let digits n = digit 1 unit n 0 0 [] in
+ String.concat "" (List.rev (digits n_digits))
+
+
let tests = "Numbers" >::: [
"add_num" >::
(fun () ->
@@ -134,51 +178,16 @@
failwith_test num_of_string ("frlshjkurty") (Failure "num_of_string");
);
- "pi digits" >::
+ "pi 100 digits" >::
(fun () ->
- (* Pi digits computed with the streaming algorithm given on pages 4, 6
- & 7 of "Unbounded Spigot Algorithms for the Digits of Pi", Jeremy
- Gibbons, August 2004. *)
- let zero = num_of_int 0
- and one = num_of_int 1
- and three = num_of_int 3
- and four = num_of_int 4
- and ten = num_of_int 10
- and neg_ten = num_of_int (-10) in
+ eq_string (pi_digits 100, IntegersTest.pi_100_digits);
+ );
+]
- (* Linear Fractional (aka Moebius) Transformations *)
- let floor_ev (q, r, s, t) x = quo_num (q */ x +/ r) (s */ x +/ t) in
- let unit = (one, zero, zero, one) in
- let comp (q, r, s, t) (q', r', s', t') =
- (q */ q' +/ r */ s', q */ r' +/ r */ t',
- s */ q' +/ t */ s', s */ r' +/ t */ t') in
-
- let next z = floor_ev z three in
- let safe z n = (n =/ (floor_ev z four)) in
- let prod z n = comp (ten, neg_ten */ n, zero, one) z in
- let cons z k =
- let den = 2 * k + 1 in
- comp z (num_of_int k, num_of_int (2 * den), zero, num_of_int den) in
-
- let rec digit k z n row col acc =
- if n > 0 then
- let y = next z in
- if safe z y then
- if col = 10 then (
- let row = row + 10 in
- digit k (prod z y) (n - 1) row 1
- ((Printf.sprintf "\t:%i\n%s" row (string_of_num y)) :: acc)
- ) else (
- digit k (prod z y) (n - 1) row (col + 1)
- ((string_of_num y) :: acc)
- )
- else digit (k + 1) (cons z k) n row col acc
- else
- (Printf.sprintf "%*s\t:%i\n" (10 - col) "" (row + col)) :: acc in
-
- let digits n = digit 1 unit n 0 0 [] in
- eq_string (String.concat "" (List.rev (digits 1000)),
- IntegersTest.pi_1000_digits);
+let bigtests = "NumbersBig" >::: [
+ "pi 1000 digits" >::
+ (fun () ->
+ eq_string (pi_digits 1000, IntegersTest.pi_1000_digits);
);
]
Modified: trunk/Toss/Solver/Num/Rationals.ml
===================================================================
--- trunk/Toss/Solver/Num/Rationals.ml 2012-02-17 02:36:56 UTC (rev 1676)
+++ trunk/Toss/Solver/Num/Rationals.ml 2012-02-18 02:07:02 UTC (rev 1677)
@@ -401,7 +401,7 @@
s1 contains one more digit than desired for the round off operation *)
if n >= 0 then begin
let s1 =
- Naturals.N.string_of_nat
+ Naturals.string_of_nat
(nat_of_big_int
(div_big_int
(base_power_big_int
@@ -478,7 +478,7 @@
div_big_int (base_power_big_int
10 k (abs_big_int r.numerator))
r.denominator) in
- Naturals.N.string_of_nat nat) in
+ Naturals.string_of_nat nat) in
if (round_futur_last_digit s 0 (String.length s))
then
let m = num_decimal_digits_int (succ msd) in
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
