/*
 * Copyright (C) 2000 ymnk<ymnk@jcraft.com>
 *               2015 Trilarion
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License as published
 * by the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */
package com.github.trilarion.sound.vorbis.jcraft.jorbis;

import java.util.logging.Logger;

/**
 *
 */
class Mdct {

    private static final Logger LOG = Logger.getLogger(Mdct.class.getName());

    int n;
    int log2n;

    float[] trig;
    int[] bitrev;

    float scale;

    void init(int n) {
        bitrev = new int[n / 4];
        trig = new float[n + n / 4];

        log2n = (int) Math.rint(Math.log(n) / Math.log(2));
        this.n = n;

        int AE = 0;
        int AO = 1;
        int BE = AE + n / 2;
        int BO = BE + 1;
        int CE = BE + n / 2;
        int CO = CE + 1;
        // trig lookups...
        for (int i = 0; i < n / 4; i++) {
            trig[AE + i * 2] = (float) Math.cos((Math.PI / n) * (4 * i));
            trig[AO + i * 2] = (float) -Math.sin((Math.PI / n) * (4 * i));
            trig[BE + i * 2] = (float) Math.cos((Math.PI / (2 * n)) * (2 * i + 1));
            trig[BO + i * 2] = (float) Math.sin((Math.PI / (2 * n)) * (2 * i + 1));
        }
        for (int i = 0; i < n / 8; i++) {
            trig[CE + i * 2] = (float) Math.cos((Math.PI / n) * (4 * i + 2));
            trig[CO + i * 2] = (float) -Math.sin((Math.PI / n) * (4 * i + 2));
        }

        {
            int mask = (1 << (log2n - 1)) - 1;
            int msb = 1 << (log2n - 2);
            for (int i = 0; i < n / 8; i++) {
                int acc = 0;
                for (int j = 0; msb >>> j != 0; j++) {
                    if (((msb >>> j) & i) != 0) {
                        acc |= 1 << j;
                    }
                }
                bitrev[i * 2] = ((~acc) & mask);
                //	bitrev[i*2]=((~acc)&mask)-1;
                bitrev[i * 2 + 1] = acc;
            }
        }
        scale = 4.f / n;
    }

    void clear() {
    }

    void forward(float[] in, float[] out) {
    }

    float[] _x = new float[1024];
    float[] _w = new float[1024];

    synchronized void backward(float[] in, float[] out) {
        if (_x.length < n / 2) {
            _x = new float[n / 2];
        }
        if (_w.length < n / 2) {
            _w = new float[n / 2];
        }
        float[] x = _x;
        float[] w = _w;
        int n2 = n >>> 1;
        int n4 = n >>> 2;
        int n8 = n >>> 3;

        // rotate + step 1
        {
            int inO = 1;
            int xO = 0;
            int A = n2;

            int i;
            for (i = 0; i < n8; i++) {
                A -= 2;
                x[xO++] = -in[inO + 2] * trig[A + 1] - in[inO] * trig[A];
                x[xO++] = in[inO] * trig[A + 1] - in[inO + 2] * trig[A];
                inO += 4;
            }

            inO = n2 - 4;

            for (i = 0; i < n8; i++) {
                A -= 2;
                x[xO++] = in[inO] * trig[A + 1] + in[inO + 2] * trig[A];
                x[xO++] = in[inO] * trig[A] - in[inO + 2] * trig[A + 1];
                inO -= 4;
            }
        }

        float[] xxx = mdct_kernel(x, w, n, n2, n4, n8);
        int xx = 0;

        // step 8
        {
            int B = n2;
            int o1 = n4, o2 = o1 - 1;
            int o3 = n4 + n2, o4 = o3 - 1;

            for (int i = 0; i < n4; i++) {
                float temp1 = (xxx[xx] * trig[B + 1] - xxx[xx + 1] * trig[B]);
                float temp2 = -(xxx[xx] * trig[B] + xxx[xx + 1] * trig[B + 1]);

                out[o1] = -temp1;
                out[o2] = temp1;
                out[o3] = temp2;
                out[o4] = temp2;

                o1++;
                o2--;
                o3++;
                o4--;
                xx += 2;
                B += 2;
            }
        }
    }

    private float[] mdct_kernel(float[] x, float[] w, int n, int n2, int n4,
            int n8) {
        // step 2

        int xA = n4;
        int xB = 0;
        int w2 = n4;
        int A = n2;

        for (int i = 0; i < n4;) {
            float x0 = x[xA] - x[xB];
            float x1;
            w[w2 + i] = x[xA++] + x[xB++];

            x1 = x[xA] - x[xB];
            A -= 4;

            w[i++] = x0 * trig[A] + x1 * trig[A + 1];
            w[i] = x1 * trig[A] - x0 * trig[A + 1];

            w[w2 + i] = x[xA++] + x[xB++];
            i++;
        }

        // step 3
        {
            for (int i = 0; i < log2n - 3; i++) {
                int k0 = n >>> (i + 2);
                int k1 = 1 << (i + 3);
                int wbase = n2 - 2;

                A = 0;
                float[] temp;

                for (int r = 0; r < (k0 >>> 2); r++) {
                    int w1 = wbase;
                    w2 = w1 - (k0 >> 1);
                    float AEv = trig[A], wA;
                    float AOv = trig[A + 1], wB;
                    wbase -= 2;

                    k0++;
                    for (int s = 0; s < (2 << i); s++) {
                        wB = w[w1] - w[w2];
                        x[w1] = w[w1] + w[w2];

                        wA = w[++w1] - w[++w2];
                        x[w1] = w[w1] + w[w2];

                        x[w2] = wA * AEv - wB * AOv;
                        x[w2 - 1] = wB * AEv + wA * AOv;

                        w1 -= k0;
                        w2 -= k0;
                    }
                    k0--;
                    A += k1;
                }

                temp = w;
                w = x;
                x = temp;
            }
        }

        // step 4, 5, 6, 7
        {
            int C = n;
            int bit = 0;
            int x1 = 0;
            int x2 = n2 - 1;

            for (int i = 0; i < n8; i++) {
                int t1 = bitrev[bit++];
                int t2 = bitrev[bit++];

                float wA = w[t1] - w[t2 + 1];
                float wB = w[t1 - 1] + w[t2];
                float wC = w[t1] + w[t2 + 1];
                float wD = w[t1 - 1] - w[t2];

                float wACE = wA * trig[C];
                float wBCE = wB * trig[C++];
                float wACO = wA * trig[C];
                float wBCO = wB * trig[C++];

                x[x1++] = (wC + wACO + wBCE) * .5f;
                x[x2--] = (-wD + wBCO - wACE) * .5f;
                x[x1++] = (wD + wBCO - wACE) * .5f;
                x[x2--] = (wC - wACO - wBCE) * .5f;
            }
        }
        return (x);
    }
}
