annotate src/fftw-3.3.3/genfft/util.ml @ 23:619f715526df sv_v2.1

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