#line 1 "scipy/fftpack/src/dst.c.src"

/*
 *****************************************************************************
 **       This file was autogenerated from a template  DO NOT EDIT!!!!      **
 **       Changes should be made to the original source (.src) file         **
 *****************************************************************************
 */

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

#include "fftpack.h"

enum normalize {
    DST_NORMALIZE_NO = 0,
    DST_NORMALIZE_ORTHONORMAL = 1
};

#line 21
extern void F_FUNC(sinti, SINTI)(int*, float*);
extern void F_FUNC(sint, SINT)(int*, float*, float*);
extern void F_FUNC(sinqi, SINQI)(int*, float*);
extern void F_FUNC(sinqb, SINQB)(int*, float*, float*);
extern void F_FUNC(sinqf, SINQF)(int*, float*, float*);
extern void dct4(float * inout, int n, int howmany, int normalize);

GEN_CACHE(dst1,(int n)
      ,float* wsave;
      ,(caches_dst1[i].n==n)
      ,caches_dst1[id].wsave = malloc(sizeof(float)*(3*n+15));
       F_FUNC(sinti, SINTI)(&n, caches_dst1[id].wsave);
      ,free(caches_dst1[id].wsave);
      ,10)

GEN_CACHE(dst2,(int n)
      ,float* wsave;
      ,(caches_dst2[i].n==n)
      ,caches_dst2[id].wsave = malloc(sizeof(float)*(3*n+15));
       F_FUNC(sinqi,SINQI)(&n,caches_dst2[id].wsave);
      ,free(caches_dst2[id].wsave);
      ,10)


void dst1(float * inout, int n, int howmany, int normalize)
{
    int i;
    float *ptr = inout;
    float *wsave = NULL;

    wsave = caches_dst1[get_cache_id_dst1(n)].wsave;

    for (i = 0; i < howmany; ++i, ptr += n) {
        F_FUNC(sint, SINT)(&n, ptr, wsave);
    }

    switch (normalize) {
        case DST_NORMALIZE_NO:
            break;
        case DST_NORMALIZE_ORTHONORMAL:
            {
                float m = 1.0 / sqrt(2 * (n+1));
                int j;
                ptr = inout;
                for (i = 0; i < howmany; ++i, ptr += n) {
                    for (j=0; j<n; j++) {
                        ptr[j] *= m;
                    }
                }
            }
            break;
        default:
            fprintf(stderr, "dst1: normalize not yet supported=%d\n",
                    normalize);
            break;
    }
}

void dst2(float * inout, int n, int howmany, int normalize)
{
    int i, j;
    float *ptr = inout;
    float *wsave = NULL;
    float n1, n2;

    wsave = caches_dst2[get_cache_id_dst2(n)].wsave;

    for (i = 0; i < howmany; ++i, ptr += n) {
        F_FUNC(sinqb, SINQB)(&n, ptr, wsave);

    }

    switch (normalize) {
        case DST_NORMALIZE_NO:
            ptr = inout;
            /* 0.5 coeff comes from fftpack defining DST as
             * 4 * sum(sin(something)), whereas most definition
             * use 2 */
            for (i = 0; i < n * howmany; ++i) {
                ptr[i] *= 0.5;
            }
            break;
        case DST_NORMALIZE_ORTHONORMAL:
            ptr = inout;
            /* 0.5 coeff comes from fftpack defining DST as
             * 4 * sum(sin(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, "dst2: normalize not yet supported=%d\n",
                    normalize);
            break;
    }
}

void dst3(float * inout, int n, int howmany, int normalize)
{
    int i, j;
    float *ptr = inout;
    float *wsave = NULL;
    float n1, n2;

    wsave = caches_dst2[get_cache_id_dst2(n)].wsave;

    switch (normalize) {
        case DST_NORMALIZE_NO:
            break;
        case DST_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, "dst3: normalize not yet supported=%d\n",
                    normalize);
            break;
    }

    ptr = inout;
    for (i = 0; i < howmany; ++i, ptr += n) {
        F_FUNC(sinqf, SINQF)(&n, ptr, wsave);

    }

}

void dst4(float * inout, int n, int howmany, int normalize)
{
    int i, j, n2=n/2;
    float *ptr = inout;
    float tmp;

    /* reverse order of inputs */
    for (i = 0; i < howmany; ++i, ptr += n) {
        for (j=0; j<n2; j++) {
            tmp = ptr[j];
            ptr[j] = ptr[n-j-1];
            ptr[n-j-1] = tmp;
        }
    }

    /* apply DCT-IV, we can use it's normalization */
    dct4(inout, n, howmany, normalize);

    /* flip sign of odd indices */
    ptr = inout;
    for (i = 0; i < howmany; ++i, ptr += n) {
        for (j=1; j<n; j+=2) {
            ptr[j] = -ptr[j];
        }
    }
}


#line 21
extern void F_FUNC(dsinti, DSINTI)(int*, double*);
extern void F_FUNC(dsint, DSINT)(int*, double*, double*);
extern void F_FUNC(dsinqi, DSINQI)(int*, double*);
extern void F_FUNC(dsinqb, DSINQB)(int*, double*, double*);
extern void F_FUNC(dsinqf, DSINQF)(int*, double*, double*);
extern void ddct4(double * inout, int n, int howmany, int normalize);

GEN_CACHE(ddst1,(int n)
      ,double* wsave;
      ,(caches_ddst1[i].n==n)
      ,caches_ddst1[id].wsave = malloc(sizeof(double)*(3*n+15));
       F_FUNC(dsinti, DSINTI)(&n, caches_ddst1[id].wsave);
      ,free(caches_ddst1[id].wsave);
      ,10)

GEN_CACHE(ddst2,(int n)
      ,double* wsave;
      ,(caches_ddst2[i].n==n)
      ,caches_ddst2[id].wsave = malloc(sizeof(double)*(3*n+15));
       F_FUNC(dsinqi,DSINQI)(&n,caches_ddst2[id].wsave);
      ,free(caches_ddst2[id].wsave);
      ,10)


void ddst1(double * inout, int n, int howmany, int normalize)
{
    int i;
    double *ptr = inout;
    double *wsave = NULL;

    wsave = caches_ddst1[get_cache_id_ddst1(n)].wsave;

    for (i = 0; i < howmany; ++i, ptr += n) {
        F_FUNC(dsint, DSINT)(&n, ptr, wsave);
    }

    switch (normalize) {
        case DST_NORMALIZE_NO:
            break;
        case DST_NORMALIZE_ORTHONORMAL:
            {
                double m = 1.0 / sqrt(2 * (n+1));
                int j;
                ptr = inout;
                for (i = 0; i < howmany; ++i, ptr += n) {
                    for (j=0; j<n; j++) {
                        ptr[j] *= m;
                    }
                }
            }
            break;
        default:
            fprintf(stderr, "dst1: normalize not yet supported=%d\n",
                    normalize);
            break;
    }
}

void ddst2(double * inout, int n, int howmany, int normalize)
{
    int i, j;
    double *ptr = inout;
    double *wsave = NULL;
    double n1, n2;

    wsave = caches_ddst2[get_cache_id_ddst2(n)].wsave;

    for (i = 0; i < howmany; ++i, ptr += n) {
        F_FUNC(dsinqb, DSINQB)(&n, ptr, wsave);

    }

    switch (normalize) {
        case DST_NORMALIZE_NO:
            ptr = inout;
            /* 0.5 coeff comes from fftpack defining DST as
             * 4 * sum(sin(something)), whereas most definition
             * use 2 */
            for (i = 0; i < n * howmany; ++i) {
                ptr[i] *= 0.5;
            }
            break;
        case DST_NORMALIZE_ORTHONORMAL:
            ptr = inout;
            /* 0.5 coeff comes from fftpack defining DST as
             * 4 * sum(sin(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, "dst2: normalize not yet supported=%d\n",
                    normalize);
            break;
    }
}

void ddst3(double * inout, int n, int howmany, int normalize)
{
    int i, j;
    double *ptr = inout;
    double *wsave = NULL;
    double n1, n2;

    wsave = caches_ddst2[get_cache_id_ddst2(n)].wsave;

    switch (normalize) {
        case DST_NORMALIZE_NO:
            break;
        case DST_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, "dst3: normalize not yet supported=%d\n",
                    normalize);
            break;
    }

    ptr = inout;
    for (i = 0; i < howmany; ++i, ptr += n) {
        F_FUNC(dsinqf, DSINQF)(&n, ptr, wsave);

    }

}

void ddst4(double * inout, int n, int howmany, int normalize)
{
    int i, j, n2=n/2;
    double *ptr = inout;
    double tmp;

    /* reverse order of inputs */
    for (i = 0; i < howmany; ++i, ptr += n) {
        for (j=0; j<n2; j++) {
            tmp = ptr[j];
            ptr[j] = ptr[n-j-1];
            ptr[n-j-1] = tmp;
        }
    }

    /* apply DCT-IV, we can use it's normalization */
    ddct4(inout, n, howmany, normalize);

    /* flip sign of odd indices */
    ptr = inout;
    for (i = 0; i < howmany; ++i, ptr += n) {
        for (j=1; j<n; j+=2) {
            ptr[j] = -ptr[j];
        }
    }
}



