001/* ***** BEGIN LICENSE BLOCK *****
002 * Version: MPL 1.1/GPL 2.0/LGPL 2.1
003 *
004 * The contents of this file are subject to the Mozilla Public License Version
005 * 1.1 (the "License"); you may not use this file except in compliance with
006 * the License. You may obtain a copy of the License at
007 * http://www.mozilla.org/MPL/
008 *
009 * Software distributed under the License is distributed on an "AS IS" basis,
010 * WITHOUT WARRANTY OF ANY KIND, either express or implied. See the License
011 * for the specific language governing rights and limitations under the
012 * License.
013 *
014 * The Original Code is JTransforms.
015 *
016 * The Initial Developer of the Original Code is
017 * Piotr Wendykier, Emory University.
018 * Portions created by the Initial Developer are Copyright (C) 2007-2009
019 * the Initial Developer. All Rights Reserved.
020 *
021 * Alternatively, the contents of this file may be used under the terms of
022 * either the GNU General Public License Version 2 or later (the "GPL"), or
023 * the GNU Lesser General Public License Version 2.1 or later (the "LGPL"),
024 * in which case the provisions of the GPL or the LGPL are applicable instead
025 * of those above. If you wish to allow use of your version of this file only
026 * under the terms of either the GPL or the LGPL, and not to allow others to
027 * use your version of this file under the terms of the MPL, indicate your
028 * decision by deleting the provisions above and replace them with the notice
029 * and other provisions required by the GPL or the LGPL. If you do not delete
030 * the provisions above, a recipient may use your version of this file under
031 * the terms of any one of the MPL, the GPL or the LGPL.
032 *
033 * ***** END LICENSE BLOCK ***** */
034
035package edu.emory.mathcs.jtransforms.dct;
036
037import java.util.concurrent.Future;
038
039import edu.emory.mathcs.jtransforms.fft.DoubleFFT_1D;
040import edu.emory.mathcs.utils.ConcurrencyUtils;
041
042/**
043 * Computes 1D Discrete Cosine Transform (DCT) of double precision data. The
044 * size of data can be an arbitrary number. This is a parallel implementation of
045 * split-radix and mixed-radix algorithms optimized for SMP systems. <br>
046 * <br>
047 * Part of the code is derived from General Purpose FFT Package written by Takuya Ooura
048 * (http://www.kurims.kyoto-u.ac.jp/~ooura/fft.html)
049 * 
050 * @author Piotr Wendykier (piotr.wendykier@gmail.com)
051 * 
052 */
053public class DoubleDCT_1D {
054
055    private int n;
056
057    private int[] ip;
058
059    private double[] w;
060
061    private int nw;
062
063    private int nc;
064
065    private boolean isPowerOfTwo = false;
066
067    private DoubleFFT_1D fft;
068
069    private static final double PI = 3.14159265358979311599796346854418516;
070
071    /**
072     * Creates new instance of DoubleDCT_1D.
073     * 
074     * @param n
075     *            size of data
076     */
077    public DoubleDCT_1D(int n) {
078        if (n < 1) {
079            throw new IllegalArgumentException("n must be greater than 0");
080        }
081        this.n = n;
082        if (ConcurrencyUtils.isPowerOf2(n)) {
083            this.isPowerOfTwo = true;
084            this.ip = new int[(int) Math.ceil(2 + (1 << (int) (Math.log(n / 2 + 0.5) / Math.log(2)) / 2))];
085            this.w = new double[n * 5 / 4];
086            nw = ip[0];
087            if (n > (nw << 2)) {
088                nw = n >> 2;
089                makewt(nw);
090            }
091            nc = ip[1];
092            if (n > nc) {
093                nc = n;
094                makect(nc, w, nw);
095            }
096        } else {
097            this.w = makect(n);
098            fft = new DoubleFFT_1D(2 * n);
099        }
100    }
101
102    /**
103     * Computes 1D forward DCT (DCT-II) leaving the result in <code>a</code>.
104     * 
105     * @param a
106     *            data to transform
107     * @param scale
108     *            if true then scaling is performed
109     */
110    public void forward(double[] a, boolean scale) {
111        forward(a, 0, scale);
112    }
113
114    /**
115     * Computes 1D forward DCT (DCT-II) leaving the result in <code>a</code>.
116     * 
117     * @param a
118     *            data to transform
119     * @param offa
120     *            index of the first element in array <code>a</code>
121     * @param scale
122     *            if true then scaling is performed
123     */
124    public void forward(final double[] a, final int offa, boolean scale) {
125        if (n == 1)
126            return;
127        if (isPowerOfTwo) {
128            double xr = a[offa + n - 1];
129            for (int j = n - 2; j >= 2; j -= 2) {
130                a[offa + j + 1] = a[offa + j] - a[offa + j - 1];
131                a[offa + j] += a[offa + j - 1];
132            }
133            a[offa + 1] = a[offa] - xr;
134            a[offa] += xr;
135            if (n > 4) {
136                rftbsub(n, a, offa, nc, w, nw);
137                cftbsub(n, a, offa, ip, nw, w);
138            } else if (n == 4) {
139                cftbsub(n, a, offa, ip, nw, w);
140            }
141            dctsub(n, a, offa, nc, w, nw);
142            if (scale) {
143                scale(Math.sqrt(2.0 / n), a, offa);
144                a[offa] = a[offa] / Math.sqrt(2.0);
145            }
146        } else {
147            final int twon = 2 * n;
148            final double[] t = new double[twon];
149            System.arraycopy(a, offa, t, 0, n);
150            int nthreads = ConcurrencyUtils.getNumberOfThreads();
151            for (int i = n; i < twon; i++) {
152                t[i] = t[twon - i - 1];
153            }
154            fft.realForward(t);
155            if ((nthreads > 1) && (n > ConcurrencyUtils.getThreadsBeginN_1D_FFT_2Threads())) {
156                nthreads = 2;
157                final int k = n / nthreads;
158                Future<?>[] futures = new Future[nthreads];
159                for (int j = 0; j < nthreads; j++) {
160                    final int firstIdx = j * k;
161                    final int lastIdx = (j == (nthreads - 1)) ? n : firstIdx + k;
162                    futures[j] = ConcurrencyUtils.submit(new Runnable() {
163                        @Override
164                                                public void run() {
165                            for (int i = firstIdx; i < lastIdx; i++) {
166                                int twoi = 2 * i;
167                                int idx = offa + i;
168                                a[idx] = w[twoi] * t[twoi] - w[twoi + 1] * t[twoi + 1];
169                            }
170
171                        }
172                    });
173                }
174                ConcurrencyUtils.waitForCompletion(futures);
175            } else {
176                for (int i = 0; i < n; i++) {
177                    int twoi = 2 * i;
178                    int idx = offa + i;
179                    a[idx] = w[twoi] * t[twoi] - w[twoi + 1] * t[twoi + 1];
180                }
181            }
182            if (scale) {
183                scale(1 / Math.sqrt(twon), a, offa);
184                a[offa] = a[offa] / Math.sqrt(2.0);
185            }
186        }
187
188    }
189
190    /**
191     * Computes 1D inverse DCT (DCT-III) leaving the result in <code>a</code>.
192     * 
193     * @param a
194     *            data to transform
195     * @param scale
196     *            if true then scaling is performed
197     */
198    public void inverse(double[] a, boolean scale) {
199        inverse(a, 0, scale);
200    }
201
202    /**
203     * Computes 1D inverse DCT (DCT-III) leaving the result in <code>a</code>.
204     * 
205     * @param a
206     *            data to transform
207     * @param offa
208     *            index of the first element in array <code>a</code>
209     * @param scale
210     *            if true then scaling is performed
211     */
212    public void inverse(final double[] a, final int offa, boolean scale) {
213        if (n == 1)
214            return;
215        if (isPowerOfTwo) {
216            double xr;
217            if (scale) {
218                scale(Math.sqrt(2.0 / n), a, offa);
219                a[offa] = a[offa] / Math.sqrt(2.0);
220            }
221            dctsub(n, a, offa, nc, w, nw);
222            if (n > 4) {
223                cftfsub(n, a, offa, ip, nw, w);
224                rftfsub(n, a, offa, nc, w, nw);
225            } else if (n == 4) {
226                cftfsub(n, a, offa, ip, nw, w);
227            }
228            xr = a[offa] - a[offa + 1];
229            a[offa] += a[offa + 1];
230            for (int j = 2; j < n; j += 2) {
231                a[offa + j - 1] = a[offa + j] - a[offa + j + 1];
232                a[offa + j] += a[offa + j + 1];
233            }
234            a[offa + n - 1] = xr;
235        } else {
236            final int twon = 2 * n;
237            if (scale) {
238                scale(Math.sqrt(twon), a, offa);
239                a[offa] = a[offa] * Math.sqrt(2.0);
240            }
241            final double[] t = new double[twon];
242            int nthreads = ConcurrencyUtils.getNumberOfThreads();
243            if ((nthreads > 1) && (n > ConcurrencyUtils.getThreadsBeginN_1D_FFT_2Threads())) {
244                nthreads = 2;
245                final int k = n / nthreads;
246                Future<?>[] futures = new Future[nthreads];
247                for (int j = 0; j < nthreads; j++) {
248                    final int firstIdx = j * k;
249                    final int lastIdx = (j == (nthreads - 1)) ? n : firstIdx + k;
250                    futures[j] = ConcurrencyUtils.submit(new Runnable() {
251                        @Override
252                                                public void run() {
253                            for (int i = firstIdx; i < lastIdx; i++) {
254                                int twoi = 2 * i;
255                                double elem = a[offa + i];
256                                t[twoi] = w[twoi] * elem;
257                                t[twoi + 1] = -w[twoi + 1] * elem;
258                            }
259                        }
260                    });
261                }
262                ConcurrencyUtils.waitForCompletion(futures);
263            } else {
264                for (int i = 0; i < n; i++) {
265                    int twoi = 2 * i;
266                    double elem = a[offa + i];
267                    t[twoi] = w[twoi] * elem;
268                    t[twoi + 1] = -w[twoi + 1] * elem;
269                }
270            }
271            fft.realInverse(t, true);
272            System.arraycopy(t, 0, a, offa, n);
273        }
274    }
275
276    /* -------- initializing routines -------- */
277
278    private double[] makect(int n) {
279        int twon = 2 * n;
280        int idx;
281        double delta = PI / twon;
282        double deltaj;
283        double[] c = new double[twon];
284        c[0] = 1;
285        for (int j = 1; j < n; j++) {
286            idx = 2 * j;
287            deltaj = delta * j;
288            c[idx] = Math.cos(deltaj);
289            c[idx + 1] = -Math.sin(deltaj);
290        }
291        return c;
292    }
293
294    private void makewt(int nw) {
295        int j, nwh, nw0, nw1;
296        double delta, wn4r, wk1r, wk1i, wk3r, wk3i;
297        double delta2, deltaj, deltaj3;
298
299        ip[0] = nw;
300        ip[1] = 1;
301        if (nw > 2) {
302            nwh = nw >> 1;
303            delta = 0.785398163397448278999490867136046290 / nwh;
304            delta2 = delta * 2;
305            wn4r = Math.cos(delta * nwh);
306            w[0] = 1;
307            w[1] = wn4r;
308            if (nwh == 4) {
309                w[2] = Math.cos(delta2);
310                w[3] = Math.sin(delta2);
311            } else if (nwh > 4) {
312                makeipt(nw);
313                w[2] = 0.5 / Math.cos(delta2);
314                w[3] = 0.5 / Math.cos(delta * 6);
315                for (j = 4; j < nwh; j += 4) {
316                    deltaj = delta * j;
317                    deltaj3 = 3 * deltaj;
318                    w[j] = Math.cos(deltaj);
319                    w[j + 1] = Math.sin(deltaj);
320                    w[j + 2] = Math.cos(deltaj3);
321                    w[j + 3] = -Math.sin(deltaj3);
322                }
323            }
324            nw0 = 0;
325            while (nwh > 2) {
326                nw1 = nw0 + nwh;
327                nwh >>= 1;
328                w[nw1] = 1;
329                w[nw1 + 1] = wn4r;
330                if (nwh == 4) {
331                    wk1r = w[nw0 + 4];
332                    wk1i = w[nw0 + 5];
333                    w[nw1 + 2] = wk1r;
334                    w[nw1 + 3] = wk1i;
335                } else if (nwh > 4) {
336                    wk1r = w[nw0 + 4];
337                    wk3r = w[nw0 + 6];
338                    w[nw1 + 2] = 0.5 / wk1r;
339                    w[nw1 + 3] = 0.5 / wk3r;
340                    for (j = 4; j < nwh; j += 4) {
341                        int idx1 = nw0 + 2 * j;
342                        int idx2 = nw1 + j;
343                        wk1r = w[idx1];
344                        wk1i = w[idx1 + 1];
345                        wk3r = w[idx1 + 2];
346                        wk3i = w[idx1 + 3];
347                        w[idx2] = wk1r;
348                        w[idx2 + 1] = wk1i;
349                        w[idx2 + 2] = wk3r;
350                        w[idx2 + 3] = wk3i;
351                    }
352                }
353                nw0 = nw1;
354            }
355        }
356    }
357
358    private void makeipt(int nw) {
359        int j, l, m, m2, p, q;
360
361        ip[2] = 0;
362        ip[3] = 16;
363        m = 2;
364        for (l = nw; l > 32; l >>= 2) {
365            m2 = m << 1;
366            q = m2 << 3;
367            for (j = m; j < m2; j++) {
368                p = ip[j] << 2;
369                ip[m + j] = p;
370                ip[m2 + j] = p + q;
371            }
372            m = m2;
373        }
374    }
375
376    private void makect(int nc, double[] c, int startc) {
377        int j, nch;
378        double delta, deltaj;
379
380        ip[1] = nc;
381        if (nc > 1) {
382            nch = nc >> 1;
383            delta = 0.785398163397448278999490867136046290 / nch;
384            c[startc] = Math.cos(delta * nch);
385            c[startc + nch] = 0.5 * c[startc];
386            for (j = 1; j < nch; j++) {
387                deltaj = delta * j;
388                c[startc + j] = 0.5 * Math.cos(deltaj);
389                c[startc + nc - j] = 0.5 * Math.sin(deltaj);
390            }
391        }
392    }
393
394    private void cftfsub(int n, double[] a, int offa, int[] ip, int nw, double[] w) {
395        if (n > 8) {
396            if (n > 32) {
397                cftf1st(n, a, offa, w, nw - (n >> 2));
398                if ((ConcurrencyUtils.getNumberOfThreads() > 1) && (n > ConcurrencyUtils.getThreadsBeginN_1D_FFT_2Threads())) {
399                    cftrec4_th(n, a, offa, nw, w);
400                } else if (n > 512) {
401                    cftrec4(n, a, offa, nw, w);
402                } else if (n > 128) {
403                    cftleaf(n, 1, a, offa, nw, w);
404                } else {
405                    cftfx41(n, a, offa, nw, w);
406                }
407                bitrv2(n, ip, a, offa);
408            } else if (n == 32) {
409                cftf161(a, offa, w, nw - 8);
410                bitrv216(a, offa);
411            } else {
412                cftf081(a, offa, w, 0);
413                bitrv208(a, offa);
414            }
415        } else if (n == 8) {
416            cftf040(a, offa);
417        } else if (n == 4) {
418            cftx020(a, offa);
419        }
420    }
421
422    private void cftbsub(int n, double[] a, int offa, int[] ip, int nw, double[] w) {
423        if (n > 8) {
424            if (n > 32) {
425                cftb1st(n, a, offa, w, nw - (n >> 2));
426                if ((ConcurrencyUtils.getNumberOfThreads() > 1) && (n > ConcurrencyUtils.getThreadsBeginN_1D_FFT_2Threads())) {
427                    cftrec4_th(n, a, offa, nw, w);
428                } else if (n > 512) {
429                    cftrec4(n, a, offa, nw, w);
430                } else if (n > 128) {
431                    cftleaf(n, 1, a, offa, nw, w);
432                } else {
433                    cftfx41(n, a, offa, nw, w);
434                }
435                bitrv2conj(n, ip, a, offa);
436            } else if (n == 32) {
437                cftf161(a, offa, w, nw - 8);
438                bitrv216neg(a, offa);
439            } else {
440                cftf081(a, offa, w, 0);
441                bitrv208neg(a, offa);
442            }
443        } else if (n == 8) {
444            cftb040(a, offa);
445        } else if (n == 4) {
446            cftx020(a, offa);
447        }
448    }
449
450    private void bitrv2(int n, int[] ip, double[] a, int offa) {
451        int j1, k1, l, m, nh, nm;
452        double xr, xi, yr, yi;
453        int idx1, idx2;
454
455        m = 1;
456        for (l = n >> 2; l > 8; l >>= 2) {
457            m <<= 1;
458        }
459        nh = n >> 1;
460        nm = 4 * m;
461        if (l == 8) {
462            for (int k = 0; k < m; k++) {
463                for (int j = 0; j < k; j++) {
464                    j1 = 4 * j + 2 * ip[m + k];
465                    k1 = 4 * k + 2 * ip[m + j];
466                    idx1 = offa + j1;
467                    idx2 = offa + k1;
468                    xr = a[idx1];
469                    xi = a[idx1 + 1];
470                    yr = a[idx2];
471                    yi = a[idx2 + 1];
472                    a[idx1] = yr;
473                    a[idx1 + 1] = yi;
474                    a[idx2] = xr;
475                    a[idx2 + 1] = xi;
476                    j1 += nm;
477                    k1 += 2 * nm;
478                    idx1 = offa + j1;
479                    idx2 = offa + k1;
480                    xr = a[idx1];
481                    xi = a[idx1 + 1];
482                    yr = a[idx2];
483                    yi = a[idx2 + 1];
484                    a[idx1] = yr;
485                    a[idx1 + 1] = yi;
486                    a[idx2] = xr;
487                    a[idx2 + 1] = xi;
488                    j1 += nm;
489                    k1 -= nm;
490                    idx1 = offa + j1;
491                    idx2 = offa + k1;
492                    xr = a[idx1];
493                    xi = a[idx1 + 1];
494                    yr = a[idx2];
495                    yi = a[idx2 + 1];
496                    a[idx1] = yr;
497                    a[idx1 + 1] = yi;
498                    a[idx2] = xr;
499                    a[idx2 + 1] = xi;
500                    j1 += nm;
501                    k1 += 2 * nm;
502                    idx1 = offa + j1;
503                    idx2 = offa + k1;
504                    xr = a[idx1];
505                    xi = a[idx1 + 1];
506                    yr = a[idx2];
507                    yi = a[idx2 + 1];
508                    a[idx1] = yr;
509                    a[idx1 + 1] = yi;
510                    a[idx2] = xr;
511                    a[idx2 + 1] = xi;
512                    j1 += nh;
513                    k1 += 2;
514                    idx1 = offa + j1;
515                    idx2 = offa + k1;
516                    xr = a[idx1];
517                    xi = a[idx1 + 1];
518                    yr = a[idx2];
519                    yi = a[idx2 + 1];
520                    a[idx1] = yr;
521                    a[idx1 + 1] = yi;
522                    a[idx2] = xr;
523                    a[idx2 + 1] = xi;
524                    j1 -= nm;
525                    k1 -= 2 * nm;
526                    idx1 = offa + j1;
527                    idx2 = offa + k1;
528                    xr = a[idx1];
529                    xi = a[idx1 + 1];
530                    yr = a[idx2];
531                    yi = a[idx2 + 1];
532                    a[idx1] = yr;
533                    a[idx1 + 1] = yi;
534                    a[idx2] = xr;
535                    a[idx2 + 1] = xi;
536                    j1 -= nm;
537                    k1 += nm;
538                    idx1 = offa + j1;
539                    idx2 = offa + k1;
540                    xr = a[idx1];
541                    xi = a[idx1 + 1];
542                    yr = a[idx2];
543                    yi = a[idx2 + 1];
544                    a[idx1] = yr;
545                    a[idx1 + 1] = yi;
546                    a[idx2] = xr;
547                    a[idx2 + 1] = xi;
548                    j1 -= nm;
549                    k1 -= 2 * nm;
550                    idx1 = offa + j1;
551                    idx2 = offa + k1;
552                    xr = a[idx1];
553                    xi = a[idx1 + 1];
554                    yr = a[idx2];
555                    yi = a[idx2 + 1];
556                    a[idx1] = yr;
557                    a[idx1 + 1] = yi;
558                    a[idx2] = xr;
559                    a[idx2 + 1] = xi;
560                    j1 += 2;
561                    k1 += nh;
562                    idx1 = offa + j1;
563                    idx2 = offa + k1;
564                    xr = a[idx1];
565                    xi = a[idx1 + 1];
566                    yr = a[idx2];
567                    yi = a[idx2 + 1];
568                    a[idx1] = yr;
569                    a[idx1 + 1] = yi;
570                    a[idx2] = xr;
571                    a[idx2 + 1] = xi;
572                    j1 += nm;
573                    k1 += 2 * nm;
574                    idx1 = offa + j1;
575                    idx2 = offa + k1;
576                    xr = a[idx1];
577                    xi = a[idx1 + 1];
578                    yr = a[idx2];
579                    yi = a[idx2 + 1];
580                    a[idx1] = yr;
581                    a[idx1 + 1] = yi;
582                    a[idx2] = xr;
583                    a[idx2 + 1] = xi;
584                    j1 += nm;
585                    k1 -= nm;
586                    idx1 = offa + j1;
587                    idx2 = offa + k1;
588                    xr = a[idx1];
589                    xi = a[idx1 + 1];
590                    yr = a[idx2];
591                    yi = a[idx2 + 1];
592                    a[idx1] = yr;
593                    a[idx1 + 1] = yi;
594                    a[idx2] = xr;
595                    a[idx2 + 1] = xi;
596                    j1 += nm;
597                    k1 += 2 * nm;
598                    idx1 = offa + j1;
599                    idx2 = offa + k1;
600                    xr = a[idx1];
601                    xi = a[idx1 + 1];
602                    yr = a[idx2];
603                    yi = a[idx2 + 1];
604                    a[idx1] = yr;
605                    a[idx1 + 1] = yi;
606                    a[idx2] = xr;
607                    a[idx2 + 1] = xi;
608                    j1 -= nh;
609                    k1 -= 2;
610                    idx1 = offa + j1;
611                    idx2 = offa + k1;
612                    xr = a[idx1];
613                    xi = a[idx1 + 1];
614                    yr = a[idx2];
615                    yi = a[idx2 + 1];
616                    a[idx1] = yr;
617                    a[idx1 + 1] = yi;
618                    a[idx2] = xr;
619                    a[idx2 + 1] = xi;
620                    j1 -= nm;
621                    k1 -= 2 * nm;
622                    idx1 = offa + j1;
623                    idx2 = offa + k1;
624                    xr = a[idx1];
625                    xi = a[idx1 + 1];
626                    yr = a[idx2];
627                    yi = a[idx2 + 1];
628                    a[idx1] = yr;
629                    a[idx1 + 1] = yi;
630                    a[idx2] = xr;
631                    a[idx2 + 1] = xi;
632                    j1 -= nm;
633                    k1 += nm;
634                    idx1 = offa + j1;
635                    idx2 = offa + k1;
636                    xr = a[idx1];
637                    xi = a[idx1 + 1];
638                    yr = a[idx2];
639                    yi = a[idx2 + 1];
640                    a[idx1] = yr;
641                    a[idx1 + 1] = yi;
642                    a[idx2] = xr;
643                    a[idx2 + 1] = xi;
644                    j1 -= nm;
645                    k1 -= 2 * nm;
646                    idx1 = offa + j1;
647                    idx2 = offa + k1;
648                    xr = a[idx1];
649                    xi = a[idx1 + 1];
650                    yr = a[idx2];
651                    yi = a[idx2 + 1];
652                    a[idx1] = yr;
653                    a[idx1 + 1] = yi;
654                    a[idx2] = xr;
655                    a[idx2 + 1] = xi;
656                }
657                k1 = 4 * k + 2 * ip[m + k];
658                j1 = k1 + 2;
659                k1 += nh;
660                idx1 = offa + j1;
661                idx2 = offa + k1;
662                xr = a[idx1];
663                xi = a[idx1 + 1];
664                yr = a[idx2];
665                yi = a[idx2 + 1];
666                a[idx1] = yr;
667                a[idx1 + 1] = yi;
668                a[idx2] = xr;
669                a[idx2 + 1] = xi;
670                j1 += nm;
671                k1 += 2 * nm;
672                idx1 = offa + j1;
673                idx2 = offa + k1;
674                xr = a[idx1];
675                xi = a[idx1 + 1];
676                yr = a[idx2];
677                yi = a[idx2 + 1];
678                a[idx1] = yr;
679                a[idx1 + 1] = yi;
680                a[idx2] = xr;
681                a[idx2 + 1] = xi;
682                j1 += nm;
683                k1 -= nm;
684                idx1 = offa + j1;
685                idx2 = offa + k1;
686                xr = a[idx1];
687                xi = a[idx1 + 1];
688                yr = a[idx2];
689                yi = a[idx2 + 1];
690                a[idx1] = yr;
691                a[idx1 + 1] = yi;
692                a[idx2] = xr;
693                a[idx2 + 1] = xi;
694                j1 -= 2;
695                k1 -= nh;
696                idx1 = offa + j1;
697                idx2 = offa + k1;
698                xr = a[idx1];
699                xi = a[idx1 + 1];
700                yr = a[idx2];
701                yi = a[idx2 + 1];
702                a[idx1] = yr;
703                a[idx1 + 1] = yi;
704                a[idx2] = xr;
705                a[idx2 + 1] = xi;
706                j1 += nh + 2;
707                k1 += nh + 2;
708                idx1 = offa + j1;
709                idx2 = offa + k1;
710                xr = a[idx1];
711                xi = a[idx1 + 1];
712                yr = a[idx2];
713                yi = a[idx2 + 1];
714                a[idx1] = yr;
715                a[idx1 + 1] = yi;
716                a[idx2] = xr;
717                a[idx2 + 1] = xi;
718                j1 -= nh - nm;
719                k1 += 2 * nm - 2;
720                idx1 = offa + j1;
721                idx2 = offa + k1;
722                xr = a[idx1];
723                xi = a[idx1 + 1];
724                yr = a[idx2];
725                yi = a[idx2 + 1];
726                a[idx1] = yr;
727                a[idx1 + 1] = yi;
728                a[idx2] = xr;
729                a[idx2 + 1] = xi;
730            }
731        } else {
732            for (int k = 0; k < m; k++) {
733                for (int j = 0; j < k; j++) {
734                    j1 = 4 * j + ip[m + k];
735                    k1 = 4 * k + ip[m + j];
736                    idx1 = offa + j1;
737                    idx2 = offa + k1;
738                    xr = a[idx1];
739                    xi = a[idx1 + 1];
740                    yr = a[idx2];
741                    yi = a[idx2 + 1];
742                    a[idx1] = yr;
743                    a[idx1 + 1] = yi;
744                    a[idx2] = xr;
745                    a[idx2 + 1] = xi;
746                    j1 += nm;
747                    k1 += nm;
748                    idx1 = offa + j1;
749                    idx2 = offa + k1;
750                    xr = a[idx1];
751                    xi = a[idx1 + 1];
752                    yr = a[idx2];
753                    yi = a[idx2 + 1];
754                    a[idx1] = yr;
755                    a[idx1 + 1] = yi;
756                    a[idx2] = xr;
757                    a[idx2 + 1] = xi;
758                    j1 += nh;
759                    k1 += 2;
760                    idx1 = offa + j1;
761                    idx2 = offa + k1;
762                    xr = a[idx1];
763                    xi = a[idx1 + 1];
764                    yr = a[idx2];
765                    yi = a[idx2 + 1];
766                    a[idx1] = yr;
767                    a[idx1 + 1] = yi;
768                    a[idx2] = xr;
769                    a[idx2 + 1] = xi;
770                    j1 -= nm;
771                    k1 -= nm;
772                    idx1 = offa + j1;
773                    idx2 = offa + k1;
774                    xr = a[idx1];
775                    xi = a[idx1 + 1];
776                    yr = a[idx2];
777                    yi = a[idx2 + 1];
778                    a[idx1] = yr;
779                    a[idx1 + 1] = yi;
780                    a[idx2] = xr;
781                    a[idx2 + 1] = xi;
782                    j1 += 2;
783                    k1 += nh;
784                    idx1 = offa + j1;
785                    idx2 = offa + k1;
786                    xr = a[idx1];
787                    xi = a[idx1 + 1];
788                    yr = a[idx2];
789                    yi = a[idx2 + 1];
790                    a[idx1] = yr;
791                    a[idx1 + 1] = yi;
792                    a[idx2] = xr;
793                    a[idx2 + 1] = xi;
794                    j1 += nm;
795                    k1 += nm;
796                    idx1 = offa + j1;
797                    idx2 = offa + k1;
798                    xr = a[idx1];
799                    xi = a[idx1 + 1];
800                    yr = a[idx2];
801                    yi = a[idx2 + 1];
802                    a[idx1] = yr;
803                    a[idx1 + 1] = yi;
804                    a[idx2] = xr;
805                    a[idx2 + 1] = xi;
806                    j1 -= nh;
807                    k1 -= 2;
808                    idx1 = offa + j1;
809                    idx2 = offa + k1;
810                    xr = a[idx1];
811                    xi = a[idx1 + 1];
812                    yr = a[idx2];
813                    yi = a[idx2 + 1];
814                    a[idx1] = yr;
815                    a[idx1 + 1] = yi;
816                    a[idx2] = xr;
817                    a[idx2 + 1] = xi;
818                    j1 -= nm;
819                    k1 -= nm;
820                    idx1 = offa + j1;
821                    idx2 = offa + k1;
822                    xr = a[idx1];
823                    xi = a[idx1 + 1];
824                    yr = a[idx2];
825                    yi = a[idx2 + 1];
826                    a[idx1] = yr;
827                    a[idx1 + 1] = yi;
828                    a[idx2] = xr;
829                    a[idx2 + 1] = xi;
830                }
831                k1 = 4 * k + ip[m + k];
832                j1 = k1 + 2;
833                k1 += nh;
834                idx1 = offa + j1;
835                idx2 = offa + k1;
836                xr = a[idx1];
837                xi = a[idx1 + 1];
838                yr = a[idx2];
839                yi = a[idx2 + 1];
840                a[idx1] = yr;
841                a[idx1 + 1] = yi;
842                a[idx2] = xr;
843                a[idx2 + 1] = xi;
844                j1 += nm;
845                k1 += nm;
846                idx1 = offa + j1;
847                idx2 = offa + k1;
848                xr = a[idx1];
849                xi = a[idx1 + 1];
850                yr = a[idx2];
851                yi = a[idx2 + 1];
852                a[idx1] = yr;
853                a[idx1 + 1] = yi;
854                a[idx2] = xr;
855                a[idx2 + 1] = xi;
856            }
857        }
858    }
859
860    private void bitrv2conj(int n, int[] ip, double[] a, int offa) {
861        int j1, k1, l, m, nh, nm;
862        double xr, xi, yr, yi;
863        int idx1, idx2;
864
865        m = 1;
866        for (l = n >> 2; l > 8; l >>= 2) {
867            m <<= 1;
868        }
869        nh = n >> 1;
870        nm = 4 * m;
871        if (l == 8) {
872            for (int k = 0; k < m; k++) {
873                for (int j = 0; j < k; j++) {
874                    j1 = 4 * j + 2 * ip[m + k];
875                    k1 = 4 * k + 2 * ip[m + j];
876                    idx1 = offa + j1;
877                    idx2 = offa + k1;
878                    xr = a[idx1];
879                    xi = -a[idx1 + 1];
880                    yr = a[idx2];
881                    yi = -a[idx2 + 1];
882                    a[idx1] = yr;
883                    a[idx1 + 1] = yi;
884                    a[idx2] = xr;
885                    a[idx2 + 1] = xi;
886                    j1 += nm;
887                    k1 += 2 * nm;
888                    idx1 = offa + j1;
889                    idx2 = offa + k1;
890                    xr = a[idx1];
891                    xi = -a[idx1 + 1];
892                    yr = a[idx2];
893                    yi = -a[idx2 + 1];
894                    a[idx1] = yr;
895                    a[idx1 + 1] = yi;
896                    a[idx2] = xr;
897                    a[idx2 + 1] = xi;
898                    j1 += nm;
899                    k1 -= nm;
900                    idx1 = offa + j1;
901                    idx2 = offa + k1;
902                    xr = a[idx1];
903                    xi = -a[idx1 + 1];
904                    yr = a[idx2];
905                    yi = -a[idx2 + 1];
906                    a[idx1] = yr;
907                    a[idx1 + 1] = yi;
908                    a[idx2] = xr;
909                    a[idx2 + 1] = xi;
910                    j1 += nm;
911                    k1 += 2 * nm;
912                    idx1 = offa + j1;
913                    idx2 = offa + k1;
914                    xr = a[idx1];
915                    xi = -a[idx1 + 1];
916                    yr = a[idx2];
917                    yi = -a[idx2 + 1];
918                    a[idx1] = yr;
919                    a[idx1 + 1] = yi;
920                    a[idx2] = xr;
921                    a[idx2 + 1] = xi;
922                    j1 += nh;
923                    k1 += 2;
924                    idx1 = offa + j1;
925                    idx2 = offa + k1;
926                    xr = a[idx1];
927                    xi = -a[idx1 + 1];
928                    yr = a[idx2];
929                    yi = -a[idx2 + 1];
930                    a[idx1] = yr;
931                    a[idx1 + 1] = yi;
932                    a[idx2] = xr;
933                    a[idx2 + 1] = xi;
934                    j1 -= nm;
935                    k1 -= 2 * nm;
936                    idx1 = offa + j1;
937                    idx2 = offa + k1;
938                    xr = a[idx1];
939                    xi = -a[idx1 + 1];
940                    yr = a[idx2];
941                    yi = -a[idx2 + 1];
942                    a[idx1] = yr;
943                    a[idx1 + 1] = yi;
944                    a[idx2] = xr;
945                    a[idx2 + 1] = xi;
946                    j1 -= nm;
947                    k1 += nm;
948                    idx1 = offa + j1;
949                    idx2 = offa + k1;
950                    xr = a[idx1];
951                    xi = -a[idx1 + 1];
952                    yr = a[idx2];
953                    yi = -a[idx2 + 1];
954                    a[idx1] = yr;
955                    a[idx1 + 1] = yi;
956                    a[idx2] = xr;
957                    a[idx2 + 1] = xi;
958                    j1 -= nm;
959                    k1 -= 2 * nm;
960                    idx1 = offa + j1;
961                    idx2 = offa + k1;
962                    xr = a[idx1];
963                    xi = -a[idx1 + 1];
964                    yr = a[idx2];
965                    yi = -a[idx2 + 1];
966                    a[idx1] = yr;
967                    a[idx1 + 1] = yi;
968                    a[idx2] = xr;
969                    a[idx2 + 1] = xi;
970                    j1 += 2;
971                    k1 += nh;
972                    idx1 = offa + j1;
973                    idx2 = offa + k1;
974                    xr = a[idx1];
975                    xi = -a[idx1 + 1];
976                    yr = a[idx2];
977                    yi = -a[idx2 + 1];
978                    a[idx1] = yr;
979                    a[idx1 + 1] = yi;
980                    a[idx2] = xr;
981                    a[idx2 + 1] = xi;
982                    j1 += nm;
983                    k1 += 2 * nm;
984                    idx1 = offa + j1;
985                    idx2 = offa + k1;
986                    xr = a[idx1];
987                    xi = -a[idx1 + 1];
988                    yr = a[idx2];
989                    yi = -a[idx2 + 1];
990                    a[idx1] = yr;
991                    a[idx1 + 1] = yi;
992                    a[idx2] = xr;
993                    a[idx2 + 1] = xi;
994                    j1 += nm;
995                    k1 -= nm;
996                    idx1 = offa + j1;
997                    idx2 = offa + k1;
998                    xr = a[idx1];
999                    xi = -a[idx1 + 1];
1000                    yr = a[idx2];
1001                    yi = -a[idx2 + 1];
1002                    a[idx1] = yr;
1003                    a[idx1 + 1] = yi;
1004                    a[idx2] = xr;
1005                    a[idx2 + 1] = xi;
1006                    j1 += nm;
1007                    k1 += 2 * nm;
1008                    idx1 = offa + j1;
1009                    idx2 = offa + k1;
1010                    xr = a[idx1];
1011                    xi = -a[idx1 + 1];
1012                    yr = a[idx2];
1013                    yi = -a[idx2 + 1];
1014                    a[idx1] = yr;
1015                    a[idx1 + 1] = yi;
1016                    a[idx2] = xr;
1017                    a[idx2 + 1] = xi;
1018                    j1 -= nh;
1019                    k1 -= 2;
1020                    idx1 = offa + j1;
1021                    idx2 = offa + k1;
1022                    xr = a[idx1];
1023                    xi = -a[idx1 + 1];
1024                    yr = a[idx2];
1025                    yi = -a[idx2 + 1];
1026                    a[idx1] = yr;
1027                    a[idx1 + 1] = yi;
1028                    a[idx2] = xr;
1029                    a[idx2 + 1] = xi;
1030                    j1 -= nm;
1031                    k1 -= 2 * nm;
1032                    idx1 = offa + j1;
1033                    idx2 = offa + k1;
1034                    xr = a[idx1];
1035                    xi = -a[idx1 + 1];
1036                    yr = a[idx2];
1037                    yi = -a[idx2 + 1];
1038                    a[idx1] = yr;
1039                    a[idx1 + 1] = yi;
1040                    a[idx2] = xr;
1041                    a[idx2 + 1] = xi;
1042                    j1 -= nm;
1043                    k1 += nm;
1044                    idx1 = offa + j1;
1045                    idx2 = offa + k1;
1046                    xr = a[idx1];
1047                    xi = -a[idx1 + 1];
1048                    yr = a[idx2];
1049                    yi = -a[idx2 + 1];
1050                    a[idx1] = yr;
1051                    a[idx1 + 1] = yi;
1052                    a[idx2] = xr;
1053                    a[idx2 + 1] = xi;
1054                    j1 -= nm;
1055                    k1 -= 2 * nm;
1056                    idx1 = offa + j1;
1057                    idx2 = offa + k1;
1058                    xr = a[idx1];
1059                    xi = -a[idx1 + 1];
1060                    yr = a[idx2];
1061                    yi = -a[idx2 + 1];
1062                    a[idx1] = yr;
1063                    a[idx1 + 1] = yi;
1064                    a[idx2] = xr;
1065                    a[idx2 + 1] = xi;
1066                }
1067                k1 = 4 * k + 2 * ip[m + k];
1068                j1 = k1 + 2;
1069                k1 += nh;
1070                idx1 = offa + j1;
1071                idx2 = offa + k1;
1072                a[idx1 - 1] = -a[idx1 - 1];
1073                xr = a[idx1];
1074                xi = -a[idx1 + 1];
1075                yr = a[idx2];
1076                yi = -a[idx2 + 1];
1077                a[idx1] = yr;
1078                a[idx1 + 1] = yi;
1079                a[idx2] = xr;
1080                a[idx2 + 1] = xi;
1081                a[idx2 + 3] = -a[idx2 + 3];
1082                j1 += nm;
1083                k1 += 2 * nm;
1084                idx1 = offa + j1;
1085                idx2 = offa + k1;
1086                xr = a[idx1];
1087                xi = -a[idx1 + 1];
1088                yr = a[idx2];
1089                yi = -a[idx2 + 1];
1090                a[idx1] = yr;
1091                a[idx1 + 1] = yi;
1092                a[idx2] = xr;
1093                a[idx2 + 1] = xi;
1094                j1 += nm;
1095                k1 -= nm;
1096                idx1 = offa + j1;
1097                idx2 = offa + k1;
1098                xr = a[idx1];
1099                xi = -a[idx1 + 1];
1100                yr = a[idx2];
1101                yi = -a[idx2 + 1];
1102                a[idx1] = yr;
1103                a[idx1 + 1] = yi;
1104                a[idx2] = xr;
1105                a[idx2 + 1] = xi;
1106                j1 -= 2;
1107                k1 -= nh;
1108                idx1 = offa + j1;
1109                idx2 = offa + k1;
1110                xr = a[idx1];
1111                xi = -a[idx1 + 1];
1112                yr = a[idx2];
1113                yi = -a[idx2 + 1];
1114                a[idx1] = yr;
1115                a[idx1 + 1] = yi;
1116                a[idx2] = xr;
1117                a[idx2 + 1] = xi;
1118                j1 += nh + 2;
1119                k1 += nh + 2;
1120                idx1 = offa + j1;
1121                idx2 = offa + k1;
1122                xr = a[idx1];
1123                xi = -a[idx1 + 1];
1124                yr = a[idx2];
1125                yi = -a[idx2 + 1];
1126                a[idx1] = yr;
1127                a[idx1 + 1] = yi;
1128                a[idx2] = xr;
1129                a[idx2 + 1] = xi;
1130                j1 -= nh - nm;
1131                k1 += 2 * nm - 2;
1132                idx1 = offa + j1;
1133                idx2 = offa + k1;
1134                a[idx1 - 1] = -a[idx1 - 1];
1135                xr = a[idx1];
1136                xi = -a[idx1 + 1];
1137                yr = a[idx2];
1138                yi = -a[idx2 + 1];
1139                a[idx1] = yr;
1140                a[idx1 + 1] = yi;
1141                a[idx2] = xr;
1142                a[idx2 + 1] = xi;
1143                a[idx2 + 3] = -a[idx2 + 3];
1144            }
1145        } else {
1146            for (int k = 0; k < m; k++) {
1147                for (int j = 0; j < k; j++) {
1148                    j1 = 4 * j + ip[m + k];
1149                    k1 = 4 * k + ip[m + j];
1150                    idx1 = offa + j1;
1151                    idx2 = offa + k1;
1152                    xr = a[idx1];
1153                    xi = -a[idx1 + 1];
1154                    yr = a[idx2];
1155                    yi = -a[idx2 + 1];
1156                    a[idx1] = yr;
1157                    a[idx1 + 1] = yi;
1158                    a[idx2] = xr;
1159                    a[idx2 + 1] = xi;
1160                    j1 += nm;
1161                    k1 += nm;
1162                    idx1 = offa + j1;
1163                    idx2 = offa + k1;
1164                    xr = a[idx1];
1165                    xi = -a[idx1 + 1];
1166                    yr = a[idx2];
1167                    yi = -a[idx2 + 1];
1168                    a[idx1] = yr;
1169                    a[idx1 + 1] = yi;
1170                    a[idx2] = xr;
1171                    a[idx2 + 1] = xi;
1172                    j1 += nh;
1173                    k1 += 2;
1174                    idx1 = offa + j1;
1175                    idx2 = offa + k1;
1176                    xr = a[idx1];
1177                    xi = -a[idx1 + 1];
1178                    yr = a[idx2];
1179                    yi = -a[idx2 + 1];
1180                    a[idx1] = yr;
1181                    a[idx1 + 1] = yi;
1182                    a[idx2] = xr;
1183                    a[idx2 + 1] = xi;
1184                    j1 -= nm;
1185                    k1 -= nm;
1186                    idx1 = offa + j1;
1187                    idx2 = offa + k1;
1188                    xr = a[idx1];
1189                    xi = -a[idx1 + 1];
1190                    yr = a[idx2];
1191                    yi = -a[idx2 + 1];
1192                    a[idx1] = yr;
1193                    a[idx1 + 1] = yi;
1194                    a[idx2] = xr;
1195                    a[idx2 + 1] = xi;
1196                    j1 += 2;
1197                    k1 += nh;
1198                    idx1 = offa + j1;
1199                    idx2 = offa + k1;
1200                    xr = a[idx1];
1201                    xi = -a[idx1 + 1];
1202                    yr = a[idx2];
1203                    yi = -a[idx2 + 1];
1204                    a[idx1] = yr;
1205                    a[idx1 + 1] = yi;
1206                    a[idx2] = xr;
1207                    a[idx2 + 1] = xi;
1208                    j1 += nm;
1209                    k1 += nm;
1210                    idx1 = offa + j1;
1211                    idx2 = offa + k1;
1212                    xr = a[idx1];
1213                    xi = -a[idx1 + 1];
1214                    yr = a[idx2];
1215                    yi = -a[idx2 + 1];
1216                    a[idx1] = yr;
1217                    a[idx1 + 1] = yi;
1218                    a[idx2] = xr;
1219                    a[idx2 + 1] = xi;
1220                    j1 -= nh;
1221                    k1 -= 2;
1222                    idx1 = offa + j1;
1223                    idx2 = offa + k1;
1224                    xr = a[idx1];
1225                    xi = -a[idx1 + 1];
1226                    yr = a[idx2];
1227                    yi = -a[idx2 + 1];
1228                    a[idx1] = yr;
1229                    a[idx1 + 1] = yi;
1230                    a[idx2] = xr;
1231                    a[idx2 + 1] = xi;
1232                    j1 -= nm;
1233                    k1 -= nm;
1234                    idx1 = offa + j1;
1235                    idx2 = offa + k1;
1236                    xr = a[idx1];
1237                    xi = -a[idx1 + 1];
1238                    yr = a[idx2];
1239                    yi = -a[idx2 + 1];
1240                    a[idx1] = yr;
1241                    a[idx1 + 1] = yi;
1242                    a[idx2] = xr;
1243                    a[idx2 + 1] = xi;
1244                }
1245                k1 = 4 * k + ip[m + k];
1246                j1 = k1 + 2;
1247                k1 += nh;
1248                idx1 = offa + j1;
1249                idx2 = offa + k1;
1250                a[idx1 - 1] = -a[idx1 - 1];
1251                xr = a[idx1];
1252                xi = -a[idx1 + 1];
1253                yr = a[idx2];
1254                yi = -a[idx2 + 1];
1255                a[idx1] = yr;
1256                a[idx1 + 1] = yi;
1257                a[idx2] = xr;
1258                a[idx2 + 1] = xi;
1259                a[idx2 + 3] = -a[idx2 + 3];
1260                j1 += nm;
1261                k1 += nm;
1262                idx1 = offa + j1;
1263                idx2 = offa + k1;
1264                a[idx1 - 1] = -a[idx1 - 1];
1265                xr = a[idx1];
1266                xi = -a[idx1 + 1];
1267                yr = a[idx2];
1268                yi = -a[idx2 + 1];
1269                a[idx1] = yr;
1270                a[idx1 + 1] = yi;
1271                a[idx2] = xr;
1272                a[idx2 + 1] = xi;
1273                a[idx2 + 3] = -a[idx2 + 3];
1274            }
1275        }
1276    }
1277
1278    private void bitrv216(double[] a, int offa) {
1279        double x1r, x1i, x2r, x2i, x3r, x3i, x4r, x4i, x5r, x5i, x7r, x7i, x8r, x8i, x10r, x10i, x11r, x11i, x12r, x12i, x13r, x13i, x14r, x14i;
1280
1281        x1r = a[offa + 2];
1282        x1i = a[offa + 3];
1283        x2r = a[offa + 4];
1284        x2i = a[offa + 5];
1285        x3r = a[offa + 6];
1286        x3i = a[offa + 7];
1287        x4r = a[offa + 8];
1288        x4i = a[offa + 9];
1289        x5r = a[offa + 10];
1290        x5i = a[offa + 11];
1291        x7r = a[offa + 14];
1292        x7i = a[offa + 15];
1293        x8r = a[offa + 16];
1294        x8i = a[offa + 17];
1295        x10r = a[offa + 20];
1296        x10i = a[offa + 21];
1297        x11r = a[offa + 22];
1298        x11i = a[offa + 23];
1299        x12r = a[offa + 24];
1300        x12i = a[offa + 25];
1301        x13r = a[offa + 26];
1302        x13i = a[offa + 27];
1303        x14r = a[offa + 28];
1304        x14i = a[offa + 29];
1305        a[offa + 2] = x8r;
1306        a[offa + 3] = x8i;
1307        a[offa + 4] = x4r;
1308        a[offa + 5] = x4i;
1309        a[offa + 6] = x12r;
1310        a[offa + 7] = x12i;
1311        a[offa + 8] = x2r;
1312        a[offa + 9] = x2i;
1313        a[offa + 10] = x10r;
1314        a[offa + 11] = x10i;
1315        a[offa + 14] = x14r;
1316        a[offa + 15] = x14i;
1317        a[offa + 16] = x1r;
1318        a[offa + 17] = x1i;
1319        a[offa + 20] = x5r;
1320        a[offa + 21] = x5i;
1321        a[offa + 22] = x13r;
1322        a[offa + 23] = x13i;
1323        a[offa + 24] = x3r;
1324        a[offa + 25] = x3i;
1325        a[offa + 26] = x11r;
1326        a[offa + 27] = x11i;
1327        a[offa + 28] = x7r;
1328        a[offa + 29] = x7i;
1329    }
1330
1331    private void bitrv216neg(double[] a, int offa) {
1332        double x1r, x1i, x2r, x2i, x3r, x3i, x4r, x4i, x5r, x5i, x6r, x6i, x7r, x7i, x8r, x8i, x9r, x9i, x10r, x10i, x11r, x11i, x12r, x12i, x13r, x13i, x14r, x14i, x15r, x15i;
1333
1334        x1r = a[offa + 2];
1335        x1i = a[offa + 3];
1336        x2r = a[offa + 4];
1337        x2i = a[offa + 5];
1338        x3r = a[offa + 6];
1339        x3i = a[offa + 7];
1340        x4r = a[offa + 8];
1341        x4i = a[offa + 9];
1342        x5r = a[offa + 10];
1343        x5i = a[offa + 11];
1344        x6r = a[offa + 12];
1345        x6i = a[offa + 13];
1346        x7r = a[offa + 14];
1347        x7i = a[offa + 15];
1348        x8r = a[offa + 16];
1349        x8i = a[offa + 17];
1350        x9r = a[offa + 18];
1351        x9i = a[offa + 19];
1352        x10r = a[offa + 20];
1353        x10i = a[offa + 21];
1354        x11r = a[offa + 22];
1355        x11i = a[offa + 23];
1356        x12r = a[offa + 24];
1357        x12i = a[offa + 25];
1358        x13r = a[offa + 26];
1359        x13i = a[offa + 27];
1360        x14r = a[offa + 28];
1361        x14i = a[offa + 29];
1362        x15r = a[offa + 30];
1363        x15i = a[offa + 31];
1364        a[offa + 2] = x15r;
1365        a[offa + 3] = x15i;
1366        a[offa + 4] = x7r;
1367        a[offa + 5] = x7i;
1368        a[offa + 6] = x11r;
1369        a[offa + 7] = x11i;
1370        a[offa + 8] = x3r;
1371        a[offa + 9] = x3i;
1372        a[offa + 10] = x13r;
1373        a[offa + 11] = x13i;
1374        a[offa + 12] = x5r;
1375        a[offa + 13] = x5i;
1376        a[offa + 14] = x9r;
1377        a[offa + 15] = x9i;
1378        a[offa + 16] = x1r;
1379        a[offa + 17] = x1i;
1380        a[offa + 18] = x14r;
1381        a[offa + 19] = x14i;
1382        a[offa + 20] = x6r;
1383        a[offa + 21] = x6i;
1384        a[offa + 22] = x10r;
1385        a[offa + 23] = x10i;
1386        a[offa + 24] = x2r;
1387        a[offa + 25] = x2i;
1388        a[offa + 26] = x12r;
1389        a[offa + 27] = x12i;
1390        a[offa + 28] = x4r;
1391        a[offa + 29] = x4i;
1392        a[offa + 30] = x8r;
1393        a[offa + 31] = x8i;
1394    }
1395
1396    private void bitrv208(double[] a, int offa) {
1397        double x1r, x1i, x3r, x3i, x4r, x4i, x6r, x6i;
1398
1399        x1r = a[offa + 2];
1400        x1i = a[offa + 3];
1401        x3r = a[offa + 6];
1402        x3i = a[offa + 7];
1403        x4r = a[offa + 8];
1404        x4i = a[offa + 9];
1405        x6r = a[offa + 12];
1406        x6i = a[offa + 13];
1407        a[offa + 2] = x4r;
1408        a[offa + 3] = x4i;
1409        a[offa + 6] = x6r;
1410        a[offa + 7] = x6i;
1411        a[offa + 8] = x1r;
1412        a[offa + 9] = x1i;
1413        a[offa + 12] = x3r;
1414        a[offa + 13] = x3i;
1415    }
1416
1417    private void bitrv208neg(double[] a, int offa) {
1418        double x1r, x1i, x2r, x2i, x3r, x3i, x4r, x4i, x5r, x5i, x6r, x6i, x7r, x7i;
1419
1420        x1r = a[offa + 2];
1421        x1i = a[offa + 3];
1422        x2r = a[offa + 4];
1423        x2i = a[offa + 5];
1424        x3r = a[offa + 6];
1425        x3i = a[offa + 7];
1426        x4r = a[offa + 8];
1427        x4i = a[offa + 9];
1428        x5r = a[offa + 10];
1429        x5i = a[offa + 11];
1430        x6r = a[offa + 12];
1431        x6i = a[offa + 13];
1432        x7r = a[offa + 14];
1433        x7i = a[offa + 15];
1434        a[offa + 2] = x7r;
1435        a[offa + 3] = x7i;
1436        a[offa + 4] = x3r;
1437        a[offa + 5] = x3i;
1438        a[offa + 6] = x5r;
1439        a[offa + 7] = x5i;
1440        a[offa + 8] = x1r;
1441        a[offa + 9] = x1i;
1442        a[offa + 10] = x6r;
1443        a[offa + 11] = x6i;
1444        a[offa + 12] = x2r;
1445        a[offa + 13] = x2i;
1446        a[offa + 14] = x4r;
1447        a[offa + 15] = x4i;
1448    }
1449
1450    private void cftf1st(int n, double[] a, int offa, double[] w, int startw) {
1451        int j0, j1, j2, j3, k, m, mh;
1452        double wn4r, csc1, csc3, wk1r, wk1i, wk3r, wk3i, wd1r, wd1i, wd3r, wd3i;
1453        double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i, y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i;
1454        int idx0, idx1, idx2, idx3, idx4, idx5;
1455        mh = n >> 3;
1456        m = 2 * mh;
1457        j1 = m;
1458        j2 = j1 + m;
1459        j3 = j2 + m;
1460        idx1 = offa + j1;
1461        idx2 = offa + j2;
1462        idx3 = offa + j3;
1463        x0r = a[offa] + a[idx2];
1464        x0i = a[offa + 1] + a[idx2 + 1];
1465        x1r = a[offa] - a[idx2];
1466        x1i = a[offa + 1] - a[idx2 + 1];
1467        x2r = a[idx1] + a[idx3];
1468        x2i = a[idx1 + 1] + a[idx3 + 1];
1469        x3r = a[idx1] - a[idx3];
1470        x3i = a[idx1 + 1] - a[idx3 + 1];
1471        a[offa] = x0r + x2r;
1472        a[offa + 1] = x0i + x2i;
1473        a[idx1] = x0r - x2r;
1474        a[idx1 + 1] = x0i - x2i;
1475        a[idx2] = x1r - x3i;
1476        a[idx2 + 1] = x1i + x3r;
1477        a[idx3] = x1r + x3i;
1478        a[idx3 + 1] = x1i - x3r;
1479        wn4r = w[startw + 1];
1480        csc1 = w[startw + 2];
1481        csc3 = w[startw + 3];
1482        wd1r = 1;
1483        wd1i = 0;
1484        wd3r = 1;
1485        wd3i = 0;
1486        k = 0;
1487        for (int j = 2; j < mh - 2; j += 4) {
1488            k += 4;
1489            idx4 = startw + k;
1490            wk1r = csc1 * (wd1r + w[idx4]);
1491            wk1i = csc1 * (wd1i + w[idx4 + 1]);
1492            wk3r = csc3 * (wd3r + w[idx4 + 2]);
1493            wk3i = csc3 * (wd3i + w[idx4 + 3]);
1494            wd1r = w[idx4];
1495            wd1i = w[idx4 + 1];
1496            wd3r = w[idx4 + 2];
1497            wd3i = w[idx4 + 3];
1498            j1 = j + m;
1499            j2 = j1 + m;
1500            j3 = j2 + m;
1501            idx1 = offa + j1;
1502            idx2 = offa + j2;
1503            idx3 = offa + j3;
1504            idx5 = offa + j;
1505            x0r = a[idx5] + a[idx2];
1506            x0i = a[idx5 + 1] + a[idx2 + 1];
1507            x1r = a[idx5] - a[idx2];
1508            x1i = a[idx5 + 1] - a[idx2 + 1];
1509            y0r = a[idx5 + 2] + a[idx2 + 2];
1510            y0i = a[idx5 + 3] + a[idx2 + 3];
1511            y1r = a[idx5 + 2] - a[idx2 + 2];
1512            y1i = a[idx5 + 3] - a[idx2 + 3];
1513            x2r = a[idx1] + a[idx3];
1514            x2i = a[idx1 + 1] + a[idx3 + 1];
1515            x3r = a[idx1] - a[idx3];
1516            x3i = a[idx1 + 1] - a[idx3 + 1];
1517            y2r = a[idx1 + 2] + a[idx3 + 2];
1518            y2i = a[idx1 + 3] + a[idx3 + 3];
1519            y3r = a[idx1 + 2] - a[idx3 + 2];
1520            y3i = a[idx1 + 3] - a[idx3 + 3];
1521            a[idx5] = x0r + x2r;
1522            a[idx5 + 1] = x0i + x2i;
1523            a[idx5 + 2] = y0r + y2r;
1524            a[idx5 + 3] = y0i + y2i;
1525            a[idx1] = x0r - x2r;
1526            a[idx1 + 1] = x0i - x2i;
1527            a[idx1 + 2] = y0r - y2r;
1528            a[idx1 + 3] = y0i - y2i;
1529            x0r = x1r - x3i;
1530            x0i = x1i + x3r;
1531            a[idx2] = wk1r * x0r - wk1i * x0i;
1532            a[idx2 + 1] = wk1r * x0i + wk1i * x0r;
1533            x0r = y1r - y3i;
1534            x0i = y1i + y3r;
1535            a[idx2 + 2] = wd1r * x0r - wd1i * x0i;
1536            a[idx2 + 3] = wd1r * x0i + wd1i * x0r;
1537            x0r = x1r + x3i;
1538            x0i = x1i - x3r;
1539            a[idx3] = wk3r * x0r + wk3i * x0i;
1540            a[idx3 + 1] = wk3r * x0i - wk3i * x0r;
1541            x0r = y1r + y3i;
1542            x0i = y1i - y3r;
1543            a[idx3 + 2] = wd3r * x0r + wd3i * x0i;
1544            a[idx3 + 3] = wd3r * x0i - wd3i * x0r;
1545            j0 = m - j;
1546            j1 = j0 + m;
1547            j2 = j1 + m;
1548            j3 = j2 + m;
1549            idx0 = offa + j0;
1550            idx1 = offa + j1;
1551            idx2 = offa + j2;
1552            idx3 = offa + j3;
1553            x0r = a[idx0] + a[idx2];
1554            x0i = a[idx0 + 1] + a[idx2 + 1];
1555            x1r = a[idx0] - a[idx2];
1556            x1i = a[idx0 + 1] - a[idx2 + 1];
1557            y0r = a[idx0 - 2] + a[idx2 - 2];
1558            y0i = a[idx0 - 1] + a[idx2 - 1];
1559            y1r = a[idx0 - 2] - a[idx2 - 2];
1560            y1i = a[idx0 - 1] - a[idx2 - 1];
1561            x2r = a[idx1] + a[idx3];
1562            x2i = a[idx1 + 1] + a[idx3 + 1];
1563            x3r = a[idx1] - a[idx3];
1564            x3i = a[idx1 + 1] - a[idx3 + 1];
1565            y2r = a[idx1 - 2] + a[idx3 - 2];
1566            y2i = a[idx1 - 1] + a[idx3 - 1];
1567            y3r = a[idx1 - 2] - a[idx3 - 2];
1568            y3i = a[idx1 - 1] - a[idx3 - 1];
1569            a[idx0] = x0r + x2r;
1570            a[idx0 + 1] = x0i + x2i;
1571            a[idx0 - 2] = y0r + y2r;
1572            a[idx0 - 1] = y0i + y2i;
1573            a[idx1] = x0r - x2r;
1574            a[idx1 + 1] = x0i - x2i;
1575            a[idx1 - 2] = y0r - y2r;
1576            a[idx1 - 1] = y0i - y2i;
1577            x0r = x1r - x3i;
1578            x0i = x1i + x3r;
1579            a[idx2] = wk1i * x0r - wk1r * x0i;
1580            a[idx2 + 1] = wk1i * x0i + wk1r * x0r;
1581            x0r = y1r - y3i;
1582            x0i = y1i + y3r;
1583            a[idx2 - 2] = wd1i * x0r - wd1r * x0i;
1584            a[idx2 - 1] = wd1i * x0i + wd1r * x0r;
1585            x0r = x1r + x3i;
1586            x0i = x1i - x3r;
1587            a[idx3] = wk3i * x0r + wk3r * x0i;
1588            a[idx3 + 1] = wk3i * x0i - wk3r * x0r;
1589            x0r = y1r + y3i;
1590            x0i = y1i - y3r;
1591            a[offa + j3 - 2] = wd3i * x0r + wd3r * x0i;
1592            a[offa + j3 - 1] = wd3i * x0i - wd3r * x0r;
1593        }
1594        wk1r = csc1 * (wd1r + wn4r);
1595        wk1i = csc1 * (wd1i + wn4r);
1596        wk3r = csc3 * (wd3r - wn4r);
1597        wk3i = csc3 * (wd3i - wn4r);
1598        j0 = mh;
1599        j1 = j0 + m;
1600        j2 = j1 + m;
1601        j3 = j2 + m;
1602        idx0 = offa + j0;
1603        idx1 = offa + j1;
1604        idx2 = offa + j2;
1605        idx3 = offa + j3;
1606        x0r = a[idx0 - 2] + a[idx2 - 2];
1607        x0i = a[idx0 - 1] + a[idx2 - 1];
1608        x1r = a[idx0 - 2] - a[idx2 - 2];
1609        x1i = a[idx0 - 1] - a[idx2 - 1];
1610        x2r = a[idx1 - 2] + a[idx3 - 2];
1611        x2i = a[idx1 - 1] + a[idx3 - 1];
1612        x3r = a[idx1 - 2] - a[idx3 - 2];
1613        x3i = a[idx1 - 1] - a[idx3 - 1];
1614        a[idx0 - 2] = x0r + x2r;
1615        a[idx0 - 1] = x0i + x2i;
1616        a[idx1 - 2] = x0r - x2r;
1617        a[idx1 - 1] = x0i - x2i;
1618        x0r = x1r - x3i;
1619        x0i = x1i + x3r;
1620        a[idx2 - 2] = wk1r * x0r - wk1i * x0i;
1621        a[idx2 - 1] = wk1r * x0i + wk1i * x0r;
1622        x0r = x1r + x3i;
1623        x0i = x1i - x3r;
1624        a[idx3 - 2] = wk3r * x0r + wk3i * x0i;
1625        a[idx3 - 1] = wk3r * x0i - wk3i * x0r;
1626        x0r = a[idx0] + a[idx2];
1627        x0i = a[idx0 + 1] + a[idx2 + 1];
1628        x1r = a[idx0] - a[idx2];
1629        x1i = a[idx0 + 1] - a[idx2 + 1];
1630        x2r = a[idx1] + a[idx3];
1631        x2i = a[idx1 + 1] + a[idx3 + 1];
1632        x3r = a[idx1] - a[idx3];
1633        x3i = a[idx1 + 1] - a[idx3 + 1];
1634        a[idx0] = x0r + x2r;
1635        a[idx0 + 1] = x0i + x2i;
1636        a[idx1] = x0r - x2r;
1637        a[idx1 + 1] = x0i - x2i;
1638        x0r = x1r - x3i;
1639        x0i = x1i + x3r;
1640        a[idx2] = wn4r * (x0r - x0i);
1641        a[idx2 + 1] = wn4r * (x0i + x0r);
1642        x0r = x1r + x3i;
1643        x0i = x1i - x3r;
1644        a[idx3] = -wn4r * (x0r + x0i);
1645        a[idx3 + 1] = -wn4r * (x0i - x0r);
1646        x0r = a[idx0 + 2] + a[idx2 + 2];
1647        x0i = a[idx0 + 3] + a[idx2 + 3];
1648        x1r = a[idx0 + 2] - a[idx2 + 2];
1649        x1i = a[idx0 + 3] - a[idx2 + 3];
1650        x2r = a[idx1 + 2] + a[idx3 + 2];
1651        x2i = a[idx1 + 3] + a[idx3 + 3];
1652        x3r = a[idx1 + 2] - a[idx3 + 2];
1653        x3i = a[idx1 + 3] - a[idx3 + 3];
1654        a[idx0 + 2] = x0r + x2r;
1655        a[idx0 + 3] = x0i + x2i;
1656        a[idx1 + 2] = x0r - x2r;
1657        a[idx1 + 3] = x0i - x2i;
1658        x0r = x1r - x3i;
1659        x0i = x1i + x3r;
1660        a[idx2 + 2] = wk1i * x0r - wk1r * x0i;
1661        a[idx2 + 3] = wk1i * x0i + wk1r * x0r;
1662        x0r = x1r + x3i;
1663        x0i = x1i - x3r;
1664        a[idx3 + 2] = wk3i * x0r + wk3r * x0i;
1665        a[idx3 + 3] = wk3i * x0i - wk3r * x0r;
1666    }
1667
1668    private void cftb1st(int n, double[] a, int offa, double[] w, int startw) {
1669        int j0, j1, j2, j3, k, m, mh;
1670        double wn4r, csc1, csc3, wk1r, wk1i, wk3r, wk3i, wd1r, wd1i, wd3r, wd3i;
1671        double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i, y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i;
1672        int idx0, idx1, idx2, idx3, idx4, idx5;
1673        mh = n >> 3;
1674        m = 2 * mh;
1675        j1 = m;
1676        j2 = j1 + m;
1677        j3 = j2 + m;
1678        idx1 = offa + j1;
1679        idx2 = offa + j2;
1680        idx3 = offa + j3;
1681
1682        x0r = a[offa] + a[idx2];
1683        x0i = -a[offa + 1] - a[idx2 + 1];
1684        x1r = a[offa] - a[idx2];
1685        x1i = -a[offa + 1] + a[idx2 + 1];
1686        x2r = a[idx1] + a[idx3];
1687        x2i = a[idx1 + 1] + a[idx3 + 1];
1688        x3r = a[idx1] - a[idx3];
1689        x3i = a[idx1 + 1] - a[idx3 + 1];
1690        a[offa] = x0r + x2r;
1691        a[offa + 1] = x0i - x2i;
1692        a[idx1] = x0r - x2r;
1693        a[idx1 + 1] = x0i + x2i;
1694        a[idx2] = x1r + x3i;
1695        a[idx2 + 1] = x1i + x3r;
1696        a[idx3] = x1r - x3i;
1697        a[idx3 + 1] = x1i - x3r;
1698        wn4r = w[startw + 1];
1699        csc1 = w[startw + 2];
1700        csc3 = w[startw + 3];
1701        wd1r = 1;
1702        wd1i = 0;
1703        wd3r = 1;
1704        wd3i = 0;
1705        k = 0;
1706        for (int j = 2; j < mh - 2; j += 4) {
1707            k += 4;
1708            idx4 = startw + k;
1709            wk1r = csc1 * (wd1r + w[idx4]);
1710            wk1i = csc1 * (wd1i + w[idx4 + 1]);
1711            wk3r = csc3 * (wd3r + w[idx4 + 2]);
1712            wk3i = csc3 * (wd3i + w[idx4 + 3]);
1713            wd1r = w[idx4];
1714            wd1i = w[idx4 + 1];
1715            wd3r = w[idx4 + 2];
1716            wd3i = w[idx4 + 3];
1717            j1 = j + m;
1718            j2 = j1 + m;
1719            j3 = j2 + m;
1720            idx1 = offa + j1;
1721            idx2 = offa + j2;
1722            idx3 = offa + j3;
1723            idx5 = offa + j;
1724            x0r = a[idx5] + a[idx2];
1725            x0i = -a[idx5 + 1] - a[idx2 + 1];
1726            x1r = a[idx5] - a[offa + j2];
1727            x1i = -a[idx5 + 1] + a[idx2 + 1];
1728            y0r = a[idx5 + 2] + a[idx2 + 2];
1729            y0i = -a[idx5 + 3] - a[idx2 + 3];
1730            y1r = a[idx5 + 2] - a[idx2 + 2];
1731            y1i = -a[idx5 + 3] + a[idx2 + 3];
1732            x2r = a[idx1] + a[idx3];
1733            x2i = a[idx1 + 1] + a[idx3 + 1];
1734            x3r = a[idx1] - a[idx3];
1735            x3i = a[idx1 + 1] - a[idx3 + 1];
1736            y2r = a[idx1 + 2] + a[idx3 + 2];
1737            y2i = a[idx1 + 3] + a[idx3 + 3];
1738            y3r = a[idx1 + 2] - a[idx3 + 2];
1739            y3i = a[idx1 + 3] - a[idx3 + 3];
1740            a[idx5] = x0r + x2r;
1741            a[idx5 + 1] = x0i - x2i;
1742            a[idx5 + 2] = y0r + y2r;
1743            a[idx5 + 3] = y0i - y2i;
1744            a[idx1] = x0r - x2r;
1745            a[idx1 + 1] = x0i + x2i;
1746            a[idx1 + 2] = y0r - y2r;
1747            a[idx1 + 3] = y0i + y2i;
1748            x0r = x1r + x3i;
1749            x0i = x1i + x3r;
1750            a[idx2] = wk1r * x0r - wk1i * x0i;
1751            a[idx2 + 1] = wk1r * x0i + wk1i * x0r;
1752            x0r = y1r + y3i;
1753            x0i = y1i + y3r;
1754            a[idx2 + 2] = wd1r * x0r - wd1i * x0i;
1755            a[idx2 + 3] = wd1r * x0i + wd1i * x0r;
1756            x0r = x1r - x3i;
1757            x0i = x1i - x3r;
1758            a[idx3] = wk3r * x0r + wk3i * x0i;
1759            a[idx3 + 1] = wk3r * x0i - wk3i * x0r;
1760            x0r = y1r - y3i;
1761            x0i = y1i - y3r;
1762            a[idx3 + 2] = wd3r * x0r + wd3i * x0i;
1763            a[idx3 + 3] = wd3r * x0i - wd3i * x0r;
1764            j0 = m - j;
1765            j1 = j0 + m;
1766            j2 = j1 + m;
1767            j3 = j2 + m;
1768            idx0 = offa + j0;
1769            idx1 = offa + j1;
1770            idx2 = offa + j2;
1771            idx3 = offa + j3;
1772            x0r = a[idx0] + a[idx2];
1773            x0i = -a[idx0 + 1] - a[idx2 + 1];
1774            x1r = a[idx0] - a[idx2];
1775            x1i = -a[idx0 + 1] + a[idx2 + 1];
1776            y0r = a[idx0 - 2] + a[idx2 - 2];
1777            y0i = -a[idx0 - 1] - a[idx2 - 1];
1778            y1r = a[idx0 - 2] - a[idx2 - 2];
1779            y1i = -a[idx0 - 1] + a[idx2 - 1];
1780            x2r = a[idx1] + a[idx3];
1781            x2i = a[idx1 + 1] + a[idx3 + 1];
1782            x3r = a[idx1] - a[idx3];
1783            x3i = a[idx1 + 1] - a[idx3 + 1];
1784            y2r = a[idx1 - 2] + a[idx3 - 2];
1785            y2i = a[idx1 - 1] + a[idx3 - 1];
1786            y3r = a[idx1 - 2] - a[idx3 - 2];
1787            y3i = a[idx1 - 1] - a[idx3 - 1];
1788            a[idx0] = x0r + x2r;
1789            a[idx0 + 1] = x0i - x2i;
1790            a[idx0 - 2] = y0r + y2r;
1791            a[idx0 - 1] = y0i - y2i;
1792            a[idx1] = x0r - x2r;
1793            a[idx1 + 1] = x0i + x2i;
1794            a[idx1 - 2] = y0r - y2r;
1795            a[idx1 - 1] = y0i + y2i;
1796            x0r = x1r + x3i;
1797            x0i = x1i + x3r;
1798            a[idx2] = wk1i * x0r - wk1r * x0i;
1799            a[idx2 + 1] = wk1i * x0i + wk1r * x0r;
1800            x0r = y1r + y3i;
1801            x0i = y1i + y3r;
1802            a[idx2 - 2] = wd1i * x0r - wd1r * x0i;
1803            a[idx2 - 1] = wd1i * x0i + wd1r * x0r;
1804            x0r = x1r - x3i;
1805            x0i = x1i - x3r;
1806            a[idx3] = wk3i * x0r + wk3r * x0i;
1807            a[idx3 + 1] = wk3i * x0i - wk3r * x0r;
1808            x0r = y1r - y3i;
1809            x0i = y1i - y3r;
1810            a[idx3 - 2] = wd3i * x0r + wd3r * x0i;
1811            a[idx3 - 1] = wd3i * x0i - wd3r * x0r;
1812        }
1813        wk1r = csc1 * (wd1r + wn4r);
1814        wk1i = csc1 * (wd1i + wn4r);
1815        wk3r = csc3 * (wd3r - wn4r);
1816        wk3i = csc3 * (wd3i - wn4r);
1817        j0 = mh;
1818        j1 = j0 + m;
1819        j2 = j1 + m;
1820        j3 = j2 + m;
1821        idx0 = offa + j0;
1822        idx1 = offa + j1;
1823        idx2 = offa + j2;
1824        idx3 = offa + j3;
1825        x0r = a[idx0 - 2] + a[idx2 - 2];
1826        x0i = -a[idx0 - 1] - a[idx2 - 1];
1827        x1r = a[idx0 - 2] - a[idx2 - 2];
1828        x1i = -a[idx0 - 1] + a[idx2 - 1];
1829        x2r = a[idx1 - 2] + a[idx3 - 2];
1830        x2i = a[idx1 - 1] + a[idx3 - 1];
1831        x3r = a[idx1 - 2] - a[idx3 - 2];
1832        x3i = a[idx1 - 1] - a[idx3 - 1];
1833        a[idx0 - 2] = x0r + x2r;
1834        a[idx0 - 1] = x0i - x2i;
1835        a[idx1 - 2] = x0r - x2r;
1836        a[idx1 - 1] = x0i + x2i;
1837        x0r = x1r + x3i;
1838        x0i = x1i + x3r;
1839        a[idx2 - 2] = wk1r * x0r - wk1i * x0i;
1840        a[idx2 - 1] = wk1r * x0i + wk1i * x0r;
1841        x0r = x1r - x3i;
1842        x0i = x1i - x3r;
1843        a[idx3 - 2] = wk3r * x0r + wk3i * x0i;
1844        a[idx3 - 1] = wk3r * x0i - wk3i * x0r;
1845        x0r = a[idx0] + a[idx2];
1846        x0i = -a[idx0 + 1] - a[idx2 + 1];
1847        x1r = a[idx0] - a[idx2];
1848        x1i = -a[idx0 + 1] + a[idx2 + 1];
1849        x2r = a[idx1] + a[idx3];
1850        x2i = a[idx1 + 1] + a[idx3 + 1];
1851        x3r = a[idx1] - a[idx3];
1852        x3i = a[idx1 + 1] - a[idx3 + 1];
1853        a[idx0] = x0r + x2r;
1854        a[idx0 + 1] = x0i - x2i;
1855        a[idx1] = x0r - x2r;
1856        a[idx1 + 1] = x0i + x2i;
1857        x0r = x1r + x3i;
1858        x0i = x1i + x3r;
1859        a[idx2] = wn4r * (x0r - x0i);
1860        a[idx2 + 1] = wn4r * (x0i + x0r);
1861        x0r = x1r - x3i;
1862        x0i = x1i - x3r;
1863        a[idx3] = -wn4r * (x0r + x0i);
1864        a[idx3 + 1] = -wn4r * (x0i - x0r);
1865        x0r = a[idx0 + 2] + a[idx2 + 2];
1866        x0i = -a[idx0 + 3] - a[idx2 + 3];
1867        x1r = a[idx0 + 2] - a[idx2 + 2];
1868        x1i = -a[idx0 + 3] + a[idx2 + 3];
1869        x2r = a[idx1 + 2] + a[idx3 + 2];
1870        x2i = a[idx1 + 3] + a[idx3 + 3];
1871        x3r = a[idx1 + 2] - a[idx3 + 2];
1872        x3i = a[idx1 + 3] - a[idx3 + 3];
1873        a[idx0 + 2] = x0r + x2r;
1874        a[idx0 + 3] = x0i - x2i;
1875        a[idx1 + 2] = x0r - x2r;
1876        a[idx1 + 3] = x0i + x2i;
1877        x0r = x1r + x3i;
1878        x0i = x1i + x3r;
1879        a[idx2 + 2] = wk1i * x0r - wk1r * x0i;
1880        a[idx2 + 3] = wk1i * x0i + wk1r * x0r;
1881        x0r = x1r - x3i;
1882        x0i = x1i - x3r;
1883        a[idx3 + 2] = wk3i * x0r + wk3r * x0i;
1884        a[idx3 + 3] = wk3i * x0i - wk3r * x0r;
1885    }
1886
1887    private void cftrec4_th(final int n, final double[] a, final int offa, final int nw, final double[] w) {
1888        int idiv4, m, nthread;
1889        int idx = 0;
1890        nthread = 2;
1891        idiv4 = 0;
1892        m = n >> 1;
1893        if (n > ConcurrencyUtils.getThreadsBeginN_1D_FFT_4Threads()) {
1894            nthread = 4;
1895            idiv4 = 1;
1896            m >>= 1;
1897        }
1898        Future<?>[] futures = new Future[nthread];
1899        final int mf = m;
1900        for (int i = 0; i < nthread; i++) {
1901            final int firstIdx = offa + i * m;
1902            if (i != idiv4) {
1903                futures[idx++] = ConcurrencyUtils.submit(new Runnable() {
1904                    @Override
1905                                        public void run() {
1906                        int isplt, k, m;
1907                        int idx1 = firstIdx + mf;
1908                        m = n;
1909                        while (m > 512) {
1910                            m >>= 2;
1911                            cftmdl1(m, a, idx1 - m, w, nw - (m >> 1));
1912                        }
1913                        cftleaf(m, 1, a, idx1 - m, nw, w);
1914                        k = 0;
1915                        int idx2 = firstIdx - m;
1916                        for (int j = mf - m; j > 0; j -= m) {
1917                            k++;
1918                            isplt = cfttree(m, j, k, a, firstIdx, nw, w);
1919                            cftleaf(m, isplt, a, idx2 + j, nw, w);
1920                        }
1921                    }
1922                });
1923            } else {
1924                futures[idx++] = ConcurrencyUtils.submit(new Runnable() {
1925                    @Override
1926                                        public void run() {
1927                        int isplt, k, m;
1928                        int idx1 = firstIdx + mf;
1929                        k = 1;
1930                        m = n;
1931                        while (m > 512) {
1932                            m >>= 2;
1933                            k <<= 2;
1934                            cftmdl2(m, a, idx1 - m, w, nw - m);
1935                        }
1936                        cftleaf(m, 0, a, idx1 - m, nw, w);
1937                        k >>= 1;
1938                        int idx2 = firstIdx - m;
1939                        for (int j = mf - m; j > 0; j -= m) {
1940                            k++;
1941                            isplt = cfttree(m, j, k, a, firstIdx, nw, w);
1942                            cftleaf(m, isplt, a, idx2 + j, nw, w);
1943                        }
1944                    }
1945                });
1946            }
1947        }
1948        ConcurrencyUtils.waitForCompletion(futures);
1949    }
1950
1951    private void cftrec4(int n, double[] a, int offa, int nw, double[] w) {
1952        int isplt, k, m;
1953        int idx1 = offa + n;
1954        m = n;
1955        while (m > 512) {
1956            m >>= 2;
1957            cftmdl1(m, a, idx1 - m, w, nw - (m >> 1));
1958        }
1959        cftleaf(m, 1, a, idx1 - m, nw, w);
1960        k = 0;
1961        int idx2 = offa - m;
1962        for (int j = n - m; j > 0; j -= m) {
1963            k++;
1964            isplt = cfttree(m, j, k, a, offa, nw, w);
1965            cftleaf(m, isplt, a, idx2 + j, nw, w);
1966        }
1967    }
1968
1969    private int cfttree(int n, int j, int k, double[] a, int offa, int nw, double[] w) {
1970        int i, isplt, m;
1971        int idx1 = offa - n;
1972        if ((k & 3) != 0) {
1973            isplt = k & 1;
1974            if (isplt != 0) {
1975                cftmdl1(n, a, idx1 + j, w, nw - (n >> 1));
1976            } else {
1977                cftmdl2(n, a, idx1 + j, w, nw - n);
1978            }
1979        } else {
1980            m = n;
1981            for (i = k; (i & 3) == 0; i >>= 2) {
1982                m <<= 2;
1983            }
1984            isplt = i & 1;
1985            int idx2 = offa + j;
1986            if (isplt != 0) {
1987                while (m > 128) {
1988                    cftmdl1(m, a, idx2 - m, w, nw - (m >> 1));
1989                    m >>= 2;
1990                }
1991            } else {
1992                while (m > 128) {
1993                    cftmdl2(m, a, idx2 - m, w, nw - m);
1994                    m >>= 2;
1995                }
1996            }
1997        }
1998        return isplt;
1999    }
2000
2001    private void cftleaf(int n, int isplt, double[] a, int offa, int nw, double[] w) {
2002        if (n == 512) {
2003            cftmdl1(128, a, offa, w, nw - 64);
2004            cftf161(a, offa, w, nw - 8);
2005            cftf162(a, offa + 32, w, nw - 32);
2006            cftf161(a, offa + 64, w, nw - 8);
2007            cftf161(a, offa + 96, w, nw - 8);
2008            cftmdl2(128, a, offa + 128, w, nw - 128);
2009            cftf161(a, offa + 128, w, nw - 8);
2010            cftf162(a, offa + 160, w, nw - 32);
2011            cftf161(a, offa + 192, w, nw - 8);
2012            cftf162(a, offa + 224, w, nw - 32);
2013            cftmdl1(128, a, offa + 256, w, nw - 64);
2014            cftf161(a, offa + 256, w, nw - 8);
2015            cftf162(a, offa + 288, w, nw - 32);
2016            cftf161(a, offa + 320, w, nw - 8);
2017            cftf161(a, offa + 352, w, nw - 8);
2018            if (isplt != 0) {
2019                cftmdl1(128, a, offa + 384, w, nw - 64);
2020                cftf161(a, offa + 480, w, nw - 8);
2021            } else {
2022                cftmdl2(128, a, offa + 384, w, nw - 128);
2023                cftf162(a, offa + 480, w, nw - 32);
2024            }
2025            cftf161(a, offa + 384, w, nw - 8);
2026            cftf162(a, offa + 416, w, nw - 32);
2027            cftf161(a, offa + 448, w, nw - 8);
2028        } else {
2029            cftmdl1(64, a, offa, w, nw - 32);
2030            cftf081(a, offa, w, nw - 8);
2031            cftf082(a, offa + 16, w, nw - 8);
2032            cftf081(a, offa + 32, w, nw - 8);
2033            cftf081(a, offa + 48, w, nw - 8);
2034            cftmdl2(64, a, offa + 64, w, nw - 64);
2035            cftf081(a, offa + 64, w, nw - 8);
2036            cftf082(a, offa + 80, w, nw - 8);
2037            cftf081(a, offa + 96, w, nw - 8);
2038            cftf082(a, offa + 112, w, nw - 8);
2039            cftmdl1(64, a, offa + 128, w, nw - 32);
2040            cftf081(a, offa + 128, w, nw - 8);
2041            cftf082(a, offa + 144, w, nw - 8);
2042            cftf081(a, offa + 160, w, nw - 8);
2043            cftf081(a, offa + 176, w, nw - 8);
2044            if (isplt != 0) {
2045                cftmdl1(64, a, offa + 192, w, nw - 32);
2046                cftf081(a, offa + 240, w, nw - 8);
2047            } else {
2048                cftmdl2(64, a, offa + 192, w, nw - 64);
2049                cftf082(a, offa + 240, w, nw - 8);
2050            }
2051            cftf081(a, offa + 192, w, nw - 8);
2052            cftf082(a, offa + 208, w, nw - 8);
2053            cftf081(a, offa + 224, w, nw - 8);
2054        }
2055    }
2056
2057    private void cftmdl1(int n, double[] a, int offa, double[] w, int startw) {
2058        int j0, j1, j2, j3, k, m, mh;
2059        double wn4r, wk1r, wk1i, wk3r, wk3i;
2060        double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
2061        int idx0, idx1, idx2, idx3, idx4, idx5;
2062
2063        mh = n >> 3;
2064        m = 2 * mh;
2065        j1 = m;
2066        j2 = j1 + m;
2067        j3 = j2 + m;
2068        idx1 = offa + j1;
2069        idx2 = offa + j2;
2070        idx3 = offa + j3;
2071        x0r = a[offa] + a[idx2];
2072        x0i = a[offa + 1] + a[idx2 + 1];
2073        x1r = a[offa] - a[idx2];
2074        x1i = a[offa + 1] - a[idx2 + 1];
2075        x2r = a[idx1] + a[idx3];
2076        x2i = a[idx1 + 1] + a[idx3 + 1];
2077        x3r = a[idx1] - a[idx3];
2078        x3i = a[idx1 + 1] - a[idx3 + 1];
2079        a[offa] = x0r + x2r;
2080        a[offa + 1] = x0i + x2i;
2081        a[idx1] = x0r - x2r;
2082        a[idx1 + 1] = x0i - x2i;
2083        a[idx2] = x1r - x3i;
2084        a[idx2 + 1] = x1i + x3r;
2085        a[idx3] = x1r + x3i;
2086        a[idx3 + 1] = x1i - x3r;
2087        wn4r = w[startw + 1];
2088        k = 0;
2089        for (int j = 2; j < mh; j += 2) {
2090            k += 4;
2091            idx4 = startw + k;
2092            wk1r = w[idx4];
2093            wk1i = w[idx4 + 1];
2094            wk3r = w[idx4 + 2];
2095            wk3i = w[idx4 + 3];
2096            j1 = j + m;
2097            j2 = j1 + m;
2098            j3 = j2 + m;
2099            idx1 = offa + j1;
2100            idx2 = offa + j2;
2101            idx3 = offa + j3;
2102            idx5 = offa + j;
2103            x0r = a[idx5] + a[idx2];
2104            x0i = a[idx5 + 1] + a[idx2 + 1];
2105            x1r = a[idx5] - a[idx2];
2106            x1i = a[idx5 + 1] - a[idx2 + 1];
2107            x2r = a[idx1] + a[idx3];
2108            x2i = a[idx1 + 1] + a[idx3 + 1];
2109            x3r = a[idx1] - a[idx3];
2110            x3i = a[idx1 + 1] - a[idx3 + 1];
2111            a[idx5] = x0r + x2r;
2112            a[idx5 + 1] = x0i + x2i;
2113            a[idx1] = x0r - x2r;
2114            a[idx1 + 1] = x0i - x2i;
2115            x0r = x1r - x3i;
2116            x0i = x1i + x3r;
2117            a[idx2] = wk1r * x0r - wk1i * x0i;
2118            a[idx2 + 1] = wk1r * x0i + wk1i * x0r;
2119            x0r = x1r + x3i;
2120            x0i = x1i - x3r;
2121            a[idx3] = wk3r * x0r + wk3i * x0i;
2122            a[idx3 + 1] = wk3r * x0i - wk3i * x0r;
2123            j0 = m - j;
2124            j1 = j0 + m;
2125            j2 = j1 + m;
2126            j3 = j2 + m;
2127            idx0 = offa + j0;
2128            idx1 = offa + j1;
2129            idx2 = offa + j2;
2130            idx3 = offa + j3;
2131            x0r = a[idx0] + a[idx2];
2132            x0i = a[idx0 + 1] + a[idx2 + 1];
2133            x1r = a[idx0] - a[idx2];
2134            x1i = a[idx0 + 1] - a[idx2 + 1];
2135            x2r = a[idx1] + a[idx3];
2136            x2i = a[idx1 + 1] + a[idx3 + 1];
2137            x3r = a[idx1] - a[idx3];
2138            x3i = a[idx1 + 1] - a[idx3 + 1];
2139            a[idx0] = x0r + x2r;
2140            a[idx0 + 1] = x0i + x2i;
2141            a[idx1] = x0r - x2r;
2142            a[idx1 + 1] = x0i - x2i;
2143            x0r = x1r - x3i;
2144            x0i = x1i + x3r;
2145            a[idx2] = wk1i * x0r - wk1r * x0i;
2146            a[idx2 + 1] = wk1i * x0i + wk1r * x0r;
2147            x0r = x1r + x3i;
2148            x0i = x1i - x3r;
2149            a[idx3] = wk3i * x0r + wk3r * x0i;
2150            a[idx3 + 1] = wk3i * x0i - wk3r * x0r;
2151        }
2152        j0 = mh;
2153        j1 = j0 + m;
2154        j2 = j1 + m;
2155        j3 = j2 + m;
2156        idx0 = offa + j0;
2157        idx1 = offa + j1;
2158        idx2 = offa + j2;
2159        idx3 = offa + j3;
2160        x0r = a[idx0] + a[idx2];
2161        x0i = a[idx0 + 1] + a[idx2 + 1];
2162        x1r = a[idx0] - a[idx2];
2163        x1i = a[idx0 + 1] - a[idx2 + 1];
2164        x2r = a[idx1] + a[idx3];
2165        x2i = a[idx1 + 1] + a[idx3 + 1];
2166        x3r = a[idx1] - a[idx3];
2167        x3i = a[idx1 + 1] - a[idx3 + 1];
2168        a[idx0] = x0r + x2r;
2169        a[idx0 + 1] = x0i + x2i;
2170        a[idx1] = x0r - x2r;
2171        a[idx1 + 1] = x0i - x2i;
2172        x0r = x1r - x3i;
2173        x0i = x1i + x3r;
2174        a[idx2] = wn4r * (x0r - x0i);
2175        a[idx2 + 1] = wn4r * (x0i + x0r);
2176        x0r = x1r + x3i;
2177        x0i = x1i - x3r;
2178        a[idx3] = -wn4r * (x0r + x0i);
2179        a[idx3 + 1] = -wn4r * (x0i - x0r);
2180    }
2181
2182    private void cftmdl2(int n, double[] a, int offa, double[] w, int startw) {
2183        int j0, j1, j2, j3, k, kr, m, mh;
2184        double wn4r, wk1r, wk1i, wk3r, wk3i, wd1r, wd1i, wd3r, wd3i;
2185        double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i, y0r, y0i, y2r, y2i;
2186        int idx0, idx1, idx2, idx3, idx4, idx5, idx6;
2187
2188        mh = n >> 3;
2189        m = 2 * mh;
2190        wn4r = w[startw + 1];
2191        j1 = m;
2192        j2 = j1 + m;
2193        j3 = j2 + m;
2194        idx1 = offa + j1;
2195        idx2 = offa + j2;
2196        idx3 = offa + j3;
2197        x0r = a[offa] - a[idx2 + 1];
2198        x0i = a[offa + 1] + a[idx2];
2199        x1r = a[offa] + a[idx2 + 1];
2200        x1i = a[offa + 1] - a[idx2];
2201        x2r = a[idx1] - a[idx3 + 1];
2202        x2i = a[idx1 + 1] + a[idx3];
2203        x3r = a[idx1] + a[idx3 + 1];
2204        x3i = a[idx1 + 1] - a[idx3];
2205        y0r = wn4r * (x2r - x2i);
2206        y0i = wn4r * (x2i + x2r);
2207        a[offa] = x0r + y0r;
2208        a[offa + 1] = x0i + y0i;
2209        a[idx1] = x0r - y0r;
2210        a[idx1 + 1] = x0i - y0i;
2211        y0r = wn4r * (x3r - x3i);
2212        y0i = wn4r * (x3i + x3r);
2213        a[idx2] = x1r - y0i;
2214        a[idx2 + 1] = x1i + y0r;
2215        a[idx3] = x1r + y0i;
2216        a[idx3 + 1] = x1i - y0r;
2217        k = 0;
2218        kr = 2 * m;
2219        for (int j = 2; j < mh; j += 2) {
2220            k += 4;
2221            idx4 = startw + k;
2222            wk1r = w[idx4];
2223            wk1i = w[idx4 + 1];
2224            wk3r = w[idx4 + 2];
2225            wk3i = w[idx4 + 3];
2226            kr -= 4;
2227            idx5 = startw + kr;
2228            wd1i = w[idx5];
2229            wd1r = w[idx5 + 1];
2230            wd3i = w[idx5 + 2];
2231            wd3r = w[idx5 + 3];
2232            j1 = j + m;
2233            j2 = j1 + m;
2234            j3 = j2 + m;
2235            idx1 = offa + j1;
2236            idx2 = offa + j2;
2237            idx3 = offa + j3;
2238            idx6 = offa + j;
2239            x0r = a[idx6] - a[idx2 + 1];
2240            x0i = a[idx6 + 1] + a[idx2];
2241            x1r = a[idx6] + a[idx2 + 1];
2242            x1i = a[idx6 + 1] - a[idx2];
2243            x2r = a[idx1] - a[idx3 + 1];
2244            x2i = a[idx1 + 1] + a[idx3];
2245            x3r = a[idx1] + a[idx3 + 1];
2246            x3i = a[idx1 + 1] - a[idx3];
2247            y0r = wk1r * x0r - wk1i * x0i;
2248            y0i = wk1r * x0i + wk1i * x0r;
2249            y2r = wd1r * x2r - wd1i * x2i;
2250            y2i = wd1r * x2i + wd1i * x2r;
2251            a[idx6] = y0r + y2r;
2252            a[idx6 + 1] = y0i + y2i;
2253            a[idx1] = y0r - y2r;
2254            a[idx1 + 1] = y0i - y2i;
2255            y0r = wk3r * x1r + wk3i * x1i;
2256            y0i = wk3r * x1i - wk3i * x1r;
2257            y2r = wd3r * x3r + wd3i * x3i;
2258            y2i = wd3r * x3i - wd3i * x3r;
2259            a[idx2] = y0r + y2r;
2260            a[idx2 + 1] = y0i + y2i;
2261            a[idx3] = y0r - y2r;
2262            a[idx3 + 1] = y0i - y2i;
2263            j0 = m - j;
2264            j1 = j0 + m;
2265            j2 = j1 + m;
2266            j3 = j2 + m;
2267            idx0 = offa + j0;
2268            idx1 = offa + j1;
2269            idx2 = offa + j2;
2270            idx3 = offa + j3;
2271            x0r = a[idx0] - a[idx2 + 1];
2272            x0i = a[idx0 + 1] + a[idx2];
2273            x1r = a[idx0] + a[idx2 + 1];
2274            x1i = a[idx0 + 1] - a[idx2];
2275            x2r = a[idx1] - a[idx3 + 1];
2276            x2i = a[idx1 + 1] + a[idx3];
2277            x3r = a[idx1] + a[idx3 + 1];
2278            x3i = a[idx1 + 1] - a[idx3];
2279            y0r = wd1i * x0r - wd1r * x0i;
2280            y0i = wd1i * x0i + wd1r * x0r;
2281            y2r = wk1i * x2r - wk1r * x2i;
2282            y2i = wk1i * x2i + wk1r * x2r;
2283            a[idx0] = y0r + y2r;
2284            a[idx0 + 1] = y0i + y2i;
2285            a[idx1] = y0r - y2r;
2286            a[idx1 + 1] = y0i - y2i;
2287            y0r = wd3i * x1r + wd3r * x1i;
2288            y0i = wd3i * x1i - wd3r * x1r;
2289            y2r = wk3i * x3r + wk3r * x3i;
2290            y2i = wk3i * x3i - wk3r * x3r;
2291            a[idx2] = y0r + y2r;
2292            a[idx2 + 1] = y0i + y2i;
2293            a[idx3] = y0r - y2r;
2294            a[idx3 + 1] = y0i - y2i;
2295        }
2296        wk1r = w[startw + m];
2297        wk1i = w[startw + m + 1];
2298        j0 = mh;
2299        j1 = j0 + m;
2300        j2 = j1 + m;
2301        j3 = j2 + m;
2302        idx0 = offa + j0;
2303        idx1 = offa + j1;
2304        idx2 = offa + j2;
2305        idx3 = offa + j3;
2306        x0r = a[idx0] - a[idx2 + 1];
2307        x0i = a[idx0 + 1] + a[idx2];
2308        x1r = a[idx0] + a[idx2 + 1];
2309        x1i = a[idx0 + 1] - a[idx2];
2310        x2r = a[idx1] - a[idx3 + 1];
2311        x2i = a[idx1 + 1] + a[idx3];
2312        x3r = a[idx1] + a[idx3 + 1];
2313        x3i = a[idx1 + 1] - a[idx3];
2314        y0r = wk1r * x0r - wk1i * x0i;
2315        y0i = wk1r * x0i + wk1i * x0r;
2316        y2r = wk1i * x2r - wk1r * x2i;
2317        y2i = wk1i * x2i + wk1r * x2r;
2318        a[idx0] = y0r + y2r;
2319        a[idx0 + 1] = y0i + y2i;
2320        a[idx1] = y0r - y2r;
2321        a[idx1 + 1] = y0i - y2i;
2322        y0r = wk1i * x1r - wk1r * x1i;
2323        y0i = wk1i * x1i + wk1r * x1r;
2324        y2r = wk1r * x3r - wk1i * x3i;
2325        y2i = wk1r * x3i + wk1i * x3r;
2326        a[idx2] = y0r - y2r;
2327        a[idx2 + 1] = y0i - y2i;
2328        a[idx3] = y0r + y2r;
2329        a[idx3 + 1] = y0i + y2i;
2330    }
2331
2332    private void cftfx41(int n, double[] a, int offa, int nw, double[] w) {
2333        if (n == 128) {
2334            cftf161(a, offa, w, nw - 8);
2335            cftf162(a, offa + 32, w, nw - 32);
2336            cftf161(a, offa + 64, w, nw - 8);
2337            cftf161(a, offa + 96, w, nw - 8);
2338        } else {
2339            cftf081(a, offa, w, nw - 8);
2340            cftf082(a, offa + 16, w, nw - 8);
2341            cftf081(a, offa + 32, w, nw - 8);
2342            cftf081(a, offa + 48, w, nw - 8);
2343        }
2344    }
2345
2346    private void cftf161(double[] a, int offa, double[] w, int startw) {
2347        double wn4r, wk1r, wk1i, x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i, y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i, y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i, y8r, y8i, y9r, y9i, y10r, y10i, y11r, y11i, y12r, y12i, y13r, y13i, y14r, y14i, y15r, y15i;
2348
2349        wn4r = w[startw + 1];
2350        wk1r = w[startw + 2];
2351        wk1i = w[startw + 3];
2352        x0r = a[offa] + a[offa + 16];
2353        x0i = a[offa + 1] + a[offa + 17];
2354        x1r = a[offa] - a[offa + 16];
2355        x1i = a[offa + 1] - a[offa + 17];
2356        x2r = a[offa + 8] + a[offa + 24];
2357        x2i = a[offa + 9] + a[offa + 25];
2358        x3r = a[offa + 8] - a[offa + 24];
2359        x3i = a[offa + 9] - a[offa + 25];
2360        y0r = x0r + x2r;
2361        y0i = x0i + x2i;
2362        y4r = x0r - x2r;
2363        y4i = x0i - x2i;
2364        y8r = x1r - x3i;
2365        y8i = x1i + x3r;
2366        y12r = x1r + x3i;
2367        y12i = x1i - x3r;
2368        x0r = a[offa + 2] + a[offa + 18];
2369        x0i = a[offa + 3] + a[offa + 19];
2370        x1r = a[offa + 2] - a[offa + 18];
2371        x1i = a[offa + 3] - a[offa + 19];
2372        x2r = a[offa + 10] + a[offa + 26];
2373        x2i = a[offa + 11] + a[offa + 27];
2374        x3r = a[offa + 10] - a[offa + 26];
2375        x3i = a[offa + 11] - a[offa + 27];
2376        y1r = x0r + x2r;
2377        y1i = x0i + x2i;
2378        y5r = x0r - x2r;
2379        y5i = x0i - x2i;
2380        x0r = x1r - x3i;
2381        x0i = x1i + x3r;
2382        y9r = wk1r * x0r - wk1i * x0i;
2383        y9i = wk1r * x0i + wk1i * x0r;
2384        x0r = x1r + x3i;
2385        x0i = x1i - x3r;
2386        y13r = wk1i * x0r - wk1r * x0i;
2387        y13i = wk1i * x0i + wk1r * x0r;
2388        x0r = a[offa + 4] + a[offa + 20];
2389        x0i = a[offa + 5] + a[offa + 21];
2390        x1r = a[offa + 4] - a[offa + 20];
2391        x1i = a[offa + 5] - a[offa + 21];
2392        x2r = a[offa + 12] + a[offa + 28];
2393        x2i = a[offa + 13] + a[offa + 29];
2394        x3r = a[offa + 12] - a[offa + 28];
2395        x3i = a[offa + 13] - a[offa + 29];
2396        y2r = x0r + x2r;
2397        y2i = x0i + x2i;
2398        y6r = x0r - x2r;
2399        y6i = x0i - x2i;
2400        x0r = x1r - x3i;
2401        x0i = x1i + x3r;
2402        y10r = wn4r * (x0r - x0i);
2403        y10i = wn4r * (x0i + x0r);
2404        x0r = x1r + x3i;
2405        x0i = x1i - x3r;
2406        y14r = wn4r * (x0r + x0i);
2407        y14i = wn4r * (x0i - x0r);
2408        x0r = a[offa + 6] + a[offa + 22];
2409        x0i = a[offa + 7] + a[offa + 23];
2410        x1r = a[offa + 6] - a[offa + 22];
2411        x1i = a[offa + 7] - a[offa + 23];
2412        x2r = a[offa + 14] + a[offa + 30];
2413        x2i = a[offa + 15] + a[offa + 31];
2414        x3r = a[offa + 14] - a[offa + 30];
2415        x3i = a[offa + 15] - a[offa + 31];
2416        y3r = x0r + x2r;
2417        y3i = x0i + x2i;
2418        y7r = x0r - x2r;
2419        y7i = x0i - x2i;
2420        x0r = x1r - x3i;
2421        x0i = x1i + x3r;
2422        y11r = wk1i * x0r - wk1r * x0i;
2423        y11i = wk1i * x0i + wk1r * x0r;
2424        x0r = x1r + x3i;
2425        x0i = x1i - x3r;
2426        y15r = wk1r * x0r - wk1i * x0i;
2427        y15i = wk1r * x0i + wk1i * x0r;
2428        x0r = y12r - y14r;
2429        x0i = y12i - y14i;
2430        x1r = y12r + y14r;
2431        x1i = y12i + y14i;
2432        x2r = y13r - y15r;
2433        x2i = y13i - y15i;
2434        x3r = y13r + y15r;
2435        x3i = y13i + y15i;
2436        a[offa + 24] = x0r + x2r;
2437        a[offa + 25] = x0i + x2i;
2438        a[offa + 26] = x0r - x2r;
2439        a[offa + 27] = x0i - x2i;
2440        a[offa + 28] = x1r - x3i;
2441        a[offa + 29] = x1i + x3r;
2442        a[offa + 30] = x1r + x3i;
2443        a[offa + 31] = x1i - x3r;
2444        x0r = y8r + y10r;
2445        x0i = y8i + y10i;
2446        x1r = y8r - y10r;
2447        x1i = y8i - y10i;
2448        x2r = y9r + y11r;
2449        x2i = y9i + y11i;
2450        x3r = y9r - y11r;
2451        x3i = y9i - y11i;
2452        a[offa + 16] = x0r + x2r;
2453        a[offa + 17] = x0i + x2i;
2454        a[offa + 18] = x0r - x2r;
2455        a[offa + 19] = x0i - x2i;
2456        a[offa + 20] = x1r - x3i;
2457        a[offa + 21] = x1i + x3r;
2458        a[offa + 22] = x1r + x3i;
2459        a[offa + 23] = x1i - x3r;
2460        x0r = y5r - y7i;
2461        x0i = y5i + y7r;
2462        x2r = wn4r * (x0r - x0i);
2463        x2i = wn4r * (x0i + x0r);
2464        x0r = y5r + y7i;
2465        x0i = y5i - y7r;
2466        x3r = wn4r * (x0r - x0i);
2467        x3i = wn4r * (x0i + x0r);
2468        x0r = y4r - y6i;
2469        x0i = y4i + y6r;
2470        x1r = y4r + y6i;
2471        x1i = y4i - y6r;
2472        a[offa + 8] = x0r + x2r;
2473        a[offa + 9] = x0i + x2i;
2474        a[offa + 10] = x0r - x2r;
2475        a[offa + 11] = x0i - x2i;
2476        a[offa + 12] = x1r - x3i;
2477        a[offa + 13] = x1i + x3r;
2478        a[offa + 14] = x1r + x3i;
2479        a[offa + 15] = x1i - x3r;
2480        x0r = y0r + y2r;
2481        x0i = y0i + y2i;
2482        x1r = y0r - y2r;
2483        x1i = y0i - y2i;
2484        x2r = y1r + y3r;
2485        x2i = y1i + y3i;
2486        x3r = y1r - y3r;
2487        x3i = y1i - y3i;
2488        a[offa] = x0r + x2r;
2489        a[offa + 1] = x0i + x2i;
2490        a[offa + 2] = x0r - x2r;
2491        a[offa + 3] = x0i - x2i;
2492        a[offa + 4] = x1r - x3i;
2493        a[offa + 5] = x1i + x3r;
2494        a[offa + 6] = x1r + x3i;
2495        a[offa + 7] = x1i - x3r;
2496    }
2497
2498    private void cftf162(double[] a, int offa, double[] w, int startw) {
2499        double wn4r, wk1r, wk1i, wk2r, wk2i, wk3r, wk3i, x0r, x0i, x1r, x1i, x2r, x2i, y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i, y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i, y8r, y8i, y9r, y9i, y10r, y10i, y11r, y11i, y12r, y12i, y13r, y13i, y14r, y14i, y15r, y15i;
2500
2501        wn4r = w[startw + 1];
2502        wk1r = w[startw + 4];
2503        wk1i = w[startw + 5];
2504        wk3r = w[startw + 6];
2505        wk3i = -w[startw + 7];
2506        wk2r = w[startw + 8];
2507        wk2i = w[startw + 9];
2508        x1r = a[offa] - a[offa + 17];
2509        x1i = a[offa + 1] + a[offa + 16];
2510        x0r = a[offa + 8] - a[offa + 25];
2511        x0i = a[offa + 9] + a[offa + 24];
2512        x2r = wn4r * (x0r - x0i);
2513        x2i = wn4r * (x0i + x0r);
2514        y0r = x1r + x2r;
2515        y0i = x1i + x2i;
2516        y4r = x1r - x2r;
2517        y4i = x1i - x2i;
2518        x1r = a[offa] + a[offa + 17];
2519        x1i = a[offa + 1] - a[offa + 16];
2520        x0r = a[offa + 8] + a[offa + 25];
2521        x0i = a[offa + 9] - a[offa + 24];
2522        x2r = wn4r * (x0r - x0i);
2523        x2i = wn4r * (x0i + x0r);
2524        y8r = x1r - x2i;
2525        y8i = x1i + x2r;
2526        y12r = x1r + x2i;
2527        y12i = x1i - x2r;
2528        x0r = a[offa + 2] - a[offa + 19];
2529        x0i = a[offa + 3] + a[offa + 18];
2530        x1r = wk1r * x0r - wk1i * x0i;
2531        x1i = wk1r * x0i + wk1i * x0r;
2532        x0r = a[offa + 10] - a[offa + 27];
2533        x0i = a[offa + 11] + a[offa + 26];
2534        x2r = wk3i * x0r - wk3r * x0i;
2535        x2i = wk3i * x0i + wk3r * x0r;
2536        y1r = x1r + x2r;
2537        y1i = x1i + x2i;
2538        y5r = x1r - x2r;
2539        y5i = x1i - x2i;
2540        x0r = a[offa + 2] + a[offa + 19];
2541        x0i = a[offa + 3] - a[offa + 18];
2542        x1r = wk3r * x0r - wk3i * x0i;
2543        x1i = wk3r * x0i + wk3i * x0r;
2544        x0r = a[offa + 10] + a[offa + 27];
2545        x0i = a[offa + 11] - a[offa + 26];
2546        x2r = wk1r * x0r + wk1i * x0i;
2547        x2i = wk1r * x0i - wk1i * x0r;
2548        y9r = x1r - x2r;
2549        y9i = x1i - x2i;
2550        y13r = x1r + x2r;
2551        y13i = x1i + x2i;
2552        x0r = a[offa + 4] - a[offa + 21];
2553        x0i = a[offa + 5] + a[offa + 20];
2554        x1r = wk2r * x0r - wk2i * x0i;
2555        x1i = wk2r * x0i + wk2i * x0r;
2556        x0r = a[offa + 12] - a[offa + 29];
2557        x0i = a[offa + 13] + a[offa + 28];
2558        x2r = wk2i * x0r - wk2r * x0i;
2559        x2i = wk2i * x0i + wk2r * x0r;
2560        y2r = x1r + x2r;
2561        y2i = x1i + x2i;
2562        y6r = x1r - x2r;
2563        y6i = x1i - x2i;
2564        x0r = a[offa + 4] + a[offa + 21];
2565        x0i = a[offa + 5] - a[offa + 20];
2566        x1r = wk2i * x0r - wk2r * x0i;
2567        x1i = wk2i * x0i + wk2r * x0r;
2568        x0r = a[offa + 12] + a[offa + 29];
2569        x0i = a[offa + 13] - a[offa + 28];
2570        x2r = wk2r * x0r - wk2i * x0i;
2571        x2i = wk2r * x0i + wk2i * x0r;
2572        y10r = x1r - x2r;
2573        y10i = x1i - x2i;
2574        y14r = x1r + x2r;
2575        y14i = x1i + x2i;
2576        x0r = a[offa + 6] - a[offa + 23];
2577        x0i = a[offa + 7] + a[offa + 22];
2578        x1r = wk3r * x0r - wk3i * x0i;
2579        x1i = wk3r * x0i + wk3i * x0r;
2580        x0r = a[offa + 14] - a[offa + 31];
2581        x0i = a[offa + 15] + a[offa + 30];
2582        x2r = wk1i * x0r - wk1r * x0i;
2583        x2i = wk1i * x0i + wk1r * x0r;
2584        y3r = x1r + x2r;
2585        y3i = x1i + x2i;
2586        y7r = x1r - x2r;
2587        y7i = x1i - x2i;
2588        x0r = a[offa + 6] + a[offa + 23];
2589        x0i = a[offa + 7] - a[offa + 22];
2590        x1r = wk1i * x0r + wk1r * x0i;
2591        x1i = wk1i * x0i - wk1r * x0r;
2592        x0r = a[offa + 14] + a[offa + 31];
2593        x0i = a[offa + 15] - a[offa + 30];
2594        x2r = wk3i * x0r - wk3r * x0i;
2595        x2i = wk3i * x0i + wk3r * x0r;
2596        y11r = x1r + x2r;
2597        y11i = x1i + x2i;
2598        y15r = x1r - x2r;
2599        y15i = x1i - x2i;
2600        x1r = y0r + y2r;
2601        x1i = y0i + y2i;
2602        x2r = y1r + y3r;
2603        x2i = y1i + y3i;
2604        a[offa] = x1r + x2r;
2605        a[offa + 1] = x1i + x2i;
2606        a[offa + 2] = x1r - x2r;
2607        a[offa + 3] = x1i - x2i;
2608        x1r = y0r - y2r;
2609        x1i = y0i - y2i;
2610        x2r = y1r - y3r;
2611        x2i = y1i - y3i;
2612        a[offa + 4] = x1r - x2i;
2613        a[offa + 5] = x1i + x2r;
2614        a[offa + 6] = x1r + x2i;
2615        a[offa + 7] = x1i - x2r;
2616        x1r = y4r - y6i;
2617        x1i = y4i + y6r;
2618        x0r = y5r - y7i;
2619        x0i = y5i + y7r;
2620        x2r = wn4r * (x0r - x0i);
2621        x2i = wn4r * (x0i + x0r);
2622        a[offa + 8] = x1r + x2r;
2623        a[offa + 9] = x1i + x2i;
2624        a[offa + 10] = x1r - x2r;
2625        a[offa + 11] = x1i - x2i;
2626        x1r = y4r + y6i;
2627        x1i = y4i - y6r;
2628        x0r = y5r + y7i;
2629        x0i = y5i - y7r;
2630        x2r = wn4r * (x0r - x0i);
2631        x2i = wn4r * (x0i + x0r);
2632        a[offa + 12] = x1r - x2i;
2633        a[offa + 13] = x1i + x2r;
2634        a[offa + 14] = x1r + x2i;
2635        a[offa + 15] = x1i - x2r;
2636        x1r = y8r + y10r;
2637        x1i = y8i + y10i;
2638        x2r = y9r - y11r;
2639        x2i = y9i - y11i;
2640        a[offa + 16] = x1r + x2r;
2641        a[offa + 17] = x1i + x2i;
2642        a[offa + 18] = x1r - x2r;
2643        a[offa + 19] = x1i - x2i;
2644        x1r = y8r - y10r;
2645        x1i = y8i - y10i;
2646        x2r = y9r + y11r;
2647        x2i = y9i + y11i;
2648        a[offa + 20] = x1r - x2i;
2649        a[offa + 21] = x1i + x2r;
2650        a[offa + 22] = x1r + x2i;
2651        a[offa + 23] = x1i - x2r;
2652        x1r = y12r - y14i;
2653        x1i = y12i + y14r;
2654        x0r = y13r + y15i;
2655        x0i = y13i - y15r;
2656        x2r = wn4r * (x0r - x0i);
2657        x2i = wn4r * (x0i + x0r);
2658        a[offa + 24] = x1r + x2r;
2659        a[offa + 25] = x1i + x2i;
2660        a[offa + 26] = x1r - x2r;
2661        a[offa + 27] = x1i - x2i;
2662        x1r = y12r + y14i;
2663        x1i = y12i - y14r;
2664        x0r = y13r - y15i;
2665        x0i = y13i + y15r;
2666        x2r = wn4r * (x0r - x0i);
2667        x2i = wn4r * (x0i + x0r);
2668        a[offa + 28] = x1r - x2i;
2669        a[offa + 29] = x1i + x2r;
2670        a[offa + 30] = x1r + x2i;
2671        a[offa + 31] = x1i - x2r;
2672    }
2673
2674    private void cftf081(double[] a, int offa, double[] w, int startw) {
2675        double wn4r, x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i, y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i, y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i;
2676
2677        wn4r = w[startw + 1];
2678        x0r = a[offa] + a[offa + 8];
2679        x0i = a[offa + 1] + a[offa + 9];
2680        x1r = a[offa] - a[offa + 8];
2681        x1i = a[offa + 1] - a[offa + 9];
2682        x2r = a[offa + 4] + a[offa + 12];
2683        x2i = a[offa + 5] + a[offa + 13];
2684        x3r = a[offa + 4] - a[offa + 12];
2685        x3i = a[offa + 5] - a[offa + 13];
2686        y0r = x0r + x2r;
2687        y0i = x0i + x2i;
2688        y2r = x0r - x2r;
2689        y2i = x0i - x2i;
2690        y1r = x1r - x3i;
2691        y1i = x1i + x3r;
2692        y3r = x1r + x3i;
2693        y3i = x1i - x3r;
2694        x0r = a[offa + 2] + a[offa + 10];
2695        x0i = a[offa + 3] + a[offa + 11];
2696        x1r = a[offa + 2] - a[offa + 10];
2697        x1i = a[offa + 3] - a[offa + 11];
2698        x2r = a[offa + 6] + a[offa + 14];
2699        x2i = a[offa + 7] + a[offa + 15];
2700        x3r = a[offa + 6] - a[offa + 14];
2701        x3i = a[offa + 7] - a[offa + 15];
2702        y4r = x0r + x2r;
2703        y4i = x0i + x2i;
2704        y6r = x0r - x2r;
2705        y6i = x0i - x2i;
2706        x0r = x1r - x3i;
2707        x0i = x1i + x3r;
2708        x2r = x1r + x3i;
2709        x2i = x1i - x3r;
2710        y5r = wn4r * (x0r - x0i);
2711        y5i = wn4r * (x0r + x0i);
2712        y7r = wn4r * (x2r - x2i);
2713        y7i = wn4r * (x2r + x2i);
2714        a[offa + 8] = y1r + y5r;
2715        a[offa + 9] = y1i + y5i;
2716        a[offa + 10] = y1r - y5r;
2717        a[offa + 11] = y1i - y5i;
2718        a[offa + 12] = y3r - y7i;
2719        a[offa + 13] = y3i + y7r;
2720        a[offa + 14] = y3r + y7i;
2721        a[offa + 15] = y3i - y7r;
2722        a[offa] = y0r + y4r;
2723        a[offa + 1] = y0i + y4i;
2724        a[offa + 2] = y0r - y4r;
2725        a[offa + 3] = y0i - y4i;
2726        a[offa + 4] = y2r - y6i;
2727        a[offa + 5] = y2i + y6r;
2728        a[offa + 6] = y2r + y6i;
2729        a[offa + 7] = y2i - y6r;
2730    }
2731
2732    private void cftf082(double[] a, int offa, double[] w, int startw) {
2733        double wn4r, wk1r, wk1i, x0r, x0i, x1r, x1i, y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i, y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i;
2734
2735        wn4r = w[startw + 1];
2736        wk1r = w[startw + 2];
2737        wk1i = w[startw + 3];
2738        y0r = a[offa] - a[offa + 9];
2739        y0i = a[offa + 1] + a[offa + 8];
2740        y1r = a[offa] + a[offa + 9];
2741        y1i = a[offa + 1] - a[offa + 8];
2742        x0r = a[offa + 4] - a[offa + 13];
2743        x0i = a[offa + 5] + a[offa + 12];
2744        y2r = wn4r * (x0r - x0i);
2745        y2i = wn4r * (x0i + x0r);
2746        x0r = a[offa + 4] + a[offa + 13];
2747        x0i = a[offa + 5] - a[offa + 12];
2748        y3r = wn4r * (x0r - x0i);
2749        y3i = wn4r * (x0i + x0r);
2750        x0r = a[offa + 2] - a[offa + 11];
2751        x0i = a[offa + 3] + a[offa + 10];
2752        y4r = wk1r * x0r - wk1i * x0i;
2753        y4i = wk1r * x0i + wk1i * x0r;
2754        x0r = a[offa + 2] + a[offa + 11];
2755        x0i = a[offa + 3] - a[offa + 10];
2756        y5r = wk1i * x0r - wk1r * x0i;
2757        y5i = wk1i * x0i + wk1r * x0r;
2758        x0r = a[offa + 6] - a[offa + 15];
2759        x0i = a[offa + 7] + a[offa + 14];
2760        y6r = wk1i * x0r - wk1r * x0i;
2761        y6i = wk1i * x0i + wk1r * x0r;
2762        x0r = a[offa + 6] + a[offa + 15];
2763        x0i = a[offa + 7] - a[offa + 14];
2764        y7r = wk1r * x0r - wk1i * x0i;
2765        y7i = wk1r * x0i + wk1i * x0r;
2766        x0r = y0r + y2r;
2767        x0i = y0i + y2i;
2768        x1r = y4r + y6r;
2769        x1i = y4i + y6i;
2770        a[offa] = x0r + x1r;
2771        a[offa + 1] = x0i + x1i;
2772        a[offa + 2] = x0r - x1r;
2773        a[offa + 3] = x0i - x1i;
2774        x0r = y0r - y2r;
2775        x0i = y0i - y2i;
2776        x1r = y4r - y6r;
2777        x1i = y4i - y6i;
2778        a[offa + 4] = x0r - x1i;
2779        a[offa + 5] = x0i + x1r;
2780        a[offa + 6] = x0r + x1i;
2781        a[offa + 7] = x0i - x1r;
2782        x0r = y1r - y3i;
2783        x0i = y1i + y3r;
2784        x1r = y5r - y7r;
2785        x1i = y5i - y7i;
2786        a[offa + 8] = x0r + x1r;
2787        a[offa + 9] = x0i + x1i;
2788        a[offa + 10] = x0r - x1r;
2789        a[offa + 11] = x0i - x1i;
2790        x0r = y1r + y3i;
2791        x0i = y1i - y3r;
2792        x1r = y5r + y7r;
2793        x1i = y5i + y7i;
2794        a[offa + 12] = x0r - x1i;
2795        a[offa + 13] = x0i + x1r;
2796        a[offa + 14] = x0r + x1i;
2797        a[offa + 15] = x0i - x1r;
2798    }
2799
2800    private void cftf040(double[] a, int offa) {
2801        double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
2802
2803        x0r = a[offa] + a[offa + 4];
2804        x0i = a[offa + 1] + a[offa + 5];
2805        x1r = a[offa] - a[offa + 4];
2806        x1i = a[offa + 1] - a[offa + 5];
2807        x2r = a[offa + 2] + a[offa + 6];
2808        x2i = a[offa + 3] + a[offa + 7];
2809        x3r = a[offa + 2] - a[offa + 6];
2810        x3i = a[offa + 3] - a[offa + 7];
2811        a[offa] = x0r + x2r;
2812        a[offa + 1] = x0i + x2i;
2813        a[offa + 2] = x1r - x3i;
2814        a[offa + 3] = x1i + x3r;
2815        a[offa + 4] = x0r - x2r;
2816        a[offa + 5] = x0i - x2i;
2817        a[offa + 6] = x1r + x3i;
2818        a[offa + 7] = x1i - x3r;
2819    }
2820
2821    private void cftb040(double[] a, int offa) {
2822        double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
2823
2824        x0r = a[offa] + a[offa + 4];
2825        x0i = a[offa + 1] + a[offa + 5];
2826        x1r = a[offa] - a[offa + 4];
2827        x1i = a[offa + 1] - a[offa + 5];
2828        x2r = a[offa + 2] + a[offa + 6];
2829        x2i = a[offa + 3] + a[offa + 7];
2830        x3r = a[offa + 2] - a[offa + 6];
2831        x3i = a[offa + 3] - a[offa + 7];
2832        a[offa] = x0r + x2r;
2833        a[offa + 1] = x0i + x2i;
2834        a[offa + 2] = x1r + x3i;
2835        a[offa + 3] = x1i - x3r;
2836        a[offa + 4] = x0r - x2r;
2837        a[offa + 5] = x0i - x2i;
2838        a[offa + 6] = x1r - x3i;
2839        a[offa + 7] = x1i + x3r;
2840    }
2841
2842    private void cftx020(double[] a, int offa) {
2843        double x0r, x0i;
2844
2845        x0r = a[offa] - a[offa + 2];
2846        x0i = a[offa + 1] - a[offa + 3];
2847        a[offa] += a[offa + 2];
2848        a[offa + 1] += a[offa + 3];
2849        a[offa + 2] = x0r;
2850        a[offa + 3] = x0i;
2851    }
2852
2853    private void rftfsub(int n, double[] a, int offa, int nc, double[] c, int startc) {
2854        int k, kk, ks, m;
2855        double wkr, wki, xr, xi, yr, yi;
2856        int idx1, idx2;
2857        m = n >> 1;
2858        ks = 2 * nc / m;
2859        kk = 0;
2860        for (int j = 2; j < m; j += 2) {
2861            k = n - j;
2862            kk += ks;
2863            wkr = 0.5 - c[startc + nc - kk];
2864            wki = c[startc + kk];
2865            idx1 = offa + j;
2866            idx2 = offa + k;
2867            xr = a[idx1] - a[idx2];
2868            xi = a[idx1 + 1] + a[idx2 + 1];
2869            yr = wkr * xr - wki * xi;
2870            yi = wkr * xi + wki * xr;
2871            a[idx1] -= yr;
2872            a[idx1 + 1] -= yi;
2873            a[idx2] += yr;
2874            a[idx2 + 1] -= yi;
2875        }
2876    }
2877
2878    private void rftbsub(int n, double[] a, int offa, int nc, double[] c, int startc) {
2879        int k, kk, ks, m;
2880        double wkr, wki, xr, xi, yr, yi;
2881        int idx1, idx2;
2882        m = n >> 1;
2883        ks = 2 * nc / m;
2884        kk = 0;
2885        for (int j = 2; j < m; j += 2) {
2886            k = n - j;
2887            kk += ks;
2888            wkr = 0.5 - c[startc + nc - kk];
2889            wki = c[startc + kk];
2890            idx1 = offa + j;
2891            idx2 = offa + k;
2892            xr = a[idx1] - a[idx2];
2893            xi = a[idx1 + 1] + a[idx2 + 1];
2894            yr = wkr * xr + wki * xi;
2895            yi = wkr * xi - wki * xr;
2896            a[idx1] -= yr;
2897            a[idx1 + 1] -= yi;
2898            a[idx2] += yr;
2899            a[idx2 + 1] -= yi;
2900        }
2901    }
2902
2903    private void dctsub(int n, double[] a, int offa, int nc, double[] c, int startc) {
2904        int k, kk, ks, m;
2905        double wkr, wki, xr;
2906        int idx0, idx1, idx2;
2907
2908        m = n >> 1;
2909        ks = nc / n;
2910        kk = 0;
2911        for (int j = 1; j < m; j++) {
2912            k = n - j;
2913            kk += ks;
2914            idx0 = startc + kk;
2915            idx1 = offa + j;
2916            idx2 = offa + k;
2917            wkr = c[idx0] - c[startc + nc - kk];
2918            wki = c[idx0] + c[startc + nc - kk];
2919            xr = wki * a[idx1] - wkr * a[idx2];
2920            a[idx1] = wkr * a[idx1] + wki * a[idx2];
2921            a[idx2] = xr;
2922        }
2923        a[offa + m] *= c[startc];
2924    }
2925
2926    private void scale(final double m, final double[] a, final int offa) {
2927        int nthreads = ConcurrencyUtils.getNumberOfThreads();
2928        if ((nthreads > 1) && (n > ConcurrencyUtils.getThreadsBeginN_1D_FFT_2Threads())) {
2929            nthreads = 2;
2930            final int k = n / nthreads;
2931            Future<?>[] futures = new Future[nthreads];
2932            for (int i = 0; i < nthreads; i++) {
2933                final int firstIdx = offa + i * k;
2934                final int lastIdx = (i == (nthreads - 1)) ? offa + n : firstIdx + k;
2935                futures[i] = ConcurrencyUtils.submit(new Runnable() {
2936                    @Override
2937                                        public void run() {
2938                        for (int i = firstIdx; i < lastIdx; i++) {
2939                            a[i] *= m;
2940                        }
2941
2942                    }
2943                });
2944            }
2945            ConcurrencyUtils.waitForCompletion(futures);
2946        } else {
2947            int firstIdx = offa;
2948            int lastIdx = offa + n;
2949            for (int i = firstIdx; i < lastIdx; i++) {
2950                a[i] *= m;
2951            }
2952        }
2953    }
2954}