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