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