Chris@10: (* Chris@10: * Copyright (c) 1997-1999 Massachusetts Institute of Technology Chris@10: * Copyright (c) 2003, 2007-11 Matteo Frigo Chris@10: * Copyright (c) 2003, 2007-11 Massachusetts Institute of Technology Chris@10: * Chris@10: * This program is free software; you can redistribute it and/or modify Chris@10: * it under the terms of the GNU General Public License as published by Chris@10: * the Free Software Foundation; either version 2 of the License, or Chris@10: * (at your option) any later version. Chris@10: * Chris@10: * This program is distributed in the hope that it will be useful, Chris@10: * but WITHOUT ANY WARRANTY; without even the implied warranty of Chris@10: * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the Chris@10: * GNU General Public License for more details. Chris@10: * Chris@10: * You should have received a copy of the GNU General Public License Chris@10: * along with this program; if not, write to the Free Software Chris@10: * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA Chris@10: * Chris@10: *) Chris@10: Chris@10: Chris@10: (* This is the part of the generator that actually computes the FFT Chris@10: in symbolic form *) Chris@10: Chris@10: open Complex Chris@10: open Util Chris@10: Chris@10: (* choose a suitable factor of n *) Chris@10: let choose_factor n = Chris@10: (* first choice: i such that gcd(i, n / i) = 1, i as big as possible *) Chris@10: let choose1 n = Chris@10: let rec loop i f = Chris@10: if (i * i > n) then f Chris@10: else if ((n mod i) == 0 && gcd i (n / i) == 1) then loop (i + 1) i Chris@10: else loop (i + 1) f Chris@10: in loop 1 1 Chris@10: Chris@10: (* second choice: the biggest factor i of n, where i < sqrt(n), if any *) Chris@10: and choose2 n = Chris@10: let rec loop i f = Chris@10: if (i * i > n) then f Chris@10: else if ((n mod i) == 0) then loop (i + 1) i Chris@10: else loop (i + 1) f Chris@10: in loop 1 1 Chris@10: Chris@10: in let i = choose1 n in Chris@10: if (i > 1) then i Chris@10: else choose2 n Chris@10: Chris@10: let is_power_of_two n = (n > 0) && ((n - 1) land n == 0) Chris@10: Chris@10: let rec dft_prime sign n input = Chris@10: let sum filter i = Chris@10: sigma 0 n (fun j -> Chris@10: let coeff = filter (exp n (sign * i * j)) Chris@10: in coeff @* (input j)) in Chris@10: let computation_even = array n (sum identity) Chris@10: and computation_odd = Chris@10: let sumr = array n (sum real) Chris@10: and sumi = array n (sum ((times Complex.i) @@ imag)) in Chris@10: array n (fun i -> Chris@10: if (i = 0) then Chris@10: (* expose some common subexpressions *) Chris@10: input 0 @+ Chris@10: sigma 1 ((n + 1) / 2) (fun j -> input j @+ input (n - j)) Chris@10: else Chris@10: let i' = min i (n - i) in Chris@10: if (i < n - i) then Chris@10: sumr i' @+ sumi i' Chris@10: else Chris@10: sumr i' @- sumi i') in Chris@10: if (n >= !Magic.rader_min) then Chris@10: dft_rader sign n input Chris@10: else if (n == 2) then Chris@10: computation_even Chris@10: else Chris@10: computation_odd Chris@10: Chris@10: Chris@10: and dft_rader sign p input = Chris@10: let half = Chris@10: let one_half = inverse_int 2 in Chris@10: times one_half Chris@10: Chris@10: and make_product n a b = Chris@10: let scale_factor = inverse_int n in Chris@10: array n (fun i -> a i @* (scale_factor @* b i)) in Chris@10: Chris@10: (* generates a convolution using ffts. (all arguments are the Chris@10: same as to gen_convolution, below) *) Chris@10: let gen_convolution_by_fft n a b addtoall = Chris@10: let fft_a = dft 1 n a Chris@10: and fft_b = dft 1 n b in Chris@10: Chris@10: let fft_ab = make_product n fft_a fft_b Chris@10: and dc_term i = if (i == 0) then addtoall else zero in Chris@10: Chris@10: let fft_ab1 = array n (fun i -> fft_ab i @+ dc_term i) Chris@10: and sum = fft_a 0 in Chris@10: let conv = dft (-1) n fft_ab1 in Chris@10: (sum, conv) Chris@10: Chris@10: (* alternate routine for convolution. Seems to work better for Chris@10: small sizes. I have no idea why. *) Chris@10: and gen_convolution_by_fft_alt n a b addtoall = Chris@10: let ap = array n (fun i -> half (a i @+ a ((n - i) mod n))) Chris@10: and am = array n (fun i -> half (a i @- a ((n - i) mod n))) Chris@10: and bp = array n (fun i -> half (b i @+ b ((n - i) mod n))) Chris@10: and bm = array n (fun i -> half (b i @- b ((n - i) mod n))) Chris@10: in Chris@10: Chris@10: let fft_ap = dft 1 n ap Chris@10: and fft_am = dft 1 n am Chris@10: and fft_bp = dft 1 n bp Chris@10: and fft_bm = dft 1 n bm in Chris@10: Chris@10: let fft_abpp = make_product n fft_ap fft_bp Chris@10: and fft_abpm = make_product n fft_ap fft_bm Chris@10: and fft_abmp = make_product n fft_am fft_bp Chris@10: and fft_abmm = make_product n fft_am fft_bm Chris@10: and sum = fft_ap 0 @+ fft_am 0 Chris@10: and dc_term i = if (i == 0) then addtoall else zero in Chris@10: Chris@10: let fft_ab1 = array n (fun i -> (fft_abpp i @+ fft_abmm i) @+ dc_term i) Chris@10: and fft_ab2 = array n (fun i -> fft_abpm i @+ fft_abmp i) in Chris@10: let conv1 = dft (-1) n fft_ab1 Chris@10: and conv2 = dft (-1) n fft_ab2 in Chris@10: let conv = array n (fun i -> Chris@10: conv1 i @+ conv2 i) in Chris@10: (sum, conv) Chris@10: Chris@10: (* generator of assignment list assigning conv to the convolution of Chris@10: a and b, all of which are of length n. addtoall is added to Chris@10: all of the elements of the result. Returns (sum, convolution) pair Chris@10: where sum is the sum of the elements of a. *) Chris@10: Chris@10: in let gen_convolution = Chris@10: if (p <= !Magic.alternate_convolution) then Chris@10: gen_convolution_by_fft_alt Chris@10: else Chris@10: gen_convolution_by_fft Chris@10: Chris@10: (* fft generator for prime n = p using Rader's algorithm for Chris@10: turning the fft into a convolution, which then can be Chris@10: performed in a variety of ways *) Chris@10: in Chris@10: let g = find_generator p in Chris@10: let ginv = pow_mod g (p - 2) p in Chris@10: let input_perm = array p (fun i -> input (pow_mod g i p)) Chris@10: and omega_perm = array p (fun i -> exp p (sign * (pow_mod ginv i p))) Chris@10: and output_perm = array p (fun i -> pow_mod ginv i p) Chris@10: in let (sum, conv) = Chris@10: (gen_convolution (p - 1) input_perm omega_perm (input 0)) Chris@10: in array p (fun i -> Chris@10: if (i = 0) then Chris@10: input 0 @+ sum Chris@10: else Chris@10: let i' = suchthat 0 (fun i' -> i = output_perm i') Chris@10: in conv i') Chris@10: Chris@10: (* our modified version of the conjugate-pair split-radix algorithm, Chris@10: which reduces the number of multiplications by rescaling the Chris@10: sub-transforms (power-of-two n's only) *) Chris@10: and newsplit sign n input = Chris@10: let rec s n k = (* recursive scale factor *) Chris@10: if n <= 4 then Chris@10: one Chris@10: else Chris@10: let k4 = (abs k) mod (n / 4) in Chris@10: let k4' = if k4 <= (n / 8) then k4 else (n/4 - k4) in Chris@10: (s (n / 4) k4') @* (real (exp n k4')) Chris@10: Chris@10: and sinv n k = (* 1 / s(n,k) *) Chris@10: if n <= 4 then Chris@10: one Chris@10: else Chris@10: let k4 = (abs k) mod (n / 4) in Chris@10: let k4' = if k4 <= (n / 8) then k4 else (n/4 - k4) in Chris@10: (sinv (n / 4) k4') @* (sec n k4') Chris@10: Chris@10: in let sdiv2 n k = (s n k) @* (sinv (2*n) k) (* s(n,k) / s(2*n,k) *) Chris@10: and sdiv4 n k = (* s(n,k) / s(4*n,k) *) Chris@10: let k4 = (abs k) mod n in Chris@10: sec (4*n) (if k4 <= (n / 2) then k4 else (n - k4)) Chris@10: Chris@10: in let t n k = (exp n k) @* (sdiv4 (n/4) k) Chris@10: Chris@10: and dft1 input = input Chris@10: and dft2 input = array 2 (fun k -> (input 0) @+ ((input 1) @* exp 2 k)) Chris@10: Chris@10: in let rec newsplit0 sign n input = Chris@10: if (n == 1) then dft1 input Chris@10: else if (n == 2) then dft2 input Chris@10: else let u = newsplit0 sign (n / 2) (fun i -> input (i*2)) Chris@10: and z = newsplitS sign (n / 4) (fun i -> input (i*4 + 1)) Chris@10: and z' = newsplitS sign (n / 4) (fun i -> input ((n + i*4 - 1) mod n)) Chris@10: and twid = array n (fun k -> s (n/4) k @* exp n (sign * k)) in Chris@10: let w = array n (fun k -> twid k @* z (k mod (n / 4))) Chris@10: and w' = array n (fun k -> conj (twid k) @* z' (k mod (n / 4))) in Chris@10: let ww = array n (fun k -> w k @+ w' k) in Chris@10: array n (fun k -> u (k mod (n / 2)) @+ ww k) Chris@10: Chris@10: and newsplitS sign n input = Chris@10: if (n == 1) then dft1 input Chris@10: else if (n == 2) then dft2 input Chris@10: else let u = newsplitS2 sign (n / 2) (fun i -> input (i*2)) Chris@10: and z = newsplitS sign (n / 4) (fun i -> input (i*4 + 1)) Chris@10: and z' = newsplitS sign (n / 4) (fun i -> input ((n + i*4 - 1) mod n)) in Chris@10: let w = array n (fun k -> t n (sign * k) @* z (k mod (n / 4))) Chris@10: and w' = array n (fun k -> conj (t n (sign * k)) @* z' (k mod (n / 4))) in Chris@10: let ww = array n (fun k -> w k @+ w' k) in Chris@10: array n (fun k -> u (k mod (n / 2)) @+ ww k) Chris@10: Chris@10: and newsplitS2 sign n input = Chris@10: if (n == 1) then dft1 input Chris@10: else if (n == 2) then dft2 input Chris@10: else let u = newsplitS4 sign (n / 2) (fun i -> input (i*2)) Chris@10: and z = newsplitS sign (n / 4) (fun i -> input (i*4 + 1)) Chris@10: and z' = newsplitS sign (n / 4) (fun i -> input ((n + i*4 - 1) mod n)) in Chris@10: let w = array n (fun k -> t n (sign * k) @* z (k mod (n / 4))) Chris@10: and w' = array n (fun k -> conj (t n (sign * k)) @* z' (k mod (n / 4))) in Chris@10: let ww = array n (fun k -> (w k @+ w' k) @* (sdiv2 n k)) in Chris@10: array n (fun k -> u (k mod (n / 2)) @+ ww k) Chris@10: Chris@10: and newsplitS4 sign n input = Chris@10: if (n == 1) then dft1 input Chris@10: else if (n == 2) then Chris@10: let f = dft2 input Chris@10: in array 2 (fun k -> (f k) @* (sinv 8 k)) Chris@10: else let u = newsplitS2 sign (n / 2) (fun i -> input (i*2)) Chris@10: and z = newsplitS sign (n / 4) (fun i -> input (i*4 + 1)) Chris@10: and z' = newsplitS sign (n / 4) (fun i -> input ((n + i*4 - 1) mod n)) in Chris@10: let w = array n (fun k -> t n (sign * k) @* z (k mod (n / 4))) Chris@10: and w' = array n (fun k -> conj (t n (sign * k)) @* z' (k mod (n / 4))) in Chris@10: let ww = array n (fun k -> w k @+ w' k) in Chris@10: array n (fun k -> (u (k mod (n / 2)) @+ ww k) @* (sdiv4 n k)) Chris@10: Chris@10: in newsplit0 sign n input Chris@10: Chris@10: and dft sign n input = Chris@10: let rec cooley_tukey sign n1 n2 input = Chris@10: let tmp1 = Chris@10: array n2 (fun i2 -> Chris@10: dft sign n1 (fun i1 -> input (i1 * n2 + i2))) in Chris@10: let tmp2 = Chris@10: array n1 (fun i1 -> Chris@10: array n2 (fun i2 -> Chris@10: exp n (sign * i1 * i2) @* tmp1 i2 i1)) in Chris@10: let tmp3 = array n1 (fun i1 -> dft sign n2 (tmp2 i1)) in Chris@10: (fun i -> tmp3 (i mod n1) (i / n1)) Chris@10: Chris@10: (* Chris@10: * This is "exponent -1" split-radix by Dan Bernstein. Chris@10: *) Chris@10: and split_radix_dit sign n input = Chris@10: let f0 = dft sign (n / 2) (fun i -> input (i * 2)) Chris@10: and f10 = dft sign (n / 4) (fun i -> input (i * 4 + 1)) Chris@10: and f11 = dft sign (n / 4) (fun i -> input ((n + i * 4 - 1) mod n)) in Chris@10: let g10 = array n (fun k -> Chris@10: exp n (sign * k) @* f10 (k mod (n / 4))) Chris@10: and g11 = array n (fun k -> Chris@10: exp n (- sign * k) @* f11 (k mod (n / 4))) in Chris@10: let g1 = array n (fun k -> g10 k @+ g11 k) in Chris@10: array n (fun k -> f0 (k mod (n / 2)) @+ g1 k) Chris@10: Chris@10: and split_radix_dif sign n input = Chris@10: let n2 = n / 2 and n4 = n / 4 in Chris@10: let x0 = array n2 (fun i -> input i @+ input (i + n2)) Chris@10: and x10 = array n4 (fun i -> input i @- input (i + n2)) Chris@10: and x11 = array n4 (fun i -> Chris@10: input (i + n4) @- input (i + n2 + n4)) in Chris@10: let x1 k i = Chris@10: exp n (k * i * sign) @* (x10 i @+ exp 4 (k * sign) @* x11 i) in Chris@10: let f0 = dft sign n2 x0 Chris@10: and f1 = array 4 (fun k -> dft sign n4 (x1 k)) in Chris@10: array n (fun k -> Chris@10: if k mod 2 = 0 then f0 (k / 2) Chris@10: else let k' = k mod 4 in f1 k' ((k - k') / 4)) Chris@10: Chris@10: and prime_factor sign n1 n2 input = Chris@10: let tmp1 = array n2 (fun i2 -> Chris@10: dft sign n1 (fun i1 -> input ((i1 * n2 + i2 * n1) mod n))) Chris@10: in let tmp2 = array n1 (fun i1 -> Chris@10: dft sign n2 (fun k2 -> tmp1 k2 i1)) Chris@10: in fun i -> tmp2 (i mod n1) (i mod n2) Chris@10: Chris@10: in let algorithm sign n = Chris@10: let r = choose_factor n in Chris@10: if List.mem n !Magic.rader_list then Chris@10: (* special cases *) Chris@10: dft_rader sign n Chris@10: else if (r == 1) then (* n is prime *) Chris@10: dft_prime sign n Chris@10: else if (gcd r (n / r)) == 1 then Chris@10: prime_factor sign r (n / r) Chris@10: else if (n mod 4 = 0 && n > 4) then Chris@10: if !Magic.newsplit && is_power_of_two n then Chris@10: newsplit sign n Chris@10: else if !Magic.dif_split_radix then Chris@10: split_radix_dif sign n Chris@10: else Chris@10: split_radix_dit sign n Chris@10: else Chris@10: cooley_tukey sign r (n / r) Chris@10: in Chris@10: array n (algorithm sign n input)