From 300f321963cf16518a5f538ae6a33a708868dac9 Mon Sep 17 00:00:00 2001
From: Ivan Ukhov <ivan.ukhov@gmail.com>
Date: Thu, 25 Jun 2015 15:40:48 -0400
Subject: [PATCH] Implement real::inverse

---
 src/complex.rs | 10 +++---
 src/lib.rs     |  8 +++--
 src/real.rs    | 83 +++++++++++++++++++++++++++++---------------------
 tests/lib.rs   |  7 +++++
 4 files changed, 66 insertions(+), 42 deletions(-)

diff --git a/src/complex.rs b/src/complex.rs
index c125ec5..e6835df 100644
--- a/src/complex.rs
+++ b/src/complex.rs
@@ -9,8 +9,7 @@ use number::c64;
 ///
 /// The number of points should be a power of two.
 pub fn forward(data: &mut [c64]) {
-    let n = data.len();
-    power_of_two!(n);
+    let n = power_of_two!(data);
     rearrange(data, n);
     perform(data, n, false);
 }
@@ -19,8 +18,7 @@ pub fn forward(data: &mut [c64]) {
 ///
 /// The number of points should be a power of two.
 pub fn inverse(data: &mut [c64], scaling: bool) {
-    let n = data.len();
-    power_of_two!(n);
+    let n = power_of_two!(data);
     rearrange(data, n);
     perform(data, n, true);
     if scaling {
@@ -46,12 +44,12 @@ fn rearrange(data: &mut [c64], n: usize) {
 fn perform(data: &mut [c64], n: usize, inverse: bool) {
     use std::f64::consts::PI;
 
-    let pi = if inverse { PI } else { -PI };
+    let sign = if inverse { 1.0 } else { -1.0 };
     let mut step = 1;
     while step < n {
         let jump = step << 1;
         let (multiplier, mut factor) = {
-            let delta = pi / step as f64;
+            let delta = sign * PI / step as f64;
             let sine = (0.5 * delta).sin();
             (c64(-2.0 * sine * sine, delta.sin()), c64(1.0, 0.0))
         };
diff --git a/src/lib.rs b/src/lib.rs
index caca504..305978d 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -7,8 +7,12 @@
 extern crate complex as number;
 
 macro_rules! power_of_two(
-    ($n:expr) => (if $n < 1 || $n & ($n - 1) != 0 {
-        panic!("expected the number of points to be a power of two");
+    ($data:expr) => ({
+        let n = $data.len();
+        if n < 1 || n & (n - 1) != 0 {
+            panic!("expected the number of points to be a power of two");
+        }
+        n
     });
 );
 
diff --git a/src/real.rs b/src/real.rs
index aeae779..073220d 100644
--- a/src/real.rs
+++ b/src/real.rs
@@ -2,6 +2,14 @@
 
 use number::{Complex, c64};
 
+macro_rules! reinterpret(
+    ($data:ident) => (unsafe {
+        use std::slice::from_raw_parts_mut;
+        let n = power_of_two!($data);
+        (from_raw_parts_mut($data.as_mut_ptr() as *mut _, n / 2), n / 2)
+    });
+);
+
 /// Perform the forward transform.
 ///
 /// The number of points should be a power of two. The data are replaced by the
@@ -15,46 +23,24 @@ use number::{Complex, c64};
 ///    Flannery, “Numerical Recipes 3rd Edition: The Art of Scientific
 ///    Computing,” Cambridge University Press, 2007.
 pub fn forward(data: &mut [f64]) {
-    use std::f64::consts::PI;
-    use std::slice::from_raw_parts_mut;
-
-    let (data, n) = unsafe {
-        let n = data.len();
-        power_of_two!(n);
-        (from_raw_parts_mut(data.as_mut_ptr() as *mut _, n / 2), n / 2)
-    };
-
+    let (data, n) = reinterpret!(data);
     ::complex::forward(data);
+    compose(data, n, false);
+}
 
-    let (multiplier, mut factor) = {
-        let delta = -PI / n as f64;
-        let sine = (0.5 * delta).sin();
-        (c64(-2.0 * sine * sine, delta.sin()), c64(0.0, 1.0))
-    };
-
-    data[0] = c64(
-        data[0].re() + data[0].im(),
-        data[0].re() - data[0].im(),
-    );
-
-    for i in 1..(n / 2) {
-        let j = n - i;
-
-        factor = multiplier * factor + factor;
-        let part1 = (data[i] + data[j].conj()) * 0.5;
-        let part2 = -(data[i] - data[j].conj()) * 0.5;
-
-        data[i] = part1 + factor * part2;
-        data[j] = (part1 - factor * part2).conj();
-    }
-
-    data[n / 2] = data[n / 2].conj();
+/// Perform the inverse transform.
+///
+/// The number of points should be a power of two. The data should be packed as
+/// described in `real::forward`.
+pub fn inverse(data: &mut [f64], scaling: bool) {
+    let (data, n) = reinterpret!(data);
+    compose(data, n, true);
+    ::complex::inverse(data, scaling);
 }
 
 /// Unpack a compressed representation produced by `real::forward`.
 pub fn unpack(data: &[f64]) -> Vec<c64> {
-    let n = data.len();
-    power_of_two!(n);
+    let n = power_of_two!(data);
 
     let mut cdata = Vec::with_capacity(n);
     unsafe { cdata.set_len(n) };
@@ -71,6 +57,35 @@ pub fn unpack(data: &[f64]) -> Vec<c64> {
     cdata
 }
 
+fn compose(data: &mut [c64], n: usize, inverse: bool) {
+    use std::f64::consts::PI;
+
+    data[0] = c64(data[0].re() + data[0].im(), data[0].re() - data[0].im());
+    if inverse {
+        data[0] = data[0] * 0.5;
+    }
+
+    let (multiplier, mut factor) = {
+        let sign = if inverse { 1.0 } else { -1.0 };
+        let delta = sign * PI / n as f64;
+        let sine = (0.5 * delta).sin();
+        (c64(-2.0 * sine * sine, delta.sin()), c64(0.0, -sign))
+    };
+
+    for i in 1..(n / 2) {
+        let j = n - i;
+
+        let part1 = (data[i] + data[j].conj()) * 0.5;
+        let part2 = -(data[i] - data[j].conj()) * 0.5;
+
+        factor = multiplier * factor + factor;
+        data[i] = part1 + factor * part2;
+        data[j] = (part1 - factor * part2).conj();
+    }
+
+    data[n / 2] = data[n / 2].conj();
+}
+
 #[cfg(test)]
 mod tests {
     use number::c64;
diff --git a/tests/lib.rs b/tests/lib.rs
index ab89000..9f6fe3f 100644
--- a/tests/lib.rs
+++ b/tests/lib.rs
@@ -36,6 +36,13 @@ fn real_forward() {
     }
 }
 
+#[test]
+fn real_inverse() {
+    let mut data = fixtures::FREQUENCY_DATA_FOR_REAL.to_vec();
+    fft::real::inverse(&mut data, true);
+    assert::close(&data, &fixtures::TIME_DATA[..], 1e-14);
+}
+
 fn as_f64<'l>(slice: &'l [c64]) -> &'l [f64] {
     unsafe {
         std::slice::from_raw_parts(slice.as_ptr() as *const _, 2 * slice.len())