1 (* From the SML
/NJ benchmark suite
. *)
7 structure Main
: BMARK
= struct
12 val printr
= print
o (Real.fmt (StringCvt.SCI(SOME
14)))
13 val printi
= print
o Int.toString
16 val PI
= 3.14159265358979323846
31 fun trace2(name
,f
) (x
,y
) =
46 val op * = trace2("*", Real.* )
47 val op - = trace2("-", Real.-)
48 val op + = trace2("+", Real.+)
53 _overload
* : ('a
* 'a
-> 'a
)
57 _overload
+ : ('a
* 'a
-> 'a
)
61 _overload
- : ('a
* 'a
-> 'a
)
66 val sin
= trace("sin", sin
)
67 val cos
= trace("cos", cos
)
92 fun find_num_points i m
=
94 then find_num_points (i
+i
) (m
+1)
96 val (n
,m
) = find_num_points
2 1
97 (* val _
= (printi n
;
106 else (update(px
, i
, 0.0);
111 print
"Use "; printi n
; print
" point fft\n"
121 val e
= tpi
/ (real n2
)
142 val r1
= sub(px
, i0
) - sub(px
, i2
)
143 val _
= update(px
, i0
, sub(px
, i0
) + sub(px
, i2
))
144 val r2
= sub(px
, i1
) - sub(px
, i3
)
145 val _
= update(px
, i1
, sub(px
, i1
) + sub(px
, i3
))
146 val s1
= sub(py
, i0
) - sub(py
, i2
)
147 val _
= update(py
, i0
, sub(py
, i0
) + sub(py
, i2
))
148 val s2
= sub(py
, i1
) - sub(py
, i3
)
149 val _
= update(py
, i1
, sub(py
, i1
) + sub(py
, i3
))
154 val _
= update(px
, i2
, r1
*cc1
- s2
*ss1
)
155 val _
= update(py
, i2
, ~s2
*cc1
- r1
*ss1
)
156 val _
= update(px
, i3
, s3
*cc3
+ r2
*ss3
)
157 val _
= update(py
, i3
, r2
*cc3
- s3
*ss3
)
163 loop_is (2 * id
- n2
+ j
) (4 * id
)
167 loop_j (j
+1) (e
* real j
)
171 loop_k (k
+1) (n2
div 2)
177 (************************************)
178 (* Last stage
, length
=2 butterfly
*)
179 (************************************)
181 let fun loop_is is id
= if is
>= n
then () else
182 let fun loop_i0 i0
= if i0
> n
then () else
185 val _
= update(px
, i0
, r1
+ sub(px
, i1
))
186 val _
= update(px
, i1
, r1
- sub(px
, i1
))
188 val _
= update(py
, i0
, r1
+ sub(py
, i1
))
189 val _
= update(py
, i1
, r1
- sub(py
, i1
))
195 loop_is (2*id
- 1) (4 * id
)
201 (*************************)
202 (* Bit reverse counter
*)
203 (*************************)
211 then (let val xt
= sub(px
, j
)
212 in update(px
, j
, sub(px
, i
)); update(px
, i
, xt
)
214 let val xt
= sub(py
, j
)
215 in update(py
, j
, sub(py
, i
)); update(py
, i
, xt
)
220 if k
< j
then loop_k (k
div 2) (j
-k
) else j
+k
221 val j
' = loop_k (n
div 2) j
233 fun abs x
= if x
>= 0.0 then x
else ~x
236 let val _
= (printi np
; print
"... ")
238 val npm
= (np
div 2) - 1
239 val pxr
= array (np
+2, 0.0)
240 val pxi
= array (np
+2, 0.0)
242 val _
= update(pxr
, 1, (enp
- 1.0) * 0.5)
243 val _
= update(pxi
, 1, 0.0)
245 val _
= update(pxr
, n2
+1, ~
0.5)
246 val _
= update(pxi
, n2
+1, 0.0)
247 fun loop_i i
= if i
> npm
then () else
249 val _
= update(pxr
, i
+1, ~
0.5)
250 val _
= update(pxr
, j
+1, ~
0.5)
252 val y
= ~
0.5*(cos(z
)/sin(z
))
253 val _
= update(pxi
, i
+1, y
)
254 val _
= update(pxi
, j
+1, ~y
)
260 (* val _
= print
"\n"
261 fun loop_i i
= if i
> 15 then () else
262 (printi i
; print
"\t";
263 printr (sub(pxr
, i
+1)); print
"\t";
264 printr (sub(pxi
, i
+1)); print
"\n"; loop_i (i
+1))
267 val _
= fft pxr pxi np
269 fun loop_i i
= if i
> 15 then () else
270 (printi i
; print
"\t";
271 printr (sub(pxr
, i
+1)); print
"\t";
272 printr (sub(pxi
, i
+1)); print
"\n"; loop_i (i
+1))
275 fun loop_i i zr zi kr ki
= if i
>= np
then (zr
,zi
) else
276 let val a
= abs(sub(pxr
, i
+1) - real i
)
278 if zr
< a
then (a
, i
) else (zr
, kr
)
279 val a
= abs(sub(pxi
, i
+1))
281 if zi
< a
then (a
, i
) else (zi
, ki
)
283 loop_i (i
+1) zr zi kr ki
285 val (zr
, zi
) = loop_i
0 0.0 0.0 0 0
286 val zm
= if abs zr
< abs zi
then zi
else zr
288 printr zm
; print
"\n"
291 fun loop_np i np
= if i
> 15 then () else
292 (test np
; loop_np (i
+1) (np
*2))
297 else (loop_np
1 256; doit (n
- 1))