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
 [smlnj] / sml / trunk / benchmarks / programs / fft / fft.sml

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

 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;