This script contains implementations of the discrete Fourier transform, dft -- an quadratic implementation close to the specification fft -- an O(n.logn) recursive implementation vft -- a sketch of an optimisation of the same as described in the accompanying paper. There are also code to exercise them and examples of its use. gj 3.ix.1989 Operator declarations ~~~~~~~~~~~~~~~~~~~~~ The odd graphics which are used as operators have to be introduced at the top of the script. Although * is used thoughout this script for `map', it has to be bound in local definitions because its use has been pre-empted by predefined numeric multiplication. Otherwise it would have been defined here by %right 90 * (*) :: (a -> b) -> [a] -> [b] (*) = map The triangle operator satisfies (f /\ x)!i = (f^i) (x!i), but is defined recursively > %right 90 /\ > (/\) :: (a -> a) -> [a] -> [a] > f /\ [] = [] > f /\ (x:xs) = x : (f /\ (f * xs)) > where (*) = map Transposition satisfies (|/| xss)!i!j = xss!j!i > %prefix |/| > (|/|) :: [[a]] -> [[a]] > |/| xss = [ [x] | x <- hd xss ], if #xss = 1 > = [ x:xs | (x,xs) <- zip (hd xss, |/| (tl xss)) ], if #xss > 1 The operator +/ will be used to stand for a various addition-reductions, and although no value is given here, it needs to be declared to be a prefix operator > %prefix +/ List catenation and its inverse ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Catenation of lists, written ++/ in the paper, is implemented by the primitive concat = foldr (++) [] for speed > %prefix ++/ > (++/) :: [[a]] -> [a] > (++/) = concat It is inverted by two functions, `by' and `into': (by n xs) is used only if #xs is divisible by n, and returns a list of lists of length n which when concatenated return xs > by :: num -> [a] -> [[a]] > by n xs = [], if xs = [] > = take n xs : by n (drop n xs), if #xs >= n and (into n xs) is used only if #xs is divisible by n, and returns a list of length n of lists which when concatenated return xs > into :: num -> [a] -> [[a]] > into n xs = by (#xs $div n) xs Instances of by and into will be represented by the \++ operator (which is written as the inverse of ++/ in the paper) here introduced by > %prefix \++ Fourier transform ~~~~~~~~~~~~~~~~~ The Fourier transform of xs is obtained by summing xs_i.w^(i.j); the transform is implemented by (dft (+/) ($igma)) where (n $igma) is a function that represents multiplication by the n-th root of unity, and ((+/) xs) represents the sum of the elements of xs. > dft :: ([a] -> a) -> (num -> a -> a) -> ([a] -> [a]) > dft (+/) ($igma) xs = (((+/) *) . (((n $igma) /\) /\) . (n $copy)) xs > where n = #xs > (*) = map and the fast transform is obtained by a divide-and-conquer strategy applied to composite width vectors, falling back on the use of dft in the other cases > fft :: ([a] -> a) -> (num -> a -> a) -> ([a] -> [a]) > fft (+/) ($igma) = f > where f xs = f' xs, if prime (#xs) > = f'' xs, otherwise > f' = dft (+/) ($igma) > f'' xs = ((++/) . (|/|) . (f *) . > (((n $igma) /\) /\) . > (|/|) . (f *) . (|/|) . (\++)) xs > where n = #xs > (\++) = into (factor n) > (*) = map Here, the factor of n chosen for division of xs is (factor n) which is the largest factor of n less than its square root, that is > factor n = (last . divides n . belowroot n) [1..n] > where divides n = filter ((= 0) . (n $mod)) > belowroot n = takewhile ((<= n) . sq) > where sq x = x * x and n is prime iff this factor is one. > prime n = factor n = 1 (I know it could be made more efficient, but I was writing for clarity!) You might try to replace the term ((|/|) . (f *) . (|/|)), according to the paper, by (fft (((+/) *) . (|/|)) ((*) . ($igma))), and indeed vft :: ([a] -> a) -> (num -> a -> a) -> ([a] -> [a]) vft (+/) ($igma) = f where f xs = f' xs, if prime (#xs) = f'' xs, otherwise f' = dft (+/) ($igma) f'' xs = ((++/) . (|/|) . (f *) . (((n $igma) /\) /\) . f''' . (\++)) xs where n = #xs (\++) = into (factor n) (*) = map f''' = vft vsum vsigma where vsum = ((+/) *) . (|/|) where (*) = map vsigma = (*) . ($igma) where (*) = map would work were it not for the fact that the recursive call leads to an ill-founded type recursion: the arguments of the recursive call of vft oblige it to have a type ([[a]] -> [a]) -> (num -> [a] -> [a]) -> ([[a]] -> [[a]]) You can make a definition which uses this algorithm for one factorisation and then falls back on the earlier algorithm: > vft :: ([a] -> a) -> (num -> a -> a) -> ([a] -> [a]) > vft (+/) ($igma) = f > where f xs = f' xs, if prime (#xs) > = f'' xs, otherwise > f' = dft (+/) ($igma) > f'' xs = ((++/) . (|/|) . (f *) . > (((n $igma) /\) /\) . > f''' . (\++)) xs > where n = #xs > (\++) = into (factor n) > (*) = map > f''' = fft vsum vsigma > where vsum = ((+/) *) . (|/|) > where (*) = map > vsigma = (*) . ($igma) > where (*) = map The large number of bindings of (*) are again an artefact of the type checking, since the uses have different types. A symbolic example ~~~~~~~~~~~~~~~~~~ The arithmetic in the example is going to be over formal polynomials in a given root of unity, that root being the smallest one to be used in a given program. The polynomial -- necessarily of finite degree -- is represented by a list of its coefficients, in decreasing order of power of the root. > poly a == [a] Multiplication of such a polynomial by the given root of unity moves the coefficient of the highest power of the root to the end of the list -- as the coefficient of unity in the result -- and shifts the other coefficients left by one place. > shift :: [a] -> [a] > shift (x:xs) = xs ++ [x] If p divides #xs, then (p $sigma xs) represents xs times the p-th principal root of unity: this is achieved by shifting the coefficients by one p-th of their length. > ($igma) :: num -> poly a -> poly a > (p $igma) = (++/) . shift . into p The sum of two formal polynomial values is the pointwise sum of the coefficients; since each of these coefficients will be represented by a list of terms, they can just be concatenated to represent addition. The sum of a list of polynomials can either be implemented by a fold of a zipwith, or more conveniently for the present purposes by the equivalent > (+/) :: [poly a] -> poly a > (+/) = ((++/) *) . (|/|) where (*) = map Demonstration code ~~~~~~~~~~~~~~~~~~ To demonstrate the symbolic Fourier transform code, it will be exercised on polynomials, of degree less than n, over lists of integers; the presence of the number i in the j-th coefficient of a polynomial (indexed of course from the right-hand end) will be taken to represent a formal variable x_i times the j-th power of w -- a formal n-th root of unity. > index == num > coeff == [index] A vector of length n, of formal polynomials of degree less than n, will be laid out by translating each of the polynomials into strings and laying them out one to a line. > showvec :: [poly coeff] -> string > showvec = layn . map showpoly Each polynomial is represented by displaying the terms with non-zero (non-empty) coefficients, and punctuating with `+' signs. > showpoly :: poly coeff -> string > showpoly = stringfold " + " . map showterm . index > where index xs > = [ (x,i) | (x,i) <- zip(reverse xs,[0..]); x ~= []] Each term is represented by the coefficient times `w' to a power. > showterm :: (coeff, num) -> string > showterm (x,i) = showcoeff x, if i = 0 > = bracket showcoeff x ++ ".w^" ++ show i, if i > 0 > where bracket f x = f x, if #x = 1 > = "(" ++ f x ++ ")", if #x > 1 Each coefficient is represented as a list of variables punctuated with `+' signs. The call of `sort' arranges that the variables appear in order of index -- it happens that no variable will appear twice in any coefficient, but in principle this function could collect the multiplicity of a variable and pass it to `showvar'. > showcoeff :: coeff -> string > showcoeff = stringfold "+" . map showvar . sort Each variable is represented by `x' followed by its index. > showvar :: index -> string > showvar i = "x" ++ show i The function `stringfold' inserts the representation of an operator in the places in the representation of an expression corresponding to a reduction of that operator. > stringfold :: string -> [string] -> string > stringfold op = foldr1 (infix op) > where infix op l r = l ++ op ++ r The input to which an n-point transform will be applied is to represent a list of n polynomials, the i-th of which represents x_i (as the coefficient of unity in a formal polynomial of degree n-1). > input :: num -> [poly num] > input n = [ nulls ++ [[i]] | i <- [1..n] ] where nulls = copy (n-1) [] Execution of the example code ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ For example, the input for an eight-point transform is demonstrated by ? (showvec . input) 8 1) x1 2) x2 3) x3 4) x4 5) x5 6) x6 7) x7 8) x8 The specification applied to this input returns ? (showvec . dft (+/) ($igma) . input) 8 1) x1+x2+x3+x4+x5+x6+x7+x8 2) x1 + x2.w^1 + x3.w^2 + x4.w^3 + x5.w^4 + x6.w^5 + x7.w^6 + x8.w^7 3) x1+x5 + (x2+x6).w^2 + (x3+x7).w^4 + (x4+x8).w^6 4) x1 + x4.w^1 + x7.w^2 + x2.w^3 + x5.w^4 + x8.w^5 + x3.w^6 + x6.w^7 5) x1+x3+x5+x7 + (x2+x4+x6+x8).w^4 6) x1 + x6.w^1 + x3.w^2 + x8.w^3 + x5.w^4 + x2.w^5 + x7.w^6 + x4.w^7 7) x1+x5 + (x4+x8).w^2 + (x3+x7).w^4 + (x2+x6).w^6 8) x1 + x8.w^1 + x7.w^2 + x6.w^3 + x5.w^4 + x4.w^5 + x3.w^6 + x2.w^7 as do the similar expressions in fft and vft. For what it is worth, on a Sun3 the costs are approximately: (showvec.dft(+/)($igma).input) 8 66 CPU seconds, 166 000 reductions (showvec.fft(+/)($igma).input) 8 4.8 seconds, 9 600 reductions (showvec.vft(+/)($igma).input) 8 5.0 seconds, 9 400 reductions NB: these costs will not scale as O(n^2) and O(n.log n) because direct implementation of ((w /\) /\) requires O(n^4) applications of `w' each of which might cost O(n). In a practical implementation it would be necessary to do something more sensible.