/* vim:syntax=c
 * vim:sw=4
 *
 * Interfaces to the DCT transforms of fftpack
 */
#include <math.h>

#include "fftpack.h"

enum normalize {
    DCT_NORMALIZE_NO = 0,
    DCT_NORMALIZE_ORTHONORMAL = 1
};

/**begin repeat

#type=float,double#
#pref=,d#
#PREF=,D#
*/
extern void F_FUNC(@pref@costi, @PREF@COSTI)(int*, @type@*);
extern void F_FUNC(@pref@cost, @PREF@COST)(int*, @type@*, @type@*);
extern void F_FUNC(@pref@cosqi, @PREF@COSQI)(int*, @type@*);
extern void F_FUNC(@pref@cosqb, @PREF@COSQB)(int*, @type@*, @type@*);
extern void F_FUNC(@pref@cosqf, @PREF@COSQF)(int*, @type@*, @type@*);

GEN_CACHE(@pref@dct1,(int n)
      ,@type@* wsave;
      ,(caches_@pref@dct1[i].n==n)
      ,caches_@pref@dct1[id].wsave = malloc(sizeof(@type@)*(3*n+15));
       F_FUNC(@pref@costi, @PREF@COSTI)(&n, caches_@pref@dct1[id].wsave);
      ,free(caches_@pref@dct1[id].wsave);
      ,10)

GEN_CACHE(@pref@dct2,(int n)
      ,@type@* wsave;
      ,(caches_@pref@dct2[i].n==n)
      ,caches_@pref@dct2[id].wsave = malloc(sizeof(@type@)*(3*n+15));
       F_FUNC(@pref@cosqi,@PREF@COSQI)(&n,caches_@pref@dct2[id].wsave);
      ,free(caches_@pref@dct2[id].wsave);
      ,10)

void @pref@dct1(@type@ * inout, int n, int howmany, int normalize)
{
    int i, j;
    @type@ *ptr = inout, n1, n2;
    @type@ *wsave = NULL;

    wsave = caches_@pref@dct1[get_cache_id_@pref@dct1(n)].wsave;

    for (i = 0; i < howmany; ++i, ptr += n) {
        F_FUNC(@pref@cost, @PREF@COST)(&n, ptr, wsave);
    }

    switch (normalize) {
        case DCT_NORMALIZE_NO:
            break;
#if 0
        case DCT_NORMALIZE_ORTHONORMAL:
            ptr = inout;
            n1 = sqrt(0.5 / (n-1));
            n2 = sqrt(1. / (n-1));
            for (i = 0; i < howmany; ++i, ptr+=n) {
                ptr[0] *= n1;
                for (j = 1; j < n-1; ++j) {
                    ptr[j] *= n2;
                }
            }
            break;
#endif
        default:
            fprintf(stderr, "dct1: normalize not yet supported=%d\n",
                    normalize);
            break;
    }
}

void @pref@dct2(@type@ * inout, int n, int howmany, int normalize)
{
    int i, j;
    @type@ *ptr = inout;
    @type@ *wsave = NULL;
    @type@ n1, n2;

    wsave = caches_@pref@dct2[get_cache_id_@pref@dct2(n)].wsave;

    for (i = 0; i < howmany; ++i, ptr += n) {
        F_FUNC(@pref@cosqb, @PREF@COSQB)(&n, ptr, wsave);

    }

    switch (normalize) {
        case DCT_NORMALIZE_NO:
            ptr = inout;
            /* 0.5 coeff comes from fftpack defining DCT as
             * 4 * sum(cos(something)), whereas most definition
             * use 2 */
            for (i = 0; i < n * howmany; ++i) {
                ptr[i] *= 0.5;
            }
            break;
        case DCT_NORMALIZE_ORTHONORMAL:
            ptr = inout;
            /* 0.5 coeff comes from fftpack defining DCT as
             * 4 * sum(cos(something)), whereas most definition
             * use 2 */
            n1 = 0.25 * sqrt(1./n);
            n2 = 0.25 * sqrt(2./n);
            for (i = 0; i < howmany; ++i, ptr+=n) {
                ptr[0] *= n1;
                for (j = 1; j < n; ++j) {
                    ptr[j] *= n2;
                }
            }
            break;
        default:
            fprintf(stderr, "dct2: normalize not yet supported=%d\n",
                    normalize);
            break;
    }
}

void @pref@dct3(@type@ * inout, int n, int howmany, int normalize)
{
    int i, j;
    @type@ *ptr = inout;
    @type@ *wsave = NULL;
    @type@ n1, n2;

    wsave = caches_@pref@dct2[get_cache_id_@pref@dct2(n)].wsave;

    switch (normalize) {
        case DCT_NORMALIZE_NO:
            break;
        case DCT_NORMALIZE_ORTHONORMAL:
            n1 = sqrt(1./n);
            n2 = sqrt(0.5/n);
            for (i = 0; i < howmany; ++i, ptr+=n) {
                ptr[0] *= n1;
                for (j = 1; j < n; ++j) {
                    ptr[j] *= n2;
                }
            }
            break;
        default:
            fprintf(stderr, "dct3: normalize not yet supported=%d\n",
                    normalize);
            break;
    }

    ptr = inout;
    for (i = 0; i < howmany; ++i, ptr += n) {
        F_FUNC(@pref@cosqf, @PREF@COSQF)(&n, ptr, wsave);

    }

}
/**end repeat**/
