Chris@82: (* Chris@82: * Copyright (c) 1997-1999 Massachusetts Institute of Technology Chris@82: * Copyright (c) 2003, 2007-14 Matteo Frigo Chris@82: * Copyright (c) 2003, 2007-14 Massachusetts Institute of Technology Chris@82: * Chris@82: * This program is free software; you can redistribute it and/or modify Chris@82: * it under the terms of the GNU General Public License as published by Chris@82: * the Free Software Foundation; either version 2 of the License, or Chris@82: * (at your option) any later version. Chris@82: * Chris@82: * This program is distributed in the hope that it will be useful, Chris@82: * but WITHOUT ANY WARRANTY; without even the implied warranty of Chris@82: * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the Chris@82: * GNU General Public License for more details. Chris@82: * Chris@82: * You should have received a copy of the GNU General Public License Chris@82: * along with this program; if not, write to the Free Software Chris@82: * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA Chris@82: * Chris@82: *) Chris@82: Chris@82: (* various utility functions *) Chris@82: open List Chris@82: open Unix Chris@82: Chris@82: (***************************************** Chris@82: * Integer operations Chris@82: *****************************************) Chris@82: (* fint the inverse of n modulo m *) Chris@82: let invmod n m = Chris@82: let rec loop i = Chris@82: if ((i * n) mod m == 1) then i Chris@82: else loop (i + 1) Chris@82: in Chris@82: loop 1 Chris@82: Chris@82: (* Yooklid's algorithm *) Chris@82: let rec gcd n m = Chris@82: if (n > m) Chris@82: then gcd m n Chris@82: else Chris@82: let r = m mod n Chris@82: in Chris@82: if (r == 0) then n Chris@82: else gcd r n Chris@82: Chris@82: (* reduce the fraction m/n to lowest terms, modulo factors of n/n *) Chris@82: let lowest_terms n m = Chris@82: if (m mod n == 0) then Chris@82: (1,0) Chris@82: else Chris@82: let nn = (abs n) in let mm = m * (n / nn) Chris@82: in let mpos = Chris@82: if (mm > 0) then (mm mod nn) Chris@82: else (mm + (1 + (abs mm) / nn) * nn) mod nn Chris@82: and d = gcd nn (abs mm) Chris@82: in (nn / d, mpos / d) Chris@82: Chris@82: (* find a generator for the multiplicative group mod p Chris@82: (where p must be prime for a generator to exist!!) *) Chris@82: Chris@82: exception No_Generator Chris@82: Chris@82: let find_generator p = Chris@82: let rec period x prod = Chris@82: if (prod == 1) then 1 Chris@82: else 1 + (period x (prod * x mod p)) Chris@82: in let rec findgen x = Chris@82: if (x == 0) then raise No_Generator Chris@82: else if ((period x x) == (p - 1)) then x Chris@82: else findgen ((x + 1) mod p) Chris@82: in findgen 1 Chris@82: Chris@82: (* raise x to a power n modulo p (requires n > 0) (in principle, Chris@82: negative powers would be fine, provided that x and p are relatively Chris@82: prime...we don't need this functionality, though) *) Chris@82: Chris@82: exception Negative_Power Chris@82: Chris@82: let rec pow_mod x n p = Chris@82: if (n == 0) then 1 Chris@82: else if (n < 0) then raise Negative_Power Chris@82: else if (n mod 2 == 0) then pow_mod (x * x mod p) (n / 2) p Chris@82: else x * (pow_mod x (n - 1) p) mod p Chris@82: Chris@82: (****************************************** Chris@82: * auxiliary functions Chris@82: ******************************************) Chris@82: let rec forall id combiner a b f = Chris@82: if (a >= b) then id Chris@82: else combiner (f a) (forall id combiner (a + 1) b f) Chris@82: Chris@82: let sum_list l = fold_right (+) l 0 Chris@82: let max_list l = fold_right (max) l (-999999) Chris@82: let min_list l = fold_right (min) l 999999 Chris@82: let count pred = fold_left Chris@82: (fun a elem -> if (pred elem) then 1 + a else a) 0 Chris@82: let remove elem = List.filter (fun e -> (e != elem)) Chris@82: let cons a b = a :: b Chris@82: let null = function Chris@82: [] -> true Chris@82: | _ -> false Chris@82: let for_list l f = List.iter f l Chris@82: let rmap l f = List.map f l Chris@82: Chris@82: (* functional composition *) Chris@82: let (@@) f g x = f (g x) Chris@82: Chris@82: let forall_flat a b = forall [] (@) a b Chris@82: Chris@82: let identity x = x Chris@82: Chris@82: let rec minimize f = function Chris@82: [] -> None Chris@82: | elem :: rest -> Chris@82: match minimize f rest with Chris@82: None -> Some elem Chris@82: | Some x -> if (f x) >= (f elem) then Some elem else Some x Chris@82: Chris@82: Chris@82: let rec find_elem condition = function Chris@82: [] -> None Chris@82: | elem :: rest -> Chris@82: if condition elem then Chris@82: Some elem Chris@82: else Chris@82: find_elem condition rest Chris@82: Chris@82: Chris@82: (* find x, x >= a, such that (p x) is true *) Chris@82: let rec suchthat a pred = Chris@82: if (pred a) then a else suchthat (a + 1) pred Chris@82: Chris@82: (* print an information message *) Chris@82: let info string = Chris@82: if !Magic.verbose then begin Chris@82: let now = Unix.times () Chris@82: and pid = Unix.getpid () in Chris@82: prerr_string ((string_of_int pid) ^ ": " ^ Chris@82: "at t = " ^ (string_of_float now.tms_utime) ^ " : "); Chris@82: prerr_string (string ^ "\n"); Chris@82: flush Pervasives.stderr; Chris@82: end Chris@82: Chris@82: (* iota n produces the list [0; 1; ...; n - 1] *) Chris@82: let iota n = forall [] cons 0 n identity Chris@82: Chris@82: (* interval a b produces the list [a; 1; ...; b - 1] *) Chris@82: let interval a b = List.map ((+) a) (iota (b - a)) Chris@82: Chris@82: (* Chris@82: * freeze a function, i.e., compute it only once on demand, and Chris@82: * cache it into an array. Chris@82: *) Chris@82: let array n f = Chris@82: let a = Array.init n (fun i -> lazy (f i)) Chris@82: in fun i -> Lazy.force a.(i) Chris@82: Chris@82: Chris@82: let rec take n l = Chris@82: match (n, l) with Chris@82: (0, _) -> [] Chris@82: | (n, (a :: b)) -> a :: (take (n - 1) b) Chris@82: | _ -> failwith "take" Chris@82: Chris@82: let rec drop n l = Chris@82: match (n, l) with Chris@82: (0, _) -> l Chris@82: | (n, (_ :: b)) -> drop (n - 1) b Chris@82: | _ -> failwith "drop" Chris@82: Chris@82: Chris@82: let either a b = Chris@82: match a with Chris@82: Some x -> x Chris@82: | _ -> b