Complete FFT example.

doc/invargent-manual.tm

@@ -1142,11 +1142,12 @@
   because InvarGenT, noticing the failure, generates an OCaml source with
   more type information, as if the <verbatim|-full_annot> option was used.
 
-  The example <verbatim|liquid_fft_simpler.gadt> contains assertions, but
-  despite them is nearly as hard as <verbatim|liquid_fft.gadt>. It needs the
-  option <verbatim|-same_with_assertions> to not switch to settings tuned for
-  cases where assertions capture the harder aspects of the invariants to
-  infer.
+  The examples <verbatim|liquid_fft_simpler.gadt> and
+  <verbatim|liquid_fft_full_asserted.gadt> contain assertions, but are nearly
+  as hard as <verbatim|liquid_fft.gadt>, <verbatim|liquid_fft_full.gadt>
+  respectively. They need the option <verbatim|-same_with_assertions> to not
+  switch to settings tuned for cases where assertions capture the harder
+  aspects of the invariants to infer.
 
   Unfortunately, inference fails for some examples regardless of parameters
   setting. We discuss them in the next section.

examples/liquid_fft.gadt

@@ -64,7 +64,7 @@ let fft px py = (* n must be a power of 2! *)
       let cc3 = cos (deref a3) in
       let ss3 = sin (deref a3) in
       asgn a (plus (deref a) e);
-      asgn a3 (plus (deref a3) e);
+      asgn a3 (plus (deref a3) e3);
       let rec loop1 i0 i1 i2 i3 id =
         if n + 1 <= i3 then () else (* out_of_bounds *)
         let g_px_i0 = array_get px i0 in

examples/liquid_fft_full.gadt

@@ -0,0 +1,189 @@
+(*
+** by: Dave Edelblute, [email protected], 05 Jan 1993
+** Translated from C to de Caml: Hongwei Xi, 07 Nov 1998
+** Modified: R. Mayer to work with hist benchmark routines.
+** Translated to InvarGenT: Lukasz Stafiniak, 30 Jan 2015
+*)
+datatype Array : type * num
+external let array_make :
+  ∀n, a [0≤n]. Num n → a → Array (a, n) = "fun a b -> Array.make a b"
+external let array_get :
+  ∀n, k, a [0≤k ∧ k+1≤n]. Array (a, n) → Num k → a = "fun a b -> Array.get a b"
+external let array_set :
+  ∀n, k, a [0≤k ∧ k+1≤n]. Array (a, n) → Num k → a → () =
+  "fun a b c -> Array.set a b c"
+external let array_length :
+  ∀n, a [0≤n]. Array (a, n) → Num n = "fun a -> Array.length a"
+
+external let n2i : ∀n. Num n → Int = "fun i -> i"
+external let div2 : ∀n. Num (2 n) → Num n = "fun i -> i / 2"
+external let div4 : ∀n. Num (4 n) → Num n = "fun i -> i / 4"
+external let n2f : ∀n. Num n → Float = "float_of_int"
+external let equal : ∀a. a → a → Bool = "fun x y -> x = y"
+external let leq : ∀a. a → a → Bool = "fun x y -> x <= y"
+external let less : ∀a. a → a → Bool = "fun x y -> x < y"
+external let ignore : ∀a. a → () = "ignore"
+
+external let abs : Float → Float = "abs_float"
+external let cos : Float → Float = "cos"
+external let sin : Float → Float = "sin"
+external let neg : Float → Float = "(~-.)"
+external let minus : Float → Float → Float = "(-.)"
+external let plus : Float → Float → Float = "(+.)"
+external let mult : Float → Float → Float = "( *. )"
+external let div : Float → Float → Float = "( /. )"
+external let fl0 : Float = "0.0"
+external let fl05 : Float = "0.5"
+external let fl1 : Float = "1.0"
+external let fl2 : Float = "2.0"
+external let fl3 : Float = "3.0"
+external let fl4 : Float = "4.0"
+external let pi : Float = "4.0 *. atan 1.0"
+external let two_pi : Float = "8.0 *. atan 1.0"
+
+datatype Bounded : num * num
+datacons Index : ∀i, k, n[n ≤ i ∧ i ≤ k].Num i ⟶ Bounded (n, k)
+
+let ffor s d body =
+  let rec loop i =
+    if i <= d then (body (Index i); loop (i + 1)) else () in
+  loop s
+
+external let ref : ∀a. a → Ref a = "fun a -> ref a"
+external let asgn : ∀a. Ref a → a → () = "fun a b -> a := b"
+external let deref : ∀a. Ref a → a = "(!)"
+
+let fft px py = (* n must be a power of 2! *)
+  let n = array_length px + (-1) in
+  let rec loop n2 n4 =
+    if n2 <= 2 then () else (* the case n2 = 2 is treated below *)
+    let e = div two_pi (n2f n2) in
+    let e3 = mult fl3 e in
+    let a = ref fl0 in
+    let a3 = ref fl0 in
+    let rec forbod j' =
+    match j' with Index j ->
+    (*for j = 1 to n4 do*)
+      let cc1 = cos (deref a) in
+      let ss1 = sin (deref a) in
+      let cc3 = cos (deref a3) in
+      let ss3 = sin (deref a3) in
+      asgn a (plus (deref a) e);
+      asgn a3 (plus (deref a3) e3);
+      let rec loop1 i0 i1 i2 i3 id =
+        if n + 1 <= i3 then () else (* out_of_bounds *)
+        let g_px_i0 = array_get px i0 in
+        let g_px_i2 = array_get px i2 in
+        let r1 = minus g_px_i0 g_px_i2 in
+        let r1' = plus g_px_i0 g_px_i2 in
+        array_set px i0 r1';
+        let g_px_i1 = array_get px i1 in
+        let g_px_i3 = array_get px i3 in
+        let r2 = minus g_px_i1 g_px_i3 in
+        let r2' = plus g_px_i1 g_px_i3 in
+        array_set px i1 r2';
+        let g_py_i0 = array_get py i0 in
+        let g_py_i2 = array_get py i2 in
+        let s1 = minus g_py_i0 g_py_i2 in
+        let s1' = plus g_py_i0 g_py_i2 in
+        array_set py i0 s1';
+        let g_py_i1 = array_get py i1 in
+        let g_py_i3 = array_get py i3 in
+        let s2 = minus g_py_i1 g_py_i3 in
+        let s2' = plus g_py_i1 g_py_i3 in
+        array_set py i1 s2';
+        let s3 = minus r1 s2 in 
+        let r1 = plus r1 s2 in
+        let s2 = minus r2 s1 in
+        let r2 = plus r2 s1 in
+        array_set px i2 (minus (mult r1 cc1) (mult s2 ss1));
+        array_set py i2 (minus (mult (neg s2) cc1) (mult r1 ss1));
+        array_set px i3 (plus (mult s3 cc3) (mult r2 ss3));
+        array_set py i3 (minus (mult r2 cc3) (mult s3 ss3));
+        loop1 (i0 + id) (i1 + id) (i2 + id) (i3 + id) id in
+      let rec loop2 is id =
+        if n <= is then ()
+        else (
+          let i1 = is + n4 in
+          let i2 = i1 + n4 in
+          let i3 = i2 + n4 in
+          loop1 is i1 i2 i3 id;
+          loop2 (2 * id + j - n2) (4 * id)) in
+      loop2 j (2 * n2) in
+    ffor 1 n4 forbod; 
+    loop (div2 n2) (div2 n4) in
+  loop n (div4 n);
+
+  let rec loop1 i0 i1 id =
+    if n + 1 <= i1 then () else
+    let r1 = array_get px i0 in
+    array_set px i0 (plus r1 (array_get px i1));
+    array_set px i1 (minus r1 (array_get px i1));
+    let r1 = array_get py i0 in
+    array_set py i0 (plus r1 (array_get py i1));
+    array_set py i1 (minus r1 (array_get py i1));
+    loop1 (i0 + id) (i1 + id) id in
+  let rec loop2 is id =
+    if n <= is then () else (
+      loop1 is (is + 1) id;
+      loop2 (2 * id + (-1)) (4 * id)) in
+  loop2 1 4;
+
+  let rec loop1 j k =
+    eif j <= k then j + k
+    else loop1 (j - k) (div2 k) in
+
+  let rec loop2 i j =
+    if n <= i then () else (
+      if j <= i then () else (
+        let xt = array_get px j in
+        array_set px j (array_get px i); array_set px i (xt);
+        let xt = array_get py j in
+        array_set py j (array_get py i); array_set py i (xt));
+      let j' = loop1 j (div2 n) in
+      loop2 (i + 1) j') in
+  loop2 1 1; n
+
+let ffttest np =
+  let enp = n2f np in
+  let n2 = div2 np in
+  let npm = n2 - 1 in
+  let pxr = array_make (np+1) fl0 in
+  let pxi = array_make (np+1) fl0 in
+  let t = div pi enp in
+  array_set pxr 1 (mult (minus enp fl1) fl05);
+  array_set pxi 1 (fl0);
+  array_set pxr (n2+1) (neg (mult fl1 fl05));
+  array_set pxi (n2+1) fl0;
+  let rec forbod i =
+    if i <= npm then
+      let j = np - i in
+      array_set pxr (i+1) (neg (mult fl1 fl05));
+      array_set pxr (j+1) (neg (mult fl1 fl05));
+      let z = mult t (n2f i) in
+      let y = mult (div (cos z) (sin z)) fl05 in
+      array_set pxi (i+1) (neg y); array_set pxi (j+1) (y)
+    else () in
+  forbod 1;
+  ignore (fft pxr pxi);
+  (* lukstafi: kr and ki are placeholders? *)
+  let rec loop i zr zi kr ki =
+    if np <= i then (zr, zi) else
+      let a = abs (minus (array_get pxr (i+1)) (n2f i)) in
+      let b = less zr a in
+      let zr = if b then a else zr in
+      let kr = eif b then i else kr in
+      let a = abs (array_get pxi (i+1)) in
+      let b = less zi a in
+      let zi = if b then a else zi in
+      let ki = eif b then i else ki in
+      loop (i+1) zr zi kr ki in
+  let zr, zi = loop 0 fl0 fl0 0 0 in
+  let zm = if less (abs zr) (abs zi) then zi else zr in
+  (*in print_float zm; print_newline ()*) zm
+
+let rec loop_np i np =
+  if 17 <= i then () else
+  ( ignore (ffttest np); loop_np (i + 1) (2 * np) )
+
+let doit _ = loop_np 4 16

examples/liquid_fft_full.gadti.target

@@ -0,0 +1,83 @@
+datatype Array : type * num
+
+external val array_make : ∀n, a[0 ≤ n]. Num n → a → Array (a, n)
+
+external val array_get :
+   ∀n, k, a[0 ≤ k ∧ k + 1 ≤ n]. Array (a, n) → Num k → a
+
+external val array_set :
+   ∀n, k, a[0 ≤ k ∧ k + 1 ≤ n]. Array (a, n) → Num k → a → ()
+
+external val array_length : ∀n, a[0 ≤ n]. Array (a, n) → Num n
+
+external val n2i : ∀n. Num n → Int
+
+external val div2 : ∀n. Num (2 n) → Num n
+
+external val div4 : ∀n. Num (4 n) → Num n
+
+external val n2f : ∀n. Num n → Float
+
+external val equal : ∀a. a → a → Bool
+
+external val leq : ∀a. a → a → Bool
+
+external val less : ∀a. a → a → Bool
+
+external val ignore : ∀a. a → ()
+
+external val abs : Float → Float
+
+external val cos : Float → Float
+
+external val sin : Float → Float
+
+external val neg : Float → Float
+
+external val minus : Float → Float → Float
+
+external val plus : Float → Float → Float
+
+external val mult : Float → Float → Float
+
+external val div : Float → Float → Float
+
+external val fl0 : Float
+
+external val fl05 : Float
+
+external val fl1 : Float
+
+external val fl2 : Float
+
+external val fl3 : Float
+
+external val fl4 : Float
+
+external val pi : Float
+
+external val two_pi : Float
+
+datatype Bounded : num * num
+
+datacons Index : ∀i, k, n[n ≤ i ∧ i ≤ k].Num i ⟶ Bounded (n, k)
+
+val ffor :
+   ∀i, j, k, n[j ≤ i ∧ k ≤ n].
+   Num n → Num j → (Bounded (k, i) → ()) → ()
+
+external val ref : ∀a. a → Ref a
+
+external val asgn : ∀a. Ref a → a → ()
+
+external val deref : ∀a. Ref a → a
+
+val fft :
+   ∀k, n[1 ≤ n ∧ n + 1 ≤ k].
+   Array (Float, n + 1) → Array (Float, k) → Num n
+
+val ffttest : ∀n[2 ≤ n]. Num n → Float
+
+val loop_np : ∀k, n[2 ≤ n]. Num k → Num n → ()
+
+val doit : ∀a. a → ()