void FastDFS(double *x, double *y, double *w, double *Wp, int N, int skip) {
    int i, Nskip;
    double *xe, *ye, *yn, *We, *Wn, re, im;
    if (N == 1) {
        for (xe = x + N*2; x < xe; x += 2, y += skip) {
            x[0] = y[0];
            x[1] = y[1];
        }
        return;
    }

    if ((N%2) != 0) {
        Nskip = N*skip;
        ye = y + Nskip, We = Wp + Nskip;
        for (Wn = Wp, i = 0; i < Nskip; i += skip, x += 2, Wn = Wp) {
            x[0] = x[1] = 0;
            for (yn = y; yn < ye; yn += skip) {
                x[0] += Wn[0]*yn[0] - Wn[1]*yn[1];
                x[1] += Wn[0]*yn[1] + Wn[1]*yn[0];
                Wn += i;
                if (Wn > We) Wn -= Nskip;
            }
        }
        return;
    }

    FastDFS(w,   y,      x,   Wp, N/2, skip*2);
    FastDFS(w+N, y+skip, x+N, Wp, N/2, skip*2);
    for (xe = x + N; x < xe; x += 2, w += 2, Wp += skip) {
        re = Wp[0]*w[N] - Wp[1]*w[N+1];
        im = Wp[0]*w[N+1] + Wp[1]*w[N];
        x[0]   = w[0] + re;
        x[1]   = w[1] + im;
        x[N]   = w[0] - re;
        x[N+1] = w[1] - im;
    }
}