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 /tests/trunk/bugs/tests.obsolete/bug1227.sml
ViewVC logotype

Annotation of /tests/trunk/bugs/tests.obsolete/bug1227.sml

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2071 - (view) (download)
Original Path: tests/trunk/bugs/tests/bug1227.sml

1 : dbm 548 (* bug1227.sml *)
2 :    
3 :     local
4 :     val sin = Real.Math.sin
5 :     val cos = Real.Math.cos
6 :     open Array
7 :     fun abs_real(x:real) = if x < 0.0 then ~x else x
8 :     fun print_real (x:real) = print(Real.toString x)
9 :    
10 :     val pi = 3.14159265358979323846
11 :    
12 :     val tpi = 2.0 * pi
13 :    
14 :     fun for(start,stop,f) =
15 :     let fun loop i = if i > stop then () else (f i; loop (i+1))
16 :     in
17 :     loop start
18 :     end
19 :    
20 :     val print_string : string -> unit = print
21 :     val print_int : int -> unit = print_string o Int.toString
22 :     val print_newline : unit -> unit = fn _ => print_string "\n"
23 :    
24 :     fun dump pxr pxi =
25 :     for(0,15,fn i =>
26 :     (print_int i;
27 :     print " ";
28 :     print_real (sub(pxr,i+1));
29 :     print " ";
30 :     print_real (sub(pxi,i+1));
31 :     print_newline()))
32 :    
33 :    
34 :     fun fft px py np =
35 :     let val i = ref 2
36 :     val m = ref 1
37 :    
38 :     val _ =
39 :     while (!i < np) do
40 :     (i := !i + !i;
41 :     m := !m + 1)
42 :     val n = !i
43 :     in
44 :     if n <> np then (
45 :     for (np+1,n,fn i=>
46 :     (update(px,i,0.0);
47 :     update(py,i,0.0)));
48 :     print_string "Use "; print_int n;
49 :     print_string " point fft"; print_newline()
50 :     ) else ();
51 :    
52 :     let val n2 = ref(n+n)
53 :     in
54 :     for(1,!m-1,fn k =>
55 :     let val _ = n2 := !n2 div 2
56 :     val n4 = !n2 div 4
57 :     val e = tpi / (real (!n2 ))
58 :     in
59 :     for(1,n4,fn j =>
60 :     let val a = e * real(j - 1)
61 :     val a3 = 3.0 * a
62 :     val cc1 = cos(a)
63 :     val ss1 = sin(a)
64 :     val cc3 = cos(a3)
65 :     val ss3 = sin(a3)
66 :     val is = ref j
67 :     val id = ref(2 * !n2)
68 :     in
69 :     while !is < n do
70 :     let val i0r = ref (!is)
71 :     in
72 :     while !i0r < n do
73 :     let val i0 = !i0r
74 :     val i1 = i0 + n4
75 :     val i2 = i1 + n4
76 :     val i3 = i2 + n4
77 :     val r1 = sub(px,i0) - sub(px,i2)
78 :     val _ = update(px,i0,sub(px,i0) + sub(px,i2))
79 :     val r2 = sub(px,i1) - sub(px,i3)
80 :     val _ = update(px,i1,sub(px,i1) + sub(px,i3))
81 :     val s1 = sub(py,i0) - sub(py,i2)
82 :     val _ = update(py,i0,sub(py,i0) + sub(py,i2))
83 :     val s2 = sub(py,i1) - sub(py,i3)
84 :     val _ = update(py,i1,sub(py,i1) + sub(py,i3))
85 :     val s3 = r1 - s2
86 :     val r1 = r1 + s2
87 :     val s2 = r2 - s1
88 :     val r2 = r2 + s1
89 :     in
90 :     update(px,i2,r1*cc1 - s2*ss1);
91 :     update(py,i2,~s2*cc1 - r1*ss1);
92 :     update(px,i3,s3*cc3 + r2*ss3);
93 :     update(py,i3,r2*cc3 - s3*ss3);
94 :     i0r := i0 + !id
95 :     end;
96 :     is := 2 * !id - !n2 + j;
97 :     id := 4 * !id
98 :     (* ; dump px py *)
99 :     end
100 :     end)
101 :     end)
102 :     end;
103 :    
104 :     (************************************)
105 :     (* Last stage, length=2 butterfly *)
106 :     (************************************)
107 :    
108 :     let val is = ref 1
109 :     val id = ref 4
110 :     in
111 :     while !is < n do
112 :     let val i0r = ref (!is)
113 :     in
114 :     while !i0r <= n do
115 :     let val i0 = !i0r
116 :     val i1 = i0 + 1
117 :     val r1 = sub(px,i0)
118 :     val _ = update(px,i0,r1 + sub(px,i1))
119 :     val _ = update(px,i1,r1 - sub(px,i1))
120 :     val r1 = sub(py,i0)
121 :     in
122 :     update(py,i0,r1 + sub(py,i1));
123 :     update(py,i1,r1 - sub(py,i1));
124 :     i0r := i0 + !id
125 :     end;
126 :     is := 2 * !id - 1;
127 :     id := 4 * !id
128 :     end
129 :     end;
130 :     (*
131 :     print "\nbutterfly\n";
132 :     dump px py;
133 :     *)
134 :     (*************************)
135 :     (* Bit reverse counter *)
136 :     (*************************)
137 :    
138 :     let val j = ref 1
139 :     in
140 :     for (1,n-1,fn i =>
141 :     (if i < !j then (
142 :     let val xt1 = sub(px,!j)
143 :     val xt2 = sub(px,i)
144 :     val _ = update(px,!j,xt2)
145 :     val _ = update(px,i,xt1)
146 :     val yt1 = sub(py,!j)
147 :     val yt2 = sub(py,i)
148 :     in
149 :     update(py,!j,yt2);
150 :     update(py,i,yt1)
151 :     end)
152 :     else ();
153 :     let val k = ref(n div 2)
154 :     in
155 :     while !k < !j do (j := !j - !k; k := !k div 2);
156 :     j := !j + !k
157 :     end));
158 :     (*
159 :     print "\nbit reverse\n";
160 :     dump px py;
161 :     *)
162 :     n
163 :     end
164 :     end;
165 :    
166 :     fun test np =
167 :     (print_int np; print_string "... ";
168 :     let val enp = real np
169 :     val npm = np div 2 - 1
170 :     val pxr = array ((np+2), 0.0)
171 :     val pxi = array ((np+2), 0.0)
172 :     val t = pi / enp
173 :     val _ = update(pxr,1,(enp - 1.0) * 0.5)
174 :     val _ = update(pxi,1,0.0)
175 :     val n2 = np div 2
176 :     val _ = update(pxr,n2+1,~0.5)
177 :     val _ = update(pxi,n2+1,0.0)
178 :     in
179 :     for (1,npm,fn i =>
180 :     let val j = np - i
181 :     val _ = update(pxr,i+1,~0.5)
182 :     val _ = update(pxr,j+1,~0.5)
183 :     val z = t * (real i)
184 :     val y = ~0.5 * (cos(z)/sin(z))
185 :     in
186 :     update(pxi,i+1,y);
187 :     update(pxi,j+1,~y)
188 :     end);
189 :     (*
190 :     print "\n"; print "before fft: \n";
191 :     dump pxr pxi;
192 :     *)
193 :     fft pxr pxi np;
194 :     (*
195 :     print "\n"; print "after fft: \n";
196 :     dump pxr pxi;
197 :     *)
198 :     let val zr = ref 0.0
199 :     val zi = ref 0.0
200 :     val kr = ref 0
201 :     val ki = ref 0
202 :     in
203 :     for (0,np-1,fn i =>
204 :     let val a = abs_real(sub(pxr,i+1) - (real i))
205 :     in
206 :     if !zr < a then
207 :     (zr := a;
208 :     kr := i)
209 :     else ();
210 :     let val a = abs_real(sub(pxi,i+1))
211 :     in
212 :     if !zi < a then
213 :     (zi := a;
214 :     ki := i)
215 :     else ()
216 :     end
217 :     end);
218 :     let val zm = if abs_real (!zr) < abs_real (!zi) then !zi else !zr
219 :     in
220 :     print_real zm; print_newline()
221 :     end
222 :     end
223 :     end)
224 :    
225 :     fun doit() =
226 :     let val np = ref 16
227 :     in for(1,13,fn i => (test (!np); np := (!np)*2))
228 :     end
229 :    
230 :     in
231 :     val _ = doit()
232 :     end

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