
/*
 * http://local.wasp.uwa.edu.au/~pbourke/miscellaneous/dft/
   Modification of Paul Bourkes FFT code by Peter Cusack
   to utilise the Microsoft complex type.

   This computes an in-place complex-to-complex FFT
   x and y are the real and imaginary arrays of 2^m points.
   dir =  1 gives forward transform
   dir = -1 gives reverse transform
*/

import Math.{sqrt, pow}
import scala.util.Using

/** Test that specialization handles tuples. Perform FFT transformation
 *  using pairs to represent complex numbers.
 */
object Test  {
  type Complex = (Double, Double)

  def swap(x: Array[Complex], i: Int, j: Int): Unit = {
    val tmp = x(i)
    x(i) = x(j)
    x(j) = tmp
  }

  def times(x: Complex, y: Complex): Complex =
    (x._1 * y._1 - x._2 * y._2, x._1 * y._2 + x._2 * y._1)

  def div(x: Complex, y: Complex): Complex = {
    val num = pow(y._1, 2) + pow(y._2, 2)
    ((x._1 * y._1 + x._2 * y._2)/num,
     (x._2 * y._1 - x._1 * y._2)/num)
  }

  def div(x: Complex, y: Long) =
    (x._1 / y, x._2 / y)

  def add(x: Complex, y: Complex) =
    (x._1 + y._1, x._2 + y._2)

  def minus(x: Complex, y: Complex) =
    (x._1 - y._1, x._2 - y._2)

  def FFT(dir: Int, m: Long, x: Array[(Double, Double)]): Unit = {
    var i, i1, i2,j, k, l, l1, l2, n = 0L
//   complex <double> tx, t1, u, c;
    var tx, t1, u, c = (0.0, 0.0)

   /*Calculate the number of points */
   n = 1
   for (i <- 0L until m)
      n <<= 1

   /* Do the bit reversal */
   i2 = n >> 1
   j = 0

   for (i <- 0L until (n - 1)) {
      if (i < j)
         swap(x, i.toInt, j.toInt);

      k = i2;

      while (k <= j) {
         j -= k;
         k >>= 1;
      }

      j += k;
   }

   /* Compute the FFT */
   // c.real(-1.0);
   // c.imag(0.0);
   c = (-1.0, 0.0)
   l2 = 1
   for (l <- 0L until m) {
     l1 = l2
     l2 <<= 1;
      // u.real(1.0);
      // u.imag(0.0);
     u = (1.0, 0.0)

     for (j <- 0L until l1) {
       for (i <- j.until(n, l2)) {
         i1 = i + l1;
         t1 = times(u, x(i1.toInt))
         x(i1.toInt) = minus(x(i.toInt), t1)
         x(i.toInt) = add(x(i.toInt), t1)
       }

       u = times(u,  c)
    }

     // c.imag(sqrt((1.0 - c.real()) / 2.0));
     c = (c._1, sqrt( (1.0 - c._1) / 2.0 ))
     // if (dir == 1)
     //    c.imag(-c.imag());
     if (dir == 1)
       c = (c._1, -c._2)

      // c.real(sqrt((1.0 + c.real()) / 2.0));
      c = (sqrt( (1.0 + c._1) / 2.0), c._2)
   }

   /* Scaling for forward transform */
   if (dir == 1) {
     for (i <- 0L until n)
       x(i.toInt) = div(x(i.toInt), n)
   }
  }

  def run(): Unit = {
    FFT(1, 16, data)
  }
  var data: Array[Complex] = null

  def inputFileName = {
    val cwd = System.getProperty("partest.cwd")
    if (cwd ne null) cwd + java.io.File.separator + "input2.txt"
    else "input2.txt"
  }

  def setUp(): Unit = {
    //print("Loading from %s.. ".format(inputFileName))
    def readData() = Using(io.Source.fromFile(inputFileName)) { f =>
      val lines = f.getLines()
      val n = lines.next().toInt
      val data = new Array[Complex](n)
      var i = 0
      for (line <- lines if line != "") {
        val pair = line.trim.split(" ")
        data(i) = (pair(0).trim.toDouble, pair(1).trim.toDouble)
        i += 1
      }
      data
    }
    data = readData().getOrElse(???)
    //println("[loaded]")
    println(s"Processing ${data.length} items")
  }

  def main(args: Array[String]): Unit = {
    setUp()
    run()

    println("Boxed doubles: " + runtime.BoxesRunTime.doubleBoxCount)
    println("Boxed ints: " + runtime.BoxesRunTime.integerBoxCount)
    println("Boxed longs: " + runtime.BoxesRunTime.longBoxCount)
  }
}
