cannam@167
|
1 (*
|
cannam@167
|
2 * Copyright (c) 1997-1999 Massachusetts Institute of Technology
|
cannam@167
|
3 * Copyright (c) 2003, 2007-14 Matteo Frigo
|
cannam@167
|
4 * Copyright (c) 2003, 2007-14 Massachusetts Institute of Technology
|
cannam@167
|
5 *
|
cannam@167
|
6 * This program is free software; you can redistribute it and/or modify
|
cannam@167
|
7 * it under the terms of the GNU General Public License as published by
|
cannam@167
|
8 * the Free Software Foundation; either version 2 of the License, or
|
cannam@167
|
9 * (at your option) any later version.
|
cannam@167
|
10 *
|
cannam@167
|
11 * This program is distributed in the hope that it will be useful,
|
cannam@167
|
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
|
cannam@167
|
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
cannam@167
|
14 * GNU General Public License for more details.
|
cannam@167
|
15 *
|
cannam@167
|
16 * You should have received a copy of the GNU General Public License
|
cannam@167
|
17 * along with this program; if not, write to the Free Software
|
cannam@167
|
18 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
|
cannam@167
|
19 *
|
cannam@167
|
20 *)
|
cannam@167
|
21
|
cannam@167
|
22 (* various utility functions *)
|
cannam@167
|
23 open List
|
cannam@167
|
24 open Unix
|
cannam@167
|
25
|
cannam@167
|
26 (*****************************************
|
cannam@167
|
27 * Integer operations
|
cannam@167
|
28 *****************************************)
|
cannam@167
|
29 (* fint the inverse of n modulo m *)
|
cannam@167
|
30 let invmod n m =
|
cannam@167
|
31 let rec loop i =
|
cannam@167
|
32 if ((i * n) mod m == 1) then i
|
cannam@167
|
33 else loop (i + 1)
|
cannam@167
|
34 in
|
cannam@167
|
35 loop 1
|
cannam@167
|
36
|
cannam@167
|
37 (* Yooklid's algorithm *)
|
cannam@167
|
38 let rec gcd n m =
|
cannam@167
|
39 if (n > m)
|
cannam@167
|
40 then gcd m n
|
cannam@167
|
41 else
|
cannam@167
|
42 let r = m mod n
|
cannam@167
|
43 in
|
cannam@167
|
44 if (r == 0) then n
|
cannam@167
|
45 else gcd r n
|
cannam@167
|
46
|
cannam@167
|
47 (* reduce the fraction m/n to lowest terms, modulo factors of n/n *)
|
cannam@167
|
48 let lowest_terms n m =
|
cannam@167
|
49 if (m mod n == 0) then
|
cannam@167
|
50 (1,0)
|
cannam@167
|
51 else
|
cannam@167
|
52 let nn = (abs n) in let mm = m * (n / nn)
|
cannam@167
|
53 in let mpos =
|
cannam@167
|
54 if (mm > 0) then (mm mod nn)
|
cannam@167
|
55 else (mm + (1 + (abs mm) / nn) * nn) mod nn
|
cannam@167
|
56 and d = gcd nn (abs mm)
|
cannam@167
|
57 in (nn / d, mpos / d)
|
cannam@167
|
58
|
cannam@167
|
59 (* find a generator for the multiplicative group mod p
|
cannam@167
|
60 (where p must be prime for a generator to exist!!) *)
|
cannam@167
|
61
|
cannam@167
|
62 exception No_Generator
|
cannam@167
|
63
|
cannam@167
|
64 let find_generator p =
|
cannam@167
|
65 let rec period x prod =
|
cannam@167
|
66 if (prod == 1) then 1
|
cannam@167
|
67 else 1 + (period x (prod * x mod p))
|
cannam@167
|
68 in let rec findgen x =
|
cannam@167
|
69 if (x == 0) then raise No_Generator
|
cannam@167
|
70 else if ((period x x) == (p - 1)) then x
|
cannam@167
|
71 else findgen ((x + 1) mod p)
|
cannam@167
|
72 in findgen 1
|
cannam@167
|
73
|
cannam@167
|
74 (* raise x to a power n modulo p (requires n > 0) (in principle,
|
cannam@167
|
75 negative powers would be fine, provided that x and p are relatively
|
cannam@167
|
76 prime...we don't need this functionality, though) *)
|
cannam@167
|
77
|
cannam@167
|
78 exception Negative_Power
|
cannam@167
|
79
|
cannam@167
|
80 let rec pow_mod x n p =
|
cannam@167
|
81 if (n == 0) then 1
|
cannam@167
|
82 else if (n < 0) then raise Negative_Power
|
cannam@167
|
83 else if (n mod 2 == 0) then pow_mod (x * x mod p) (n / 2) p
|
cannam@167
|
84 else x * (pow_mod x (n - 1) p) mod p
|
cannam@167
|
85
|
cannam@167
|
86 (******************************************
|
cannam@167
|
87 * auxiliary functions
|
cannam@167
|
88 ******************************************)
|
cannam@167
|
89 let rec forall id combiner a b f =
|
cannam@167
|
90 if (a >= b) then id
|
cannam@167
|
91 else combiner (f a) (forall id combiner (a + 1) b f)
|
cannam@167
|
92
|
cannam@167
|
93 let sum_list l = fold_right (+) l 0
|
cannam@167
|
94 let max_list l = fold_right (max) l (-999999)
|
cannam@167
|
95 let min_list l = fold_right (min) l 999999
|
cannam@167
|
96 let count pred = fold_left
|
cannam@167
|
97 (fun a elem -> if (pred elem) then 1 + a else a) 0
|
cannam@167
|
98 let remove elem = List.filter (fun e -> (e != elem))
|
cannam@167
|
99 let cons a b = a :: b
|
cannam@167
|
100 let null = function
|
cannam@167
|
101 [] -> true
|
cannam@167
|
102 | _ -> false
|
cannam@167
|
103 let for_list l f = List.iter f l
|
cannam@167
|
104 let rmap l f = List.map f l
|
cannam@167
|
105
|
cannam@167
|
106 (* functional composition *)
|
cannam@167
|
107 let (@@) f g x = f (g x)
|
cannam@167
|
108
|
cannam@167
|
109 let forall_flat a b = forall [] (@) a b
|
cannam@167
|
110
|
cannam@167
|
111 let identity x = x
|
cannam@167
|
112
|
cannam@167
|
113 let rec minimize f = function
|
cannam@167
|
114 [] -> None
|
cannam@167
|
115 | elem :: rest ->
|
cannam@167
|
116 match minimize f rest with
|
cannam@167
|
117 None -> Some elem
|
cannam@167
|
118 | Some x -> if (f x) >= (f elem) then Some elem else Some x
|
cannam@167
|
119
|
cannam@167
|
120
|
cannam@167
|
121 let rec find_elem condition = function
|
cannam@167
|
122 [] -> None
|
cannam@167
|
123 | elem :: rest ->
|
cannam@167
|
124 if condition elem then
|
cannam@167
|
125 Some elem
|
cannam@167
|
126 else
|
cannam@167
|
127 find_elem condition rest
|
cannam@167
|
128
|
cannam@167
|
129
|
cannam@167
|
130 (* find x, x >= a, such that (p x) is true *)
|
cannam@167
|
131 let rec suchthat a pred =
|
cannam@167
|
132 if (pred a) then a else suchthat (a + 1) pred
|
cannam@167
|
133
|
cannam@167
|
134 (* print an information message *)
|
cannam@167
|
135 let info string =
|
cannam@167
|
136 if !Magic.verbose then begin
|
cannam@167
|
137 let now = Unix.times ()
|
cannam@167
|
138 and pid = Unix.getpid () in
|
cannam@167
|
139 prerr_string ((string_of_int pid) ^ ": " ^
|
cannam@167
|
140 "at t = " ^ (string_of_float now.tms_utime) ^ " : ");
|
cannam@167
|
141 prerr_string (string ^ "\n");
|
cannam@167
|
142 flush Pervasives.stderr;
|
cannam@167
|
143 end
|
cannam@167
|
144
|
cannam@167
|
145 (* iota n produces the list [0; 1; ...; n - 1] *)
|
cannam@167
|
146 let iota n = forall [] cons 0 n identity
|
cannam@167
|
147
|
cannam@167
|
148 (* interval a b produces the list [a; 1; ...; b - 1] *)
|
cannam@167
|
149 let interval a b = List.map ((+) a) (iota (b - a))
|
cannam@167
|
150
|
cannam@167
|
151 (*
|
cannam@167
|
152 * freeze a function, i.e., compute it only once on demand, and
|
cannam@167
|
153 * cache it into an array.
|
cannam@167
|
154 *)
|
cannam@167
|
155 let array n f =
|
cannam@167
|
156 let a = Array.init n (fun i -> lazy (f i))
|
cannam@167
|
157 in fun i -> Lazy.force a.(i)
|
cannam@167
|
158
|
cannam@167
|
159
|
cannam@167
|
160 let rec take n l =
|
cannam@167
|
161 match (n, l) with
|
cannam@167
|
162 (0, _) -> []
|
cannam@167
|
163 | (n, (a :: b)) -> a :: (take (n - 1) b)
|
cannam@167
|
164 | _ -> failwith "take"
|
cannam@167
|
165
|
cannam@167
|
166 let rec drop n l =
|
cannam@167
|
167 match (n, l) with
|
cannam@167
|
168 (0, _) -> l
|
cannam@167
|
169 | (n, (_ :: b)) -> drop (n - 1) b
|
cannam@167
|
170 | _ -> failwith "drop"
|
cannam@167
|
171
|
cannam@167
|
172
|
cannam@167
|
173 let either a b =
|
cannam@167
|
174 match a with
|
cannam@167
|
175 Some x -> x
|
cannam@167
|
176 | _ -> b
|