gt: 	j = j - k ;
context: /**************************************************************************
  *                                                                         *
  *             Java Grande Forum Benchmark Suite - Version 2.0             *
  *                                                                         *
  *                            produced by                                  *
  *                                                                         *
  *                  Java Grande Benchmarking Project                       *
  *                                                                         *
  *                                at                                       *
  *                                                                         *
  *                Edinburgh Parallel Computing Centre                      *
  *                                                                         *
  *                email: epcc-javagrande@epcc.ed.ac.uk                     *
  *                                                                         *
  *      adapted from SciMark 2.0, author Roldan Pozo (pozo@cam.nist.gov)   *
  *             see below for previous history of this code                 *
  *                                                                         *
  *      This version copyright (c) The University of Edinburgh, 1999.      *
  *                         All rights reserved.                            *
  *                                                                         *
  **************************************************************************/
  
  
  package fft;
  import jgfutil.*; 
  
  /** Computes FFT's of complex, double precision data where n is an integer power of 2.
    * This appears to be slower than the Radix2 method,
    * but the code is smaller and simpler, and it requires no extra storage.
    * <P>
    *
    * @author Bruce R. Miller bruce.miller@nist.gov,
    * @author Derived from GSL (Gnu Scientific Library), 
    * @author GSL's FFT Code by Brian Gough bjg@vvv.lanl.gov
    */
  
    /* See {@link ComplexDoubleFFT ComplexDoubleFFT} for details of data layout.
     */
  
  public class FFT {
  
    public static double JDKtotal = 0.0;
    public static double JDKtotali = 0.0;
  
    /** Compute Fast Fourier Transform of (complex) data, in place.*/
    public static void transform (double data[]) {
      int JDKrange;
  
      JGFInstrumentor.startTimer("Section2:FFT:Kernel");
  
      transform_internal(data, -1); 
  
      JGFInstrumentor.stopTimer("Section2:FFT:Kernel");
  
      JDKrange = data.length;
      for (int i=0; i<JDKrange; i++) {
       JDKtotal += data[i];
      }
  
    }
  
    /** Compute Inverse Fast Fourier Transform of (complex) data, in place.*/
    public static void inverse (double data[]) {
  
      JGFInstrumentor.startTimer("Section2:FFT:Kernel");
  
      transform_internal(data, +1);  
      // Normalize
      int nd=data.length;
      int n =nd/2;
      double norm=1/((double) n);
      for(int i=0; i<nd; i++)
        data[i] *= norm;
  
      JGFInstrumentor.stopTimer("Section2:FFT:Kernel");
  
      for(int i=0; i<nd; i++) {
        JDKtotali += data[i];
      }
  
    }
  
    /** Accuracy check on FFT of data. Make a copy of data, Compute the FFT, then
      * the inverse and compare to the original.  Returns the rms difference.*/
    public static double test(double data[]){
      int nd = data.length;
      // Make duplicate for comparison
      double copy[] = new double[nd];
      System.arraycopy(data,0,copy,0,nd);
      // Transform & invert
      transform(data);
      inverse(data);
      // Compute RMS difference.
      double diff = 0.0;
      for(int i=0; i<nd; i++) {
        double d = data[i]-copy[i];
        diff += d*d; }
      return Math.sqrt(diff/nd); }
  
    /** Make a random array of n (complex) elements. */
    public static double[] makeRandom(int n){
      int nd = 2*n;
      double data[] = new double[nd];
      for(int i=0; i<nd; i++)
        data[i]= Math.random();
      return data; }
  
    /** Simple Test routine. */
    public static void main(String args[]){
      if (args.length == 0) {
        int n = 1024;
        System.out.println("n="+n+" => RMS Error="+test(makeRandom(n))); }
      for(int i=0; i<args.length; i++) {
        int n = Integer.parseInt(args[i]);
        System.out.println("n="+n+" => RMS Error="+test(makeRandom(n))); }
    }
    /* ______________________________________________________________________ */
  
    protected static int log2 (int n){
      int log = 0;
      for(int k=1; k < n; k *= 2, log++);
      if (n != (1 << log))
        throw new Error("FFT: Data length is not a power of 2!: "+n);
      return log; }
  
    protected static void transform_internal (double data[], int direction) {
      int n = data.length/2;
      if (n == 1) return;		// Identity operation!
      int logn = log2(n);
  
      /* bit reverse the input data for decimation in time algorithm */
      bitreverse(data) ;
  
      /* apply fft recursion */
      for (int bit = 0, dual = 1; bit < logn; bit++, dual *= 2) {
        double w_real = 1.0;
        double w_imag = 0.0;
  
        double theta = 2.0 * direction * Math.PI / (2.0 * (double) dual);
        double s = Math.sin(theta);
        double t = Math.sin(theta / 2.0);
        double s2 = 2.0 * t * t;
  
        /* a = 0 */
        for (int b = 0; b < n; b += 2 * dual) {
  	int i = 2*b ;
  	int j = 2*(b + dual);
  
  	double wd_real = data[j] ;
  	double wd_imag = data[j+1] ;
  	  
  	data[j]   = data[i]   - wd_real;
  	data[j+1] = data[i+1] - wd_imag;
  	data[i]  += wd_real;
  	data[i+1]+= wd_imag;
        }
        
        /* a = 1 .. (dual-1) */
        for (int a = 1; a < dual; a++) {
  	/* trignometric recurrence for w-> exp(i theta) w */
  	{
  	  double tmp_real = w_real - s * w_imag - s2 * w_real;
  	  double tmp_imag = w_imag + s * w_real - s2 * w_imag;
  	  w_real = tmp_real;
  	  w_imag = tmp_imag;
  	}
  	for (int b = 0; b < n; b += 2 * dual) {
  	  int i = 2*(b + a);
  	  int j = 2*(b + a + dual);
  
  	  double z1_real = data[j];
  	  double z1_imag = data[j+1];
  	      
  	  double wd_real = w_real * z1_real - w_imag * z1_imag;
  	  double wd_imag = w_real * z1_imag + w_imag * z1_real;
  
  	  data[j]   = data[i]   - wd_real;
  	  data[j+1] = data[i+1] - wd_imag;
  	  data[i]  += wd_real;
  	  data[i+1]+= wd_imag;
  	}
        }
      }
    }
  
  
    protected static void bitreverse(double data[]) {
      /* This is the Goldrader bit-reversal algorithm */
      int n=data.length/2;
      for (int i = 0, j=0; i < n - 1; i++) {
        int ii = 2*i;
        int jj = 2*j;
        int k = n / 2 ;
        if (i < j) {
  	double tmp_real    = data[ii];
  	double tmp_imag    = data[ii+1];
  	data[ii]   = data[jj];
  	data[ii+1] = data[jj+1];
  	data[jj]   = tmp_real;
  	data[jj+1] = tmp_imag; }
  
        while (k <= j) {
  