Remove deprecated _scc_ (#8170)
[nofib.git] / parallel / OLD / NESL / fft.lhs
1 Time-stamp: <Wed May 22 1996 19:30:00 Stardate: [-31]7543.85 hwloidl>
2
3 Computing a Fast Fourier Transformation.
4
5 Based on the NESL code presented in:
6    Programming Parallel Algorithms
7    by Guy E. Blelloch 
8    in CACM 39(3), March 1996
9    URL: http://www.cs.cmu.edu/afs/cs.cmu.edu/project/scandal/public/www/nesl/alg-numerical.html
10
11 The following description is taken from the original NESL code:
12
13   This is a simple parallel version of the standard sequential fast Fourier
14   transform algorithm. The algorithm does O(n log n) work and has O(log n)
15   depth.
16   
17   In the code we first give a general FFT that works with any commutative
18   ring. As well as the data, this FFT takes as arguments the n roots of unity
19   (an array w) and the +, * functions on the ring (add, mult). We then
20   specialize the algorithm to a particular commutative ring---the complex
21   numbers. Multiline text area Submit
22
23 NB: THIS IS UNTESTED CODE %--  HWL
24
25 \begin{code}
26 #if defined(GRAN)
27 import Strategies
28 #endif
29 #if defined(ARGS)
30 import GranRandom
31 import LibSystem       -- for getArgs
32 #endif
33
34 even_elts []     = []
35 even_elts [x]    = [x]
36 even_elts xs@(x:_) = x : (even_elts (drop 2 xs))
37
38 odd_elts [] = []
39 odd_elts (x:xs) = even_elts xs
40 \end{code}
41
42 The main function, parameterised with function for adding and multiplying 
43 two elements (with Haskell's class system that shouldn't be necessary).
44
45 \begin{code}
46 -- fft :: Integral a => [a] -> [a] -> [a]
47 {- SPECIALISE fft :: [Complex Double] -> [Complex Double] -> [Complex Double] -}
48 fft :: [Complex Double] -> [Complex Double] -> [Complex Double]
49 fft a w  
50   | length a <= 1  =  a
51   | otherwise      = let r0 = {-# SCC "head" #-} fft (even_elts a) (even_elts w)
52                          r1 = {-# SCC "head" #-} fft (odd_elts a) (even_elts w)
53                          z =  {-# SCC "zip3" #-}  zip3 (r0++r0) (r1++r1) w
54                      in 
55 #if defined(GRAN)
56                         parList rnf r0 `par`    
57                         parList rnf r1 `par`
58                         parList rwhnf z  `par`
59                         parList rnf     
60                         --_parGlobal_ 11# 11# 0# 0#  (rnf r0) $
61                         --_parGlobal_ 12# 12# 0# 0#  (rnf r1) $
62                         --_parGlobal_ 13# 13# 0# 0#  (rnf z) $
63 #endif
64                         [ r0'+r1'*w' 
65                         | (r0',r1',w') <- z ]
66
67 complex_fft :: [Complex Double] -> [Complex Double]
68 complex_fft a =
69   let 
70       c :: Double  
71       c = (2.0*pi)/(fromIntegral (length a))
72       w = {-# SCC "w" #-}  [ (cos (c*(fromIntegral i)) :+ sin (c*(fromIntegral i)) ) 
73                            | i <- [0..length a] ]
74       -- add = \ (ar,ai) (br,bi) -> (ar+br,ai+bi)
75       -- mult = \ (ar,ai) (br,bi) -> (ar*br-ai*bi,ar*bi+ai*br)
76   in (rnf w) `seq` fft a w 
77 \end{code}
78
79 Test data.
80
81 \begin{code}
82 x0 :: [Complex Double]
83 x0 = [(2.0 :+ 0.0),(-1.0 :+ 1.0),(0.0 :+ 0.0),(-1.0 :+ -1.0)]
84
85 -- 64 elems with values between 0 and 99
86 x1 :: [Complex Double]
87 x1 = [ (67.0 :+ -17.0), (80.0 :+ 53.0), (-58.0 :+ 45.0), (16.0 :+ 96.0), (-15.0 :+ -4.0), (83.0 :+ 91.0), (-31.0 :+ -44.0), (96.0 :+ -37.0), (-63.0 :+ -99.0), (70.0 :+ 1.0), (-12.0 :+ 38.0), (-57.0 :+ 5.0), (-84.0 :+ -94.0), (31.0 :+ -47.0), (70.0 :+ -11.0), (41.0 :+ 24.0), (-99.0 :+ -74.0), (-59.0 :+ -41.0), (58.0 :+ -83.0), (27.0 :+ 25.0), (-38.0 :+ -32.0), (78.0 :+ -18.0), (-95.0 :+ 37.0), (68.0 :+ -42.0), (-61.0 :+ 92.0), (-34.0 :+ 16.0), (-66.0 :+ -52.0), (39.0 :+ 20.0), (-1.0 :+ -45.0), (69.0 :+ 47.0), (39.0 :+ -78.0), (-78.0 :+ -12.0), (81.0 :+ -96.0), (-88.0 :+ 20.0), (31.0 :+ 20.0), (-22.0 :+ -46.0), (-54.0 :+ 59.0), (53.0 :+ 16.0), (-96.0 :+ 55.0), (44.0 :+ -77.0), (48.0 :+ -49.0), (83.0 :+ -6.0), (77.0 :+ 49.0), (20.0 :+ 46.0), (21.0 :+ 9.0), (-40.0 :+ 78.0), (-57.0 :+ -65.0), (-12.0 :+ 21.0), (7.0 :+ -74.0), (36.0 :+ 87.0), (-86.0 :+ 0.0), (-6.0 :+ 21.0), (-15.0 :+ -2.0), (45.0 :+ -26.0), (18.0 :+ -65.0), (12.0 :+ -39.0), (-36.0 :+ -68.0), (49.0 :+ 6.0), (95.0 :+ -92.0), (74.0 :+ -44.0), (-30.0 :+ -95.0), (-26.0 :+ 86.0), (84.0 :+ -61.0), (-64.0 :+ -28.0) ]
88
89 \end{code}
90
91 Main function.
92
93 \begin{code}
94 args_to_IntList a = if length a < 2
95                       then error "Usage: fft <lst-length> <max-elem>\n"
96                       else map (\ a1 -> fst ((readDec a1) !! 0)) a
97
98 #if defined(ARGS)
99 munch_args =    getArgs >>= \a ->
100                 return (args_to_IntList a) >>= \[n,m] ->
101                 getRandomDoubles (fromIntegral m) >>= \ random_list -> 
102                 let 
103                   (l1, random_list') = splitAt n random_list
104                   (l2, random_list'') = splitAt n random_list'
105                   x = zipWith (\ x y -> (x :+ y)) l1 l2
106                 in
107                 return (x)
108
109 #else
110 munch_args = return (x1)
111 #endif
112
113 #ifdef PRINT
114 main = munch_args >>= \ x -> print (complex_fft x)
115 #else
116 main = munch_args >>= \ x -> (rnf $ complex_fft x) `seq` putStr "Done\n"
117 #endif
118
119 \end{code}
120
121 -----------------------------------------------------------------------------
122
123 This is the original NESL code:
124
125 function fft(a,w,add,mult) =
126 if #a == 1 then a
127 else
128   let r = {fft(b, even_elts(w), add, mult):
129            b in [even_elts(a),odd_elts(a)]}
130   in {add(a, mult(b, w)):
131       a in r[0] ++ r[0]; 
132       b in r[1] ++ r[1];
133       w in w};
134
135 function complex_fft(a) =
136 let 
137     c = 2.*pi/float(#a);
138     w = {cos(c*float(i)),sin(c*float(i)) : i in [0:#a]};
139     add = ((ar,ai),(br,bi)) => (ar+br,ai+bi);
140     mult = ((ar,ai),(br,bi)) => (ar*br-ai*bi,ar*bi+ai*br);
141 in fft(a,w,add,mult);
142
143 complex_fft([(2.,0.),(-1.,1.),(0.,0.),(-1.,-1.)]);