Robust/src/Benchmarks/Prefetch/2DFFT/java/FFT1d.java

   1 //Title:        1-d mixed radix FFT.
   2 //Version:
   3 //Copyright:    Copyright (c) 1998
   4 //Author:       Dongyan Wang
   5 //Company:      University of Wisconsin-Milwaukee.
   6 //Description:
   7 //  The number of DFT is factorized.
   8 //
   9 // Some short FFTs, such as length 2, 3, 4, 5, 8, 10, are used
  10 // to improve the speed.
  11 //
  12 // Prime factors are processed using DFT. In the future, we can
  13 // improve this part.
  14 // Note: there is no limit how large the prime factor can be,
  15 // because for a set of data of an image, the length can be
  16 // random, ie. an image can have size 263 x 300, where 263 is
  17 // a large prime factor.
  18 //
  19 // A permute() function is used to make sure FFT can be calculated
  20 // in place.
  21 //
  22 // A triddle() function is used to perform the FFT.
  23 //
  24 // This program is for FFT of complex data, if the input is real,
  25 // the program can be further improved. Because I want to use the
  26 // same program to do IFFT, whose input is often complex, so I
  27 // still use this program.
  28 //
  29 // To save the memory and improve the speed, double data are used
  30 // instead of double, but I do have a double version transforms.fft.
  31 //
  32 // Factorize() is done in constructor, transforms.fft() is needed to be
  33 // called to do FFT, this is good for use in fft2d, then
  34 // factorize() is not needed for each row/column of data, since
  35 // each row/column of a matrix has the same length.
  36 //
  37
  38
  39 public class FFT1d {
  40   // Maximum numbers of factors allowed.
  41   //private static final int MaxFactorsNumber = 30;
  42   private static final int MaxFactorsNumber = 37;
  43
  44   // cos2to3PI = cos(2*pi/3), using for 3 point FFT.
  45   // cos(2*PI/3) is not -1.5
  46   private static final double cos2to3PI = -1.5000f;
  47   // sin2to3PI = sin(2*pi/3), using for 3 point FFT.
  48   private static final double sin2to3PI = 8.6602540378444E-01f;
  49
  50   // TwotoFivePI   = 2*pi/5.
  51   // c51, c52, c53, c54, c55 are used in fft5().
  52   // c51 =(cos(TwotoFivePI)+cos(2*TwotoFivePI))/2-1.
  53   private static final double c51 = -1.25f;
  54   // c52 =(cos(TwotoFivePI)-cos(2*TwotoFivePI))/2.
  55   private static final double c52 = 5.5901699437495E-01f;
  56   // c53 = -sin(TwotoFivePI).
  57   private static final double c53 = -9.5105651629515E-01f;
  58   // c54 =-(sin(TwotoFivePI)+sin(2*TwotoFivePI)).
  59   private static final double c54 = -1.5388417685876E+00f;
  60   // c55 =(sin(TwotoFivePI)-sin(2*TwotoFivePI)).
  61   private static final double c55 = 3.6327126400268E-01f;
  62
  63   // OnetoSqrt2 = 1/sqrt(2), used in fft8().
  64   private static final double OnetoSqrt2 = 7.0710678118655E-01f;
  65
  66   private static int lastRadix = 0;
  67
  68   int N;              // length of N point FFT.
  69   int NumofFactors;   // Number of factors of N.
  70   static final int maxFactor = 20;      // Maximum factor of N.
  71
  72   int factors[];      // Factors of N processed in the current stage.
  73   int sofar[];        // Finished factors before the current stage.
  74   int remain[];       // Finished factors after the current stage.
  75
  76   double inputRe[],  inputIm[];   // Input  of FFT.
  77   double temRe[],    temIm[];     // Intermediate result of FFT.
  78   double outputRe[], outputIm[];  // Output of FFT.
  79   //static boolean factorsWerePrinted = false;
  80   boolean factorsWerePrinted = false;
  81
  82   // Constructor: FFT of Complex data.
  83   public FFT1d(int N) {
  84     this.N = N;
  85     outputRe = new double[N];
  86     outputIm = new double[N];
  87
  88     factorize();
  89     //printFactors();
  90
  91     // Allocate memory for intermediate result of FFT.
  92     temRe = new double[maxFactor];
  93     temIm = new double[maxFactor];
  94   }
  95
  96   public void fft(double inputRe[], double inputIm[]) {
  97     // First make sure inputRe & inputIm are of the same length.
  98     if (inputRe.length != N || inputIm.length != N) {
  99       System.out.println("Error: the length of real part & imaginary part " +
 100           "of the input to 1-d FFT are different");
 101       return;
 102     } else {
 103       this.inputRe = inputRe;
 104       this.inputIm = inputIm;
 105
 106       permute();
 107       //System.out.println("ready to twiddle");
 108
 109       for (int factorIndex = 0; factorIndex < NumofFactors; factorIndex++)
 110         twiddle(factorIndex);
 111       //System.out.println("ready to copy");
 112
 113       // Copy the output[] data to input[], so the output can be
 114       // returned in the input array.
 115       for (int i = 0; i < N; i++) {
 116         inputRe[i] = outputRe[i];
 117         inputIm[i] = outputIm[i];
 118       }
 119
 120     }
 121   }
 122
 123   public void printFactors() {
 124     if (factorsWerePrinted) return;
 125     factorsWerePrinted = true;
 126     //System.out.println("factors.length = " + factors.length + "\n");
 127     for (int i = 0; i < factors.length; i++)
 128       System.out.println("factors[i] = " + factors[i]);
 129   }
 130
 131   private void factorize() {
 132     int radices[] = {2, 3, 4, 5, 8, 10};
 133     int temFactors[] = new int[MaxFactorsNumber];
 134
 135     // 1 - point FFT, no need to factorize N.
 136     if (N == 1) {
 137       temFactors[0] = 1;
 138       NumofFactors = 1;
 139     }
 140
 141     // N - point FFT, N is needed to be factorized.
 142     int n = N;
 143     int index = 0;    // index of temFactors.
 144     int i = radices.length - 1;
 145
 146     while ((n > 1) && (i >= 0)) {
 147       if ((n % radices[i]) == 0) {
 148         n /= radices[i];
 149         temFactors[index++] = radices[i];
 150       } else
 151         i--;
 152     }
 153
 154     // Substitute 2x8 with 4x4.
 155     // index>0, in the case only one prime factor, such as N=263.
 156     if ((index > 0) && (temFactors[index - 1] == 2))
 157       for (i = index - 2; i >= 0; i--)
 158         if (temFactors[i] == 8) {
 159           temFactors[index - 1] = temFactors[i] = 4;
 160           // break out of for loop, because only one '2' will exist in
 161           // temFactors, so only one substitutation is needed.
 162           break;
 163         }
 164
 165     if (n > 1) {
 166       for (int k = 2; k < Math.sqrt(n) + 1; k++)
 167         while ((n % k) == 0) {
 168           n /= k;
 169           temFactors[index++] = k;
 170         }
 171       if (n > 1) {
 172         temFactors[index++] = n;
 173       }
 174     }
 175     NumofFactors = index;
 176     /*
 177     if(temFactors[NumofFactors-1] > 10)
 178        maxFactor = n;
 179     else
 180        maxFactor = 10;
 181        */
 182
 183     // Inverse temFactors and store factors into factors[].
 184     factors = new int[NumofFactors];
 185     for (i = 0; i < NumofFactors; i++) {
 186       factors[i] = temFactors[NumofFactors - i - 1];
 187     }
 188
 189     // Calculate sofar[], remain[].
 190     // sofar[]  : finished factors before the current stage.
 191     // factors[]: factors of N processed in the current stage.
 192     // remain[] : finished factors after the current stage.
 193     sofar = new int[NumofFactors];
 194     remain = new int[NumofFactors];
 195
 196     remain[0] = N / factors[0];
 197     sofar[0] = 1;
 198     for (i = 1; i < NumofFactors; i++) {
 199       sofar[i] = sofar[i - 1] * factors[i - 1];
 200       remain[i] = remain[i - 1] / factors[i];
 201     }
 202   }   // End of function factorize().
 203
 204   private void permute() {
 205     int count[] = new int[MaxFactorsNumber];
 206     int j;
 207     int k = 0;
 208
 209     for (int i = 0; i < N - 1; i++) {
 210       outputRe[i] = inputRe[k];
 211       outputIm[i] = inputIm[k];
 212       j = 0;
 213       k = k + remain[j];
 214       count[0] = count[0] + 1;
 215       while (count[j] >= factors[j]) {
 216         count[j] = 0;
 217         k = k - (j == 0?N:remain[j - 1]) + remain[j + 1];
 218         j++;
 219         count[j] = count[j] + 1;
 220       }
 221     }
 222     outputRe[N - 1] = inputRe[N - 1];
 223     outputIm[N - 1] = inputIm[N - 1];
 224   }   // End of function permute().
 225
 226   private void twiddle(int factorIndex) {
 227     // Get factor data.
 228     int sofarRadix = sofar[factorIndex];
 229     int radix = factors[factorIndex];
 230     int remainRadix = remain[factorIndex];
 231
 232     double tem;   // Temporary variable to do data exchange.
 233
 234     double W = 2 * (double) Math.PI / (sofarRadix * radix);
 235     double cosW = (double) Math.cos(W);
 236     double sinW = -(double) Math.sin(W);
 237
 238     double twiddleRe[] = new double[radix];
 239     double twiddleIm[] = new double[radix];
 240     double twRe = 1.0f, twIm = 0f;
 241
 242     //Initialize twiddle addBk.address variables.
 243     int dataOffset = 0, groupOffset = 0, address = 0;
 244
 245     for (int dataNo = 0; dataNo < sofarRadix; dataNo++) {
 246       //System.out.println("datano="+dataNo);
 247       if (sofarRadix > 1) {
 248         twiddleRe[0] = 1.0f;
 249         twiddleIm[0] = 0.0f;
 250         twiddleRe[1] = twRe;
 251         twiddleIm[1] = twIm;
 252         for (int i = 2; i < radix; i++) {
 253
 254
 255           twiddleRe[i] = twRe * twiddleRe[i - 1] - twIm * twiddleIm[i - 1];
 256           twiddleIm[i] = twIm * twiddleRe[i - 1] + twRe * twiddleIm[i - 1];
 257         }
 258         tem = cosW * twRe - sinW * twIm;
 259         twIm = sinW * twRe + cosW * twIm;
 260         twRe = tem;
 261       }
 262       for (int groupNo = 0; groupNo < remainRadix; groupNo++) {
 263         //System.out.println("groupNo="+groupNo);
 264         if ((sofarRadix > 1) && (dataNo > 0)) {
 265           temRe[0] = outputRe[address];
 266           temIm[0] = outputIm[address];
 267           int blockIndex = 1;
 268           do {
 269             address = address + sofarRadix;
 270             temRe[blockIndex] = twiddleRe[blockIndex] * outputRe[address] -
 271               twiddleIm[blockIndex] * outputIm[address];
 272             temIm[blockIndex] = twiddleRe[blockIndex] * outputIm[address] +
 273               twiddleIm[blockIndex] * outputRe[address];
 274             blockIndex++;
 275           } while (blockIndex < radix);
 276         } else
 277           for (int i = 0; i < radix; i++) {
 278             //System.out.println("temRe.length="+temRe.length);
 279             //System.out.println("i = "+i);
 280             temRe[i] = outputRe[address];
 281             temIm[i] = outputIm[address];
 282             address += sofarRadix;
 283           }
 284         //System.out.println("radix="+radix);
 285         switch (radix) {
 286           case 2:
 287             tem = temRe[0] + temRe[1];
 288             temRe[1] = temRe[0] - temRe[1];
 289             temRe[0] = tem;
 290             tem = temIm[0] + temIm[1];
 291             temIm[1] = temIm[0] - temIm[1];
 292             temIm[0] = tem;
 293             break;
 294           case 3:
 295             double t1Re = temRe[1] + temRe[2];
 296             double t1Im = temIm[1] + temIm[2];
 297             temRe[0] = temRe[0] + t1Re;
 298             temIm[0] = temIm[0] + t1Im;
 299
 300             double m1Re = cos2to3PI * t1Re;
 301             double m1Im = cos2to3PI * t1Im;
 302             double m2Re = sin2to3PI * (temIm[1] - temIm[2]);
 303             double m2Im = sin2to3PI * (temRe[2] - temRe[1]);
 304             double s1Re = temRe[0] + m1Re;
 305             double s1Im = temIm[0] + m1Im;
 306
 307             temRe[1] = s1Re + m2Re;
 308             temIm[1] = s1Im + m2Im;
 309             temRe[2] = s1Re - m2Re;
 310             temIm[2] = s1Im - m2Im;
 311             break;
 312           case 4:
 313             fft4(temRe, temIm);
 314             break;
 315           case 5:
 316             fft5(temRe, temIm);
 317             break;
 318           case 8:
 319             fft8();
 320             break;
 321           case 10:
 322             fft10();
 323             break;
 324           default  :
 325             fftPrime(radix);
 326             break;
 327         }
 328         address = groupOffset;
 329         for (int i = 0; i < radix; i++) {
 330           outputRe[address] = temRe[i];
 331           outputIm[address] = temIm[i];
 332           address += sofarRadix;
 333         }
 334         groupOffset += sofarRadix * radix;
 335         address = groupOffset;
 336       }
 337       groupOffset = ++dataOffset;
 338       address = groupOffset;
 339     }
 340   } // End of function twiddle().
 341
 342   // The two arguments dataRe[], dataIm[] are mainly for using in fft8();
 343   private void fft4(double dataRe[], double dataIm[]) {
 344     double t1Re,t1Im, t2Re,t2Im;
 345     double m2Re,m2Im, m3Re,m3Im;
 346
 347     t1Re = dataRe[0] + dataRe[2];
 348     t1Im = dataIm[0] + dataIm[2];
 349     t2Re = dataRe[1] + dataRe[3];
 350     t2Im = dataIm[1] + dataIm[3];
 351
 352     m2Re = dataRe[0] - dataRe[2];
 353     m2Im = dataIm[0] - dataIm[2];
 354     m3Re = dataIm[1] - dataIm[3];
 355     m3Im = dataRe[3] - dataRe[1];
 356
 357     dataRe[0] = t1Re + t2Re;
 358     dataIm[0] = t1Im + t2Im;
 359     dataRe[2] = t1Re - t2Re;
 360     dataIm[2] = t1Im - t2Im;
 361     dataRe[1] = m2Re + m3Re;
 362     dataIm[1] = m2Im + m3Im;
 363     dataRe[3] = m2Re - m3Re;
 364     dataIm[3] = m2Im - m3Im;
 365   }   // End of function fft4().
 366
 367   // The two arguments dataRe[], dataIm[] are mainly for using in fft10();
 368   private void fft5(double dataRe[], double dataIm[]) {
 369     double t1Re,t1Im, t2Re,t2Im, t3Re,t3Im, t4Re,t4Im, t5Re,t5Im;
 370     double m1Re,m1Im, m2Re,m2Im, m3Re,m3Im, m4Re,m4Im, m5Re,m5Im;
 371     double s1Re,s1Im, s2Re,s2Im, s3Re,s3Im, s4Re,s4Im, s5Re,s5Im;
 372
 373     t1Re = dataRe[1] + dataRe[4];
 374     t1Im = dataIm[1] + dataIm[4];
 375     t2Re = dataRe[2] + dataRe[3];
 376     t2Im = dataIm[2] + dataIm[3];
 377     t3Re = dataRe[1] - dataRe[4];
 378     t3Im = dataIm[1] - dataIm[4];
 379     t4Re = dataRe[3] - dataRe[2];
 380     t4Im = dataIm[3] - dataIm[2];
 381     t5Re = t1Re + t2Re;
 382     t5Im = t1Im + t2Im;
 383
 384     dataRe[0] = dataRe[0] + t5Re;
 385     dataIm[0] = dataIm[0] + t5Im;
 386
 387     m1Re = c51 * t5Re;
 388     m1Im = c51 * t5Im;
 389     m2Re = c52 * (t1Re - t2Re);
 390     m2Im = c52 * (t1Im - t2Im);
 391     m3Re = -c53 * (t3Im + t4Im);
 392     m3Im = c53 * (t3Re + t4Re);
 393     m4Re = -c54 * t4Im;
 394     m4Im = c54 * t4Re;
 395     m5Re = -c55 * t3Im;
 396     m5Im = c55 * t3Re;
 397
 398     s3Re = m3Re - m4Re;
 399     s3Im = m3Im - m4Im;
 400     s5Re = m3Re + m5Re;
 401     s5Im = m3Im + m5Im;
 402     s1Re = dataRe[0] + m1Re;
 403     s1Im = dataIm[0] + m1Im;
 404     s2Re = s1Re + m2Re;
 405     s2Im = s1Im + m2Im;
 406     s4Re = s1Re - m2Re;
 407     s4Im = s1Im - m2Im;
 408
 409     dataRe[1] = s2Re + s3Re;
 410     dataIm[1] = s2Im + s3Im;
 411     dataRe[2] = s4Re + s5Re;
 412     dataIm[2] = s4Im + s5Im;
 413     dataRe[3] = s4Re - s5Re;
 414     dataIm[3] = s4Im - s5Im;
 415     dataRe[4] = s2Re - s3Re;
 416     dataIm[4] = s2Im - s3Im;
 417   }   // End of function fft5().
 418
 419   private void fft8() {
 420     double data1Re[] = new double[4];
 421     double data1Im[] = new double[4];
 422     double data2Re[] = new double[4];
 423     double data2Im[] = new double[4];
 424     double tem;
 425
 426     // To improve the speed, use direct assaignment instead for loop here.
 427     data1Re[0] = temRe[0];
 428     data2Re[0] = temRe[1];
 429     data1Re[1] = temRe[2];
 430     data2Re[1] = temRe[3];
 431     data1Re[2] = temRe[4];
 432     data2Re[2] = temRe[5];
 433     data1Re[3] = temRe[6];
 434     data2Re[3] = temRe[7];
 435
 436     data1Im[0] = temIm[0];
 437     data2Im[0] = temIm[1];
 438     data1Im[1] = temIm[2];
 439     data2Im[1] = temIm[3];
 440     data1Im[2] = temIm[4];
 441     data2Im[2] = temIm[5];
 442     data1Im[3] = temIm[6];
 443     data2Im[3] = temIm[7];
 444
 445     fft4(data1Re, data1Im);
 446     fft4(data2Re, data2Im);
 447
 448     tem = OnetoSqrt2 * (data2Re[1] + data2Im[1]);
 449     data2Im[1] = OnetoSqrt2 * (data2Im[1] - data2Re[1]);
 450     data2Re[1] = tem;
 451     tem = data2Im[2];
 452     data2Im[2] = -data2Re[2];
 453     data2Re[2] = tem;
 454     tem = OnetoSqrt2 * (data2Im[3] - data2Re[3]);
 455     data2Im[3] = -OnetoSqrt2 * (data2Re[3] + data2Im[3]);
 456     data2Re[3] = tem;
 457
 458     temRe[0] = data1Re[0] + data2Re[0];
 459     temRe[4] = data1Re[0] - data2Re[0];
 460     temRe[1] = data1Re[1] + data2Re[1];
 461     temRe[5] = data1Re[1] - data2Re[1];
 462     temRe[2] = data1Re[2] + data2Re[2];
 463     temRe[6] = data1Re[2] - data2Re[2];
 464     temRe[3] = data1Re[3] + data2Re[3];
 465     temRe[7] = data1Re[3] - data2Re[3];
 466
 467     temIm[0] = data1Im[0] + data2Im[0];
 468     temIm[4] = data1Im[0] - data2Im[0];
 469     temIm[1] = data1Im[1] + data2Im[1];
 470     temIm[5] = data1Im[1] - data2Im[1];
 471     temIm[2] = data1Im[2] + data2Im[2];
 472     temIm[6] = data1Im[2] - data2Im[2];
 473     temIm[3] = data1Im[3] + data2Im[3];
 474     temIm[7] = data1Im[3] - data2Im[3];
 475   }   // End of function fft8().
 476
 477   private void fft10() {
 478     double data1Re[] = new double[5];
 479     double data1Im[] = new double[5];
 480     double data2Re[] = new double[5];
 481     double data2Im[] = new double[5];
 482
 483     // To improve the speed, use direct assaignment instead for loop here.
 484     data1Re[0] = temRe[0];
 485     data2Re[0] = temRe[5];
 486     data1Re[1] = temRe[2];
 487     data2Re[1] = temRe[7];
 488     data1Re[2] = temRe[4];
 489     data2Re[2] = temRe[9];
 490     data1Re[3] = temRe[6];
 491     data2Re[3] = temRe[1];
 492     data1Re[4] = temRe[8];
 493     data2Re[4] = temRe[3];
 494     data1Im[0] = temIm[0];
 495     data2Im[0] = temIm[5];
 496     data1Im[1] = temIm[2];
 497     data2Im[1] = temIm[7];
 498     data1Im[2] = temIm[4];
 499     data2Im[2] = temIm[9];
 500     data1Im[3] = temIm[6];
 501     data2Im[3] = temIm[1];
 502     data1Im[4] = temIm[8];
 503     data2Im[4] = temIm[3];
 504
 505     fft5(data1Re, data1Im);
 506     fft5(data2Re, data2Im);
 507
 508     temRe[0] = data1Re[0] + data2Re[0];
 509     temRe[5] = data1Re[0] - data2Re[0];
 510     temRe[6] = data1Re[1] + data2Re[1];
 511     temRe[1] = data1Re[1] - data2Re[1];
 512     temRe[2] = data1Re[2] + data2Re[2];
 513     temRe[7] = data1Re[2] - data2Re[2];
 514     temRe[8] = data1Re[3] + data2Re[3];
 515     temRe[3] = data1Re[3] - data2Re[3];
 516     temRe[4] = data1Re[4] + data2Re[4];
 517     temRe[9] = data1Re[4] - data2Re[4];
 518
 519     temIm[0] = data1Im[0] + data2Im[0];
 520     temIm[5] = data1Im[0] - data2Im[0];
 521     temIm[6] = data1Im[1] + data2Im[1];
 522     temIm[1] = data1Im[1] - data2Im[1];
 523     temIm[2] = data1Im[2] + data2Im[2];
 524     temIm[7] = data1Im[2] - data2Im[2];
 525     temIm[8] = data1Im[3] + data2Im[3];
 526     temIm[3] = data1Im[3] - data2Im[3];
 527     temIm[4] = data1Im[4] + data2Im[4];
 528     temIm[9] = data1Im[4] - data2Im[4];
 529   }   // End of function fft10().
 530
 531   public double sqrt(double d) {
 532     return Math.sqrt(d);
 533   }
 534
 535   private void fftPrime(int radix) {
 536     // Initial WRe, WIm.
 537     double W = 2 * (double) Math.PI / radix;
 538     double cosW = (double) Math.cos(W);
 539     double sinW = -(double) Math.sin(W);
 540     double WRe[] = new double[radix];
 541     double WIm[] = new double[radix];
 542
 543     WRe[0] = 1;
 544     WIm[0] = 0;
 545     WRe[1] = cosW;
 546     WIm[1] = sinW;
 547
 548     for (int i = 2; i < radix; i++) {
 549       WRe[i] = cosW * WRe[i - 1] - sinW * WIm[i - 1];
 550       WIm[i] = sinW * WRe[i - 1] + cosW * WIm[i - 1];
 551     }
 552
 553     // FFT of prime length data, using DFT, can be improved in the future.
 554     double rere, reim, imre, imim;
 555     int j, k;
 556     int max = (radix + 1) / 2;
 557
 558     double tem1Re[] = new double[max];
 559     double tem1Im[] = new double[max];
 560     double tem2Re[] = new double[max];
 561     double tem2Im[] = new double[max];
 562
 563     for (j = 1; j < max; j++) {
 564       tem1Re[j] = temRe[j] + temRe[radix - j];
 565       tem1Im[j] = temIm[j] - temIm[radix - j];
 566       tem2Re[j] = temRe[j] - temRe[radix - j];
 567       tem2Im[j] = temIm[j] + temIm[radix - j];
 568     }
 569
 570     for (j = 1; j < max; j++) {
 571       temRe[j] = temRe[0];
 572       temIm[j] = temIm[0];
 573       temRe[radix - j] = temRe[0];
 574       temIm[radix - j] = temIm[0];
 575       k = j;
 576       for (int i = 1; i < max; i++) {
 577         rere = WRe[k] * tem1Re[i];
 578         imim = WIm[k] * tem1Im[i];
 579         reim = WRe[k] * tem2Im[i];
 580         imre = WIm[k] * tem2Re[i];
 581
 582         temRe[radix - j] += rere + imim;
 583         temIm[radix - j] += reim - imre;
 584         temRe[j] += rere - imim;
 585         temIm[j] += reim + imre;
 586
 587         k = k + j;
 588         if (k >= radix)
 589           k = k - radix;
 590       }
 591     }
 592     for (j = 1; j < max; j++) {
 593       temRe[0] = temRe[0] + tem1Re[j];
 594       temIm[0] = temIm[0] + tem2Im[j];
 595     }
 596   }   // End of function fftPrime().
 597
 598 } // End of class FFT2d