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

# SCM Repository

[smlnj] View of /tests/trunk/bugs/tests.obsolete/bug1227.sml
 [smlnj] / tests / trunk / bugs / tests.obsolete / bug1227.sml

# View of /tests/trunk/bugs/tests.obsolete/bug1227.sml

Tue Oct 31 02:58:12 2006 UTC (12 years, 8 months ago) by blume
Original Path: tests/trunk/bugs/tests/bug1227.sml
File size: 5071 byte(s)
`moved tests trunk`
```(* bug1227.sml *)

local
val sin = Real.Math.sin
val cos = Real.Math.cos
open Array
fun abs_real(x:real) = if x < 0.0 then ~x else x
fun print_real (x:real) = print(Real.toString x)

val pi = 3.14159265358979323846

val tpi = 2.0 * pi

fun for(start,stop,f) =
let fun loop i = if i > stop then () else (f i; loop (i+1))
in
loop start
end

val print_string : string -> unit = print
val print_int : int -> unit = print_string o Int.toString
val print_newline : unit -> unit = fn _ => print_string "\n"

fun dump pxr pxi =
for(0,15,fn i =>
(print_int i;
print " ";
print_real (sub(pxr,i+1));
print " ";
print_real (sub(pxi,i+1));
print_newline()))

fun fft px py np =
let val i = ref 2
val m = ref 1

val _ =
while (!i < np) do
(i := !i + !i;
m := !m + 1)
val n = !i
in
if n <> np then (
for (np+1,n,fn i=>
(update(px,i,0.0);
update(py,i,0.0)));
print_string "Use "; print_int n;
print_string " point fft"; print_newline()
) else ();

let val n2 = ref(n+n)
in
for(1,!m-1,fn k =>
let val _ = n2 := !n2 div 2
val n4 = !n2 div 4
val e  = tpi / (real (!n2 ))
in
for(1,n4,fn j =>
let val a = e * real(j - 1)
val a3 = 3.0 * a
val cc1 = cos(a)
val ss1 = sin(a)
val cc3 = cos(a3)
val ss3 = sin(a3)
val is = ref j
val id = ref(2 * !n2)
in
while !is < n do
let val i0r = ref (!is)
in
while !i0r < n do
let val i0 = !i0r
val i1 = i0 + n4
val i2 = i1 + n4
val i3 = i2 + n4
val r1 = sub(px,i0) - sub(px,i2)
val _ = update(px,i0,sub(px,i0) + sub(px,i2))
val r2 = sub(px,i1) - sub(px,i3)
val _ = update(px,i1,sub(px,i1) + sub(px,i3))
val s1 = sub(py,i0) - sub(py,i2)
val _ = update(py,i0,sub(py,i0) + sub(py,i2))
val s2 = sub(py,i1) - sub(py,i3)
val _ = update(py,i1,sub(py,i1) + sub(py,i3))
val s3 = r1 - s2
val r1 = r1 + s2
val s2 = r2 - s1
val r2 = r2 + s1
in
update(px,i2,r1*cc1 - s2*ss1);
update(py,i2,~s2*cc1 - r1*ss1);
update(px,i3,s3*cc3 + r2*ss3);
update(py,i3,r2*cc3 - s3*ss3);
i0r := i0 + !id
end;
is := 2 * !id - !n2 + j;
id := 4 * !id
(* ; dump px py *)
end
end)
end)
end;

(************************************)
(*  Last stage, length=2 butterfly  *)
(************************************)

let val is = ref 1
val id = ref 4
in
while !is < n do
let val i0r = ref (!is)
in
while !i0r <= n do
let val i0 = !i0r
val i1 = i0 + 1
val r1 = sub(px,i0)
val _ = update(px,i0,r1 + sub(px,i1))
val _ = update(px,i1,r1 - sub(px,i1))
val r1 = sub(py,i0)
in
update(py,i0,r1 + sub(py,i1));
update(py,i1,r1 - sub(py,i1));
i0r := i0 + !id
end;
is := 2 * !id - 1;
id := 4 * !id
end
end;
(*
print "\nbutterfly\n";
dump px py;
*)
(*************************)
(*  Bit reverse counter  *)
(*************************)

let val j = ref 1
in
for (1,n-1,fn i =>
(if i < !j then (
let val xt1 = sub(px,!j)
val xt2 = sub(px,i)
val _ = update(px,!j,xt2)
val _ = update(px,i,xt1)
val yt1 = sub(py,!j)
val yt2 = sub(py,i)
in
update(py,!j,yt2);
update(py,i,yt1)
end)
else ();
let val k = ref(n div 2)
in
while !k < !j do (j := !j - !k; k := !k div 2);
j := !j + !k
end));
(*
print "\nbit reverse\n";
dump px py;
*)
n
end
end;

fun test np =
(print_int np; print_string "... ";
let val enp = real np
val npm = np div 2 - 1
val pxr = array ((np+2), 0.0)
val pxi = array ((np+2), 0.0)
val t = pi / enp
val _ = update(pxr,1,(enp - 1.0) * 0.5)
val _ = update(pxi,1,0.0)
val n2 = np div 2
val _ = update(pxr,n2+1,~0.5)
val _ = update(pxi,n2+1,0.0)
in
for (1,npm,fn i =>
let val j = np - i
val _ = update(pxr,i+1,~0.5)
val _ = update(pxr,j+1,~0.5)
val z = t * (real i)
val y = ~0.5 * (cos(z)/sin(z))
in
update(pxi,i+1,y);
update(pxi,j+1,~y)
end);
(*
print "\n"; print "before fft: \n";
dump pxr pxi;
*)
fft pxr pxi np;
(*
print "\n"; print "after fft: \n";
dump pxr pxi;
*)
let val zr = ref 0.0
val zi = ref 0.0
val kr = ref 0
val ki = ref 0
in
for (0,np-1,fn i =>
let val a = abs_real(sub(pxr,i+1) - (real i))
in
if !zr < a then
(zr := a;
kr := i)
else ();
let val a = abs_real(sub(pxi,i+1))
in
if !zi < a then
(zi := a;
ki := i)
else ()
end
end);
let val zm = if abs_real (!zr) < abs_real (!zi) then !zi else !zr
in
print_real zm; print_newline()
end
end
end)

fun doit() =
let val np = ref 16
in for(1,13,fn i => (test (!np); np := (!np)*2))
end

in
val _ = doit()
end
```