Home My Page Projects Code Snippets Project Openings SML/NJ
Summary Activity Forums Tracker Lists Tasks Docs Surveys News SCM Files

SCM Repository

[smlnj] Annotation of /sml/trunk/benchmarks/programs/fft/fft.sml
ViewVC logotype

Annotation of /sml/trunk/benchmarks/programs/fft/fft.sml

Parent Directory Parent Directory | Revision Log Revision Log


Revision 193 - (view) (download)

1 : monnier 193 structure Main: BMARK = struct
2 :    
3 :     local
4 :     open Array Math
5 :    
6 :     val printr = print o Real.toString
7 :     val printi = print o Int.toString
8 :     in
9 :    
10 :     val PI = 3.14159265358979323846
11 :    
12 :     val tpi = 2.0 * PI
13 :    
14 :     fun fft px py np =
15 :     let fun find_num_points i m =
16 :     if i < np then find_num_points (i+i) (m+1) else (i,m)
17 :     val (n,m) = find_num_points 2 1 in
18 :     if n <> np then
19 :     let fun loop i = if i > n then () else (
20 :     update(px, i, 0.0);
21 :     update(py, i, 0.0);
22 :     loop (i+1))
23 :     in
24 :     loop (np+1);
25 :     print "Use "; printi n; print " point fft\n"
26 :     end
27 :     else ();
28 :    
29 :     let fun loop_k k n2 = if k >= m then () else
30 :     let val n4 = n2 div 4
31 :     val e = tpi / (real n2)
32 :     fun loop_j j a = if j > n4 then () else
33 :     let val a3 = 3.0 * a
34 :     val cc1 = cos(a)
35 :     val ss1 = sin(a)
36 :     val cc3 = cos(a3)
37 :     val ss3 = sin(a3)
38 :     fun loop_is is id = if is >= n then () else
39 :     let fun loop_i0 i0 = if i0 >= n then () else
40 :     let val i1 = i0 + n4
41 :     val i2 = i1 + n4
42 :     val i3 = i2 + n4
43 :     val r1 = sub(px, i0) - sub(px, i2)
44 :     val _ = update(px, i0, sub(px, i0) + sub(px, i2))
45 :     val r2 = sub(px, i1) - sub(px, i3)
46 :     val _ = update(px, i1, sub(px, i1) + sub(px, i3))
47 :     val s1 = sub(py, i0) - sub(py, i2)
48 :     val _ = update(py, i0, sub(py, i0) + sub(py, i2))
49 :     val s2 = sub(py, i1) - sub(py, i3)
50 :     val _ = update(py, i1, sub(py, i1) + sub(py, i3))
51 :     val s3 = r1 - s2
52 :     val r1 = r1 + s2
53 :     val s2 = r2 - s1
54 :     val r2 = r2 + s1
55 :     val _ = update(px, i2, r1*cc1 - s2*ss1)
56 :     val _ = update(py, i2, ~s2*cc1 - r1*ss1)
57 :     val _ = update(px, i3, s3*cc3 + r2*ss3)
58 :     val _ = update(py, i3, r2*cc3 - s3*ss3)
59 :     in
60 :     loop_i0 (i0 + id)
61 :     end
62 :     in
63 :     loop_i0 is;
64 :     loop_is (2 * id - n2 + j) (4 * id)
65 :     end
66 :     in
67 :     loop_is j (2 * n2);
68 :     loop_j (j+1) (e * real j)
69 :     end
70 :     in
71 :     loop_j 1 0.0;
72 :     loop_k (k+1) (n2 div 2)
73 :     end
74 :     in
75 :     loop_k 1 n
76 :     end;
77 :    
78 :     (************************************)
79 :     (* Last stage, length=2 butterfly *)
80 :     (************************************)
81 :    
82 :     let fun loop_is is id = if is >= n then () else
83 :     let fun loop_i0 i0 = if i0 > n then () else
84 :     let val i1 = i0 + 1
85 :     val r1 = sub(px, i0)
86 :     val _ = update(px, i0, r1 + sub(px, i1))
87 :     val _ = update(px, i1, r1 - sub(px, i1))
88 :     val r1 = sub(py, i0)
89 :     val _ = update(py, i0, r1 + sub(py, i1))
90 :     val _ = update(py, i1, r1 - sub(py, i1))
91 :     in
92 :     loop_i0 (i0 + id)
93 :     end
94 :     in
95 :     loop_i0 is;
96 :     loop_is (2*id - 1) (4 * id)
97 :     end
98 :     in
99 :     loop_is 1 4
100 :     end;
101 :    
102 :     (*************************)
103 :     (* Bit reverse counter *)
104 :     (*************************)
105 :    
106 :     let fun loop_i i j = if i >= n then () else
107 :     (if i < j then
108 :     (let val xt = sub(px, j)
109 :     in update(px, j, sub(px, i)); update(px, i, xt)
110 :     end;
111 :     let val xt = sub(py, j)
112 :     in update(py, j, sub(py, i)); update(py, i, xt)
113 :     end)
114 :     else ();
115 :     let fun loop_k k j =
116 :     if k < j then loop_k (k div 2) (j-k) else j+k
117 :     val j' = loop_k (n div 2) j
118 :     in
119 :     loop_i (i+1) j'
120 :     end)
121 :     in
122 :     loop_i 1 1
123 :     end;
124 :    
125 :     n
126 :     end
127 :    
128 :     fun abs x = if x >= 0.0 then x else ~x
129 :    
130 :     fun test np =
131 :     let val _ = (printi np; print "... ")
132 :     val enp = real np
133 :     val npm = (np div 2) - 1
134 :     val pxr = array (np+2, 0.0)
135 :     val pxi = array (np+2, 0.0)
136 :     val t = PI / enp
137 :     val _ = update(pxr, 1, (enp - 1.0) * 0.5)
138 :     val _ = update(pxi, 1, 0.0)
139 :     val n2 = np div 2
140 :     val _ = update(pxr, n2+1, ~0.5)
141 :     val _ = update(pxi, n2+1, 0.0)
142 :     fun loop_i i = if i > npm then () else
143 :     let val j = np - i
144 :     val _ = update(pxr, i+1, ~0.5)
145 :     val _ = update(pxr, j+1, ~0.5)
146 :     val z = t * real i
147 :     val y = ~0.5*(cos(z)/sin(z))
148 :     val _ = update(pxi, i+1, y)
149 :     val _ = update(pxi, j+1, ~y)
150 :     in
151 :     loop_i (i+1)
152 :     end
153 :     val _ = loop_i 1
154 :     (***
155 :     val _ = print "\n"
156 :     fun loop_i i = if i > 15 then () else
157 :     (print i; print "\t";
158 :     print (sub(pxr, i+1)); print "\t";
159 :     print (sub(pxi, i+1)); print "\n"; loop_i (i+1))
160 :     val _ = loop_i 0
161 :     ***)
162 :     val _ = fft pxr pxi np
163 :     (***
164 :     fun loop_i i = if i > 15 then () else
165 :     (print i; print "\t";
166 :     print (sub(pxr, i+1)); print "\t";
167 :     print (sub(pxi, i+1)); print "\n"; loop_i (i+1))
168 :     val _ = loop_i 0
169 :     ***)
170 :     fun loop_i i zr zi kr ki = if i >= np then (zr,zi) else
171 :     let val a = abs(sub(pxr, i+1) - real i)
172 :     val (zr, kr) =
173 :     if zr < a then (a, i) else (zr, kr)
174 :     val a = abs(sub(pxi, i+1))
175 :     val (zi, ki) =
176 :     if zi < a then (a, i) else (zi, ki)
177 :     in
178 :     loop_i (i+1) zr zi kr ki
179 :     end
180 :     val (zr, zi) = loop_i 0 0.0 0.0 0 0
181 :     val zm = if abs zr < abs zi then zi else zr
182 :     in
183 :     printr zm; print "\n"
184 :     end
185 :    
186 :     fun loop_np i np = if i > 13 then () else
187 :     (test np; loop_np (i+1) (np*2))
188 :    
189 :     fun doit () = loop_np 1 16
190 :    
191 :     fun testit outstream = doit()
192 :    
193 :     end
194 :     end;

root@smlnj-gforge.cs.uchicago.edu
ViewVC Help
Powered by ViewVC 1.0.0