void vector_hartley_transform_simple(
    double * v, long size)                     /* size must be a power of 2 */
{
    long old_len = 0;                        /* initialize to avoid warning */
    long len;
    long new_len = 0;                        /* initialize to avoid warning */

    vector_bitreverse_indices(v,size); /* v[b0 b1 .. bn] <-> v[bn .. b1 b0] */

    for (len = 1; len < size; old_len = len, len = new_len)
    {
        long i, j;
        double hi, hj;
        new_len = 2 * len;              /* build transform of double length */

        for (i = 0; i < size; i += new_len)             /* for all blocks */
        {
            j = i + len;                                /* special case: */
            hi = v[i]; hj = v[j];
            v[i] = hi + hj;                             /* f = 0 */
            v[j] = hi - hj;                             /* f = PI */
        }
        if (len < 2) continue;
        for (i = old_len; i < size; i += new_len)       /* for all blocks */
        {
            j = i + len;                                /* special case: */
            hi = v[i]; hj = v[j];
            v[i] = hi + hj;                             /* f = PI/2 */
            v[j] = hi - hj;                             /* f = 3 * PI/2 */
        }
        if (len >= 4)
        {
            double d = MATH_2_MUL_PI / new_len;
            double a = 2.0 * M_SQR(sin(d * 0.5));       /* initialize trig. */
            double b = sin(d);                          /* recurrence */
            double cos_t = 1.0;
            double sin_t = 0.0;
            long f;
            for (f = 1; f < old_len; f++)               /* for all freqs in */
            {                                           /*   the first quad */
                double one = a * cos_t + b * sin_t;     /* trig. recurrence */
                double two = a * sin_t - b * cos_t;     
                cos_t -= one;                           /*   cos (t + d) */
                sin_t -= two;                           /*   sin (t + d) */
                for (i = f, j = len - f;                /* for all blocks */
                     i < size;
                     i += new_len, j += new_len)
                {
                    long k = i + len;
                    long l = j + len;
                    one = cos_t * v[k] + sin_t * v[l];
                    two = cos_t * v[l] - sin_t * v[k];
                    hi = v[i]; hj = v[j];
                    v[i] = hi + one; v[k] = hi - one;   /* all four quads */
                    v[j] = hj - two; v[l] = hj + two;   /* (i,j,k,l) */
                }
            }
        }
    }
}