PageRenderTime 61ms CodeModel.GetById 30ms RepoModel.GetById 1ms app.codeStats 0ms

/coq-8.3pl2/lib/bigint.ml

https://bitbucket.org/yoshihiro503/coq2sml
OCaml | 411 lines | 303 code | 56 blank | 52 comment | 128 complexity | 2e1c8436d59f5730ea391dcfd94633c5 MD5 | raw file
Possible License(s): LGPL-2.1
  1. (************************************************************************)
  2. (* v * The Coq Proof Assistant / The Coq Development Team *)
  3. (* <O___,, * INRIA - CNRS - LIX - LRI - PPS - Copyright 1999-2010 *)
  4. (* \VV/ **************************************************************)
  5. (* // * This file is distributed under the terms of the *)
  6. (* * GNU Lesser General Public License Version 2.1 *)
  7. (************************************************************************)
  8. (* $Id: bigint.ml 13323 2010-07-24 15:57:30Z herbelin $ *)
  9. (*i*)
  10. open Pp
  11. (*i*)
  12. (***************************************************)
  13. (* Basic operations on (unbounded) integer numbers *)
  14. (***************************************************)
  15. (* An integer is canonically represented as an array of k-digits blocs.
  16. 0 is represented by the empty array and -1 by the singleton [|-1|].
  17. The first bloc is in the range ]0;10^k[ for positive numbers.
  18. The first bloc is in the range ]-10^k;-1[ for negative ones.
  19. All other blocs are numbers in the range [0;10^k[.
  20. Negative numbers are represented using 2's complementation. For instance,
  21. with 4-digits blocs, [-9655;6789] denotes -96543211
  22. *)
  23. (* The base is a power of 10 in order to facilitate the parsing and printing
  24. of numbers in digital notation.
  25. All functions, to the exception of to_string and of_string should work
  26. with an arbitrary base, even if not a power of 10.
  27. In practice, we set k=4 so that no overflow in ocaml machine words
  28. (i.e. the interval [-2^30;2^30-1]) occur when multiplying two
  29. numbers less than (10^k)
  30. *)
  31. (* The main parameters *)
  32. let size =
  33. let rec log10 n = if n < 10 then 0 else 1 + log10 (n / 10) in
  34. (log10 max_int) / 2
  35. let format_size =
  36. (* How to parametrize a printf format *)
  37. if size = 4 then Printf.sprintf "%04d"
  38. else fun n ->
  39. let rec aux j l n =
  40. if j=size then l else aux (j+1) (string_of_int (n mod 10) :: l) (n/10)
  41. in String.concat "" (aux 0 [] n)
  42. (* The base is 10^size *)
  43. let base =
  44. let rec exp10 = function 0 -> 1 | n -> 10 * exp10 (n-1) in exp10 size
  45. (* Basic numbers *)
  46. let zero = [||]
  47. let neg_one = [|-1|]
  48. (* Sign of an integer *)
  49. let is_strictly_neg n = n<>[||] && n.(0) < 0
  50. let is_strictly_pos n = n<>[||] && n.(0) > 0
  51. let is_neg_or_zero n = n=[||] or n.(0) < 0
  52. let is_pos_or_zero n = n=[||] or n.(0) > 0
  53. let normalize_pos n =
  54. let k = ref 0 in
  55. while !k < Array.length n & n.(!k) = 0 do incr k done;
  56. Array.sub n !k (Array.length n - !k)
  57. let normalize_neg n =
  58. let k = ref 1 in
  59. while !k < Array.length n & n.(!k) = base - 1 do incr k done;
  60. let n' = Array.sub n !k (Array.length n - !k) in
  61. if Array.length n' = 0 then [|-1|] else (n'.(0) <- n'.(0) - base; n')
  62. let rec normalize n =
  63. if Array.length n = 0 then n else
  64. if n.(0) = -1 then normalize_neg n else normalize_pos n
  65. let neg m =
  66. if m = zero then zero else
  67. let n = Array.copy m in
  68. let i = ref (Array.length m - 1) in
  69. while !i > 0 & n.(!i) = 0 do decr i done;
  70. if !i > 0 then begin
  71. n.(!i) <- base - n.(!i); decr i;
  72. while !i > 0 do n.(!i) <- base - 1 - n.(!i); decr i done;
  73. n.(0) <- - n.(0) - 1;
  74. if n.(0) < -1 then (n.(0) <- n.(0) + base; Array.append [| -1 |] n) else
  75. if n.(0) = - base then (n.(0) <- 0; Array.append [| -1 |] n)
  76. else normalize n
  77. end else (n.(0) <- - n.(0); n)
  78. let push_carry r j =
  79. let j = ref j in
  80. while !j > 0 & r.(!j) < 0 do
  81. r.(!j) <- r.(!j) + base; decr j; r.(!j) <- r.(!j) - 1
  82. done;
  83. while !j > 0 & r.(!j) >= base do
  84. r.(!j) <- r.(!j) - base; decr j; r.(!j) <- r.(!j) + 1
  85. done;
  86. if r.(0) >= base then (r.(0) <- r.(0) - base; Array.append [| 1 |] r)
  87. else if r.(0) < -base then (r.(0) <- r.(0) + 2*base; Array.append [| -2 |] r)
  88. else if r.(0) = -base then (r.(0) <- 0; Array.append [| -1 |] r)
  89. else normalize r
  90. let add_to r a j =
  91. if a = zero then r else begin
  92. for i = Array.length r - 1 downto j+1 do
  93. r.(i) <- r.(i) + a.(i-j);
  94. if r.(i) >= base then (r.(i) <- r.(i) - base; r.(i-1) <- r.(i-1) + 1)
  95. done;
  96. r.(j) <- r.(j) + a.(0);
  97. push_carry r j
  98. end
  99. let add n m =
  100. let d = Array.length n - Array.length m in
  101. if d > 0 then add_to (Array.copy n) m d else add_to (Array.copy m) n (-d)
  102. let sub_to r a j =
  103. if a = zero then r else begin
  104. for i = Array.length r - 1 downto j+1 do
  105. r.(i) <- r.(i) - a.(i-j);
  106. if r.(i) < 0 then (r.(i) <- r.(i) + base; r.(i-1) <- r.(i-1) - 1)
  107. done;
  108. r.(j) <- r.(j) - a.(0);
  109. push_carry r j
  110. end
  111. let sub n m =
  112. let d = Array.length n - Array.length m in
  113. if d >= 0 then sub_to (Array.copy n) m d
  114. else let r = neg m in add_to r n (Array.length r - Array.length n)
  115. let rec mult m n =
  116. if m = zero or n = zero then zero else
  117. let l = Array.length m + Array.length n in
  118. let r = Array.create l 0 in
  119. for i = Array.length m - 1 downto 0 do
  120. for j = Array.length n - 1 downto 0 do
  121. let p = m.(i) * n.(j) + r.(i+j+1) in
  122. let (q,s) =
  123. if p < 0
  124. then (p + 1) / base - 1, (p + 1) mod base + base - 1
  125. else p / base, p mod base in
  126. r.(i+j+1) <- s;
  127. if q <> 0 then r.(i+j) <- r.(i+j) + q;
  128. done
  129. done;
  130. normalize r
  131. let rec less_than_same_size m n i j =
  132. i < Array.length m &&
  133. (m.(i) < n.(j) or (m.(i) = n.(j) && less_than_same_size m n (i+1) (j+1)))
  134. let less_than m n =
  135. if is_strictly_neg m then
  136. is_pos_or_zero n or Array.length m > Array.length n
  137. or (Array.length m = Array.length n && less_than_same_size m n 0 0)
  138. else
  139. is_strictly_pos n && (Array.length m < Array.length n or
  140. (Array.length m = Array.length n && less_than_same_size m n 0 0))
  141. let equal m n = (m = n)
  142. let less_than_shift_pos k m n =
  143. (Array.length m - k < Array.length n)
  144. or (Array.length m - k = Array.length n && less_than_same_size m n k 0)
  145. let rec can_divide k m d i =
  146. (i = Array.length d) or
  147. (m.(k+i) > d.(i)) or
  148. (m.(k+i) = d.(i) && can_divide k m d (i+1))
  149. (* computes m - d * q * base^(|m|-k) in-place on positive numbers *)
  150. let sub_mult m d q k =
  151. if q <> 0 then
  152. for i = Array.length d - 1 downto 0 do
  153. let v = d.(i) * q in
  154. m.(k+i) <- m.(k+i) - v mod base;
  155. if m.(k+i) < 0 then (m.(k+i) <- m.(k+i) + base; m.(k+i-1) <- m.(k+i-1) -1);
  156. if v >= base then m.(k+i-1) <- m.(k+i-1) - v / base;
  157. done
  158. let euclid m d =
  159. let isnegm, m =
  160. if is_strictly_neg m then (-1),neg m else 1,Array.copy m in
  161. let isnegd, d = if is_strictly_neg d then (-1),neg d else 1,d in
  162. if d = zero then raise Division_by_zero;
  163. let q,r =
  164. if less_than m d then (zero,m) else
  165. let ql = Array.length m - Array.length d in
  166. let q = Array.create (ql+1) 0 in
  167. let i = ref 0 in
  168. while not (less_than_shift_pos !i m d) do
  169. if m.(!i)=0 then incr i else
  170. if can_divide !i m d 0 then begin
  171. let v =
  172. if Array.length d > 1 && d.(0) <> m.(!i) then
  173. (m.(!i) * base + m.(!i+1)) / (d.(0) * base + d.(1) + 1)
  174. else
  175. m.(!i) / d.(0) in
  176. q.(!i) <- q.(!i) + v;
  177. sub_mult m d v !i
  178. end else begin
  179. let v = (m.(!i) * base + m.(!i+1)) / (d.(0) + 1) in
  180. q.(!i) <- q.(!i) + v / base;
  181. sub_mult m d (v / base) !i;
  182. q.(!i+1) <- q.(!i+1) + v mod base;
  183. if q.(!i+1) >= base then
  184. (q.(!i+1) <- q.(!i+1)-base; q.(!i) <- q.(!i)+1);
  185. sub_mult m d (v mod base) (!i+1)
  186. end
  187. done;
  188. (normalize q, normalize m) in
  189. (if isnegd * isnegm = -1 then neg q else q),
  190. (if isnegm = -1 then neg r else r)
  191. (* Parsing/printing ordinary 10-based numbers *)
  192. let of_string s =
  193. let isneg = String.length s > 1 & s.[0] = '-' in
  194. let n = if isneg then 1 else 0 in
  195. let d = ref n in
  196. while !d < String.length s && s.[!d] = '0' do incr d done;
  197. if !d = String.length s then zero else
  198. let r = (String.length s - !d) mod size in
  199. let h = String.sub s (!d) r in
  200. if !d = String.length s - 1 && isneg && h="1" then neg_one else
  201. let e = if h<>"" then 1 else 0 in
  202. let l = (String.length s - !d) / size in
  203. let a = Array.create (l + e + n) 0 in
  204. if isneg then begin
  205. a.(0) <- (-1);
  206. let carry = ref 0 in
  207. for i=l downto 1 do
  208. let v = int_of_string (String.sub s ((i-1)*size + !d +r) size)+ !carry in
  209. if v <> 0 then (a.(i+e)<- base - v; carry := 1) else carry := 0
  210. done;
  211. if e=1 then a.(1) <- base - !carry - int_of_string h;
  212. end
  213. else begin
  214. if e=1 then a.(0) <- int_of_string h;
  215. for i=1 to l do
  216. a.(i+e-1) <- int_of_string (String.sub s ((i-1)*size + !d + r) size)
  217. done
  218. end;
  219. a
  220. let to_string_pos sgn n =
  221. if Array.length n = 0 then "0" else
  222. sgn ^
  223. String.concat ""
  224. (string_of_int n.(0) :: List.map format_size (List.tl (Array.to_list n)))
  225. let to_string n =
  226. if is_strictly_neg n then to_string_pos "-" (neg n)
  227. else to_string_pos "" n
  228. (******************************************************************)
  229. (* Optimized operations on (unbounded) integer numbers *)
  230. (* integers smaller than base are represented as machine integers *)
  231. (******************************************************************)
  232. type bigint = Obj.t
  233. let ints_of_int n =
  234. if n >= base then [| n / base; n mod base |]
  235. else if n <= - base then [| n / base - 1; n mod base + base |]
  236. else if n = 0 then [| |] else [| n |]
  237. let big_of_int n =
  238. if n >= base then Obj.repr [| n / base; n mod base |]
  239. else if n <= - base then Obj.repr [| n / base - 1; n mod base + base |]
  240. else Obj.repr n
  241. let big_of_ints n =
  242. let n = normalize n in
  243. if n = zero then Obj.repr 0 else
  244. if Array.length n = 1 then Obj.repr n.(0) else
  245. Obj.repr n
  246. let coerce_to_int = (Obj.magic : Obj.t -> int)
  247. let coerce_to_ints = (Obj.magic : Obj.t -> int array)
  248. let ints_of_z n =
  249. if Obj.is_int n then ints_of_int (coerce_to_int n)
  250. else coerce_to_ints n
  251. let app_pair f (m, n) =
  252. (f m, f n)
  253. let add m n =
  254. if Obj.is_int m & Obj.is_int n
  255. then big_of_int (coerce_to_int m + coerce_to_int n)
  256. else big_of_ints (add (ints_of_z m) (ints_of_z n))
  257. let sub m n =
  258. if Obj.is_int m & Obj.is_int n
  259. then big_of_int (coerce_to_int m - coerce_to_int n)
  260. else big_of_ints (sub (ints_of_z m) (ints_of_z n))
  261. let mult m n =
  262. if Obj.is_int m & Obj.is_int n
  263. then big_of_int (coerce_to_int m * coerce_to_int n)
  264. else big_of_ints (mult (ints_of_z m) (ints_of_z n))
  265. let euclid m n =
  266. if Obj.is_int m & Obj.is_int n
  267. then app_pair big_of_int
  268. (coerce_to_int m / coerce_to_int n, coerce_to_int m mod coerce_to_int n)
  269. else app_pair big_of_ints (euclid (ints_of_z m) (ints_of_z n))
  270. let less_than m n =
  271. if Obj.is_int m & Obj.is_int n
  272. then coerce_to_int m < coerce_to_int n
  273. else less_than (ints_of_z m) (ints_of_z n)
  274. let neg n =
  275. if Obj.is_int n then big_of_int (- (coerce_to_int n))
  276. else big_of_ints (neg (ints_of_z n))
  277. let of_string m = big_of_ints (of_string m)
  278. let to_string m = to_string (ints_of_z m)
  279. let zero = big_of_int 0
  280. let one = big_of_int 1
  281. let sub_1 n = sub n one
  282. let add_1 n = add n one
  283. let two = big_of_int 2
  284. let mult_2 n = add n n
  285. let div2_with_rest n =
  286. let (q,b) = euclid n two in
  287. (q, b = one)
  288. let is_strictly_neg n = is_strictly_neg (ints_of_z n)
  289. let is_strictly_pos n = is_strictly_pos (ints_of_z n)
  290. let is_neg_or_zero n = is_neg_or_zero (ints_of_z n)
  291. let is_pos_or_zero n = is_pos_or_zero (ints_of_z n)
  292. let pr_bigint n = str (to_string n)
  293. (* spiwack: computes n^m *)
  294. (* The basic idea of the algorithm is that n^(2m) = (n^2)^m *)
  295. (* In practice the algorithm performs :
  296. k*n^0 = k
  297. k*n^(2m) = k*(n*n)^m
  298. k*n^(2m+1) = (n*k)*(n*n)^m *)
  299. let pow =
  300. let rec pow_aux odd_rest n m = (* odd_rest is the k from above *)
  301. if is_neg_or_zero m then
  302. odd_rest
  303. else
  304. let (quo,rem) = div2_with_rest m in
  305. pow_aux
  306. ((* [if m mod 2 = 1]*)
  307. if rem then
  308. mult n odd_rest
  309. else
  310. odd_rest )
  311. (* quo = [m/2] *)
  312. (mult n n) quo
  313. in
  314. pow_aux one
  315. (* Testing suite *)
  316. let check () =
  317. let numbers = [
  318. "1";"2";"99";"100";"101";"9999";"10000";"10001";
  319. "999999";"1000000";"1000001";"99999999";"100000000";"100000001";
  320. "1234";"5678";"12345678";"987654321";
  321. "-1";"-2";"-99";"-100";"-101";"-9999";"-10000";"-10001";
  322. "-999999";"-1000000";"-1000001";"-99999999";"-100000000";"-100000001";
  323. "-1234";"-5678";"-12345678";"-987654321";"0"
  324. ]
  325. in
  326. let eucl n m =
  327. let n' = abs_float n and m' = abs_float m in
  328. let q' = floor (n' /. m') in let r' = n' -. m' *. q' in
  329. (if n *. m < 0. & q' <> 0. then -. q' else q'),
  330. (if n < 0. then -. r' else r') in
  331. let round f = floor (abs_float f +. 0.5) *. (if f < 0. then -1. else 1.) in
  332. let i = ref 0 in
  333. let compare op n n' =
  334. incr i;
  335. let s = Printf.sprintf "%30s" (to_string n) in
  336. let s' = Printf.sprintf "% 30.0f" (round n') in
  337. if s <> s' then Printf.printf "%s: %s <> %s\n" op s s' in
  338. List.iter (fun a -> List.iter (fun b ->
  339. let n = of_string a and m = of_string b in
  340. let n' = float_of_string a and m' = float_of_string b in
  341. let a = add n m and a' = n' +. m' in
  342. let s = sub n m and s' = n' -. m' in
  343. let p = mult n m and p' = n' *. m' in
  344. let q,r = try euclid n m with Division_by_zero -> zero,zero
  345. and q',r' = eucl n' m' in
  346. compare "+" a a';
  347. compare "-" s s';
  348. compare "*" p p';
  349. compare "/" q q';
  350. compare "%" r r') numbers) numbers;
  351. Printf.printf "%i tests done\n" !i