#include "tensorTypes.h"
#include <time.h>
#include <stdlib.h>

// there should be a special cases that do NOT require temp space to be allocated for all primes *smaller* than DFTP_GENERIC_SIZE
#define DFTP_GENERIC_SIZE 11

#ifdef STATS
int crtRqCtr = 0;
int crtInvRqCtr = 0;
int crtCCtr = 0;
int crtInvCCtr = 0;

struct timespec crttime1 = {0,0};
struct timespec crttime2 = {0,0};
struct timespec crttime3 = {0,0};
struct timespec crttime4 = {0,0};

struct timespec crtInvRqTime = {0,0};
struct timespec crtCTime = {0,0};
struct timespec crtInvCTime = {0,0};
#endif

hDim_t bitrev (PrimeExponent pe, hDim_t j) {
    hShort_t e;
    hDim_t p = pe.prime;
    hDim_t tempj = j;
    hDim_t acc = 0;

    for(e = pe.exponent-1; e >= 0; e--) {
        div_t qr = div(tempj,p);
        acc += qr.rem * ipow(p,e);
        tempj = qr.quot;
    }
    return acc;
}

void crtTwiddleRq (hInt_t* y, hShort_t tupSize, hDim_t lts, hDim_t rts, 
                   PrimeExponent pe, hInt_t* ru, hInt_t q)
{
    hDim_t p = pe.prime;
    hShort_t e = pe.exponent;
    
#ifdef DEBUG_MODE
    ASSERT(e != 0);
#endif
    pe.exponent -= 1; // used for an argument to bitrev
    
    if(p == 2)
    {
        hDim_t mprime = 1<<(e-1);
        hDim_t blockDim = rts*mprime; // size of block in block diagonal tensor matrix

        for(hDim_t i0 = 1; i0 < mprime; i0++) // loops over i/(p-1) for i = 0..(m'-1), we can skip i0 = 0
        {
            hDim_t temp2 = i0*rts;
            hInt_t twid = ru[bitrev(pe, i0)*tupSize];

            for(hDim_t blockIdx = 0; blockIdx < lts; blockIdx++)
            {
                hDim_t temp3 = blockIdx*blockDim + temp2;
                for(hDim_t modOffset = 0; modOffset < rts; modOffset++)
                {
                    hDim_t idx = (temp3 + modOffset)*tupSize;
                    y[idx] = (y[idx]*twid) % q;
                }
            }
        }
    }
    else // This loop is faster, probably due to the division in the loop above.
    // cilk also slows it down
    {
        hDim_t mprime = ipow(p,e-1);
        hDim_t blockDim = rts*(p-1)*mprime; // size of block in block diagonal tensor matrix
        
        for(hDim_t i0 = 1; i0 < mprime; i0++) // loops over i/(p-1) for i = 0..(m'-1), we can skip i0 = 0
        {
            hDim_t temp1 = i0*(p-1);
            for(hDim_t i1 = 0; i1 < (p-1); i1++) // loops over i%(p-1) for i = 0..(m'-1)
            {        
                hDim_t temp2 = (temp1+i1)*rts;
                hInt_t twid = ru[bitrev(pe, i0)*(i1+1)*tupSize];

                for(hDim_t blockIdx = 0; blockIdx < lts; blockIdx++)
                {
                    hDim_t temp3 = blockIdx*blockDim + temp2;
                    for(hDim_t modOffset = 0; modOffset < rts; modOffset++)
                    {
                        hDim_t idx = (temp3 + modOffset)*tupSize;
                        y[idx] = (y[idx]*twid) % q;
                    }
                }
            }
        }
    }
}

// dim is power of p
void dftptwidRq (hInt_t* y, hShort_t tupSize, hDim_t lts, hDim_t rts, 
                 PrimeExponent pe, hDim_t dim, hDim_t rustride, hInt_t* ru, hInt_t q)
{
    hDim_t idx;
    hDim_t p = pe.prime;

    pe.exponent -= 1; // used for an argument to bitrev

    if(p == 2) {
        hDim_t mprime = dim>>1; // divides evenly
        hDim_t temp1 = rts*dim; // for use in computing [modified] tensorOffset
        for(hDim_t i0 = 1; i0 < mprime; i0++) // loops over i/p for i = 0..(dim-1), but we skip i0=0
        {
            hDim_t temp3 = rts*(i0*p+1);
            hInt_t twid = ru[bitrev(pe,i0)*rustride*tupSize];

            for(hDim_t blockOffset = 0; blockOffset < lts; blockOffset++)
            {
                hDim_t temp2 = blockOffset*temp1 + temp3;
                for(hDim_t modOffset = 0; modOffset < rts; modOffset++)
                {
                    idx = (temp2 + modOffset)*tupSize;
                    y[idx] = (y[idx]*twid) % q;
                }
            }
        }
    }
    else
    {
        hDim_t mprime = dim/p; // divides evenly
        hDim_t temp1 = rts*dim; // for use in computing [modified] tensorOffset
        for(hDim_t i0 = 1; i0 < mprime; i0++) // loops over i/p for i = 0..(dim-1), but we skip i0=0
        {
            for(hDim_t i1 = 1; i1 < p; i1++) // loops over i%p for i = 0..(dim-1), but we skip i1=0
            {
                hDim_t temp3 = rts*(i0*p+i1);
                hInt_t twid = ru[bitrev(pe,i0)*i1*rustride*tupSize];

                for(hDim_t blockOffset = 0; blockOffset < lts; blockOffset++)
                {
                    hDim_t temp2 = blockOffset*temp1 + temp3;
                    for(hDim_t modOffset = 0; modOffset < rts; modOffset++)
                    {
                        idx = (temp2 + modOffset)*tupSize;
                        y[idx] = (y[idx]*twid) % q;
                    }
                }
            }
        }
    }
}

//implied length of ru is rustride*p
//implied length of tempSpace is p, if p is not a special case
// temp is allowed to be NULL if p < DFTP_GENERIC_SIZE
void dftpRq (hInt_t* y, hShort_t tupSize, hDim_t lts, hDim_t rts, 
             hDim_t p, hDim_t rustride, hInt_t* ru, hInt_t* tempSpace, hInt_t q)
{
    hDim_t tensorOffset;

    if(p == 2)
    {
        hDim_t temp1 = rts<<1;

        for(hDim_t blockOffset = 0; blockOffset < lts; blockOffset++)
        {
            hDim_t temp2 = blockOffset*temp1;
            for(hDim_t modOffset = 0; modOffset < rts; modOffset++)
            {
                tensorOffset = temp2 + modOffset;
                hInt_t u = y[tensorOffset*tupSize];
                hInt_t t = y[(tensorOffset+rts)*tupSize];
                y[tensorOffset*tupSize] = (u + t) % q;
                y[(tensorOffset+rts)*tupSize] = (u - t) % q;
            }
        }
    }
    else if(p == 3)
    {
        hInt_t ru1 = ru[rustride*tupSize];
        hInt_t ru2 = ru[(rustride<<1)*tupSize];
        hDim_t temp1 = rts*3;

        for(hDim_t blockOffset = 0; blockOffset < lts; blockOffset++)
        {
            hDim_t temp2 = blockOffset*temp1;
            for(hDim_t modOffset = 0; modOffset < rts; modOffset++)
            {
                tensorOffset = temp2 + modOffset;
                hInt_t y1, y2, y3;
                y1 = y[tensorOffset*tupSize];
                y2 = y[(tensorOffset+rts)*tupSize];
                y3 = y[(tensorOffset+(rts<<1))*tupSize];
                //q is <32 bits, so we can do 3 additions without overflow
                y[tensorOffset*tupSize]          = (y1 + y2 + y3) % q;
                y[(tensorOffset+rts)*tupSize]      = (y1 + ((ru1*y2) % q) + ((ru2*y3) % q)) % q;
                y[(tensorOffset+(rts<<1))*tupSize] = (y1 + ((ru2*y2) % q) + ((ru1*y3) % q)) % q;
            }   
        }

    }
    else if(p == 5)
    {
        hDim_t temp1 = rts*5;
        hInt_t ru1 = ru[rustride*tupSize];
        hInt_t ru2 = ru[(rustride<<1)*tupSize];
        hInt_t ru3 = ru[(rustride*3)*tupSize];
        hInt_t ru4 = ru[(rustride<<2)*tupSize];

        for(hDim_t blockOffset = 0; blockOffset < lts; blockOffset++)
        {
            hDim_t temp2 = blockOffset*temp1;
            for(hDim_t modOffset = 0; modOffset < rts; modOffset++)
            {
                tensorOffset = temp2 + modOffset;
                hInt_t y1, y2, y3, y4, y5;
                y1 = y[tensorOffset*tupSize];
                y2 = y[(tensorOffset+rts)*tupSize];
                y3 = y[(tensorOffset+(rts<<1))*tupSize];
                y4 = y[(tensorOffset+3*rts)*tupSize];
                y5 = y[(tensorOffset+(rts<<2))*tupSize];
                y[tensorOffset*tupSize]          = (y1 + y2 + y3 + y4 + y5) % q;
                y[(tensorOffset+rts)*tupSize]      = (y1 + ((ru1*y2) % q) + ((ru2*y3) % q) + ((ru3*y4) % q) + ((ru4*y5) % q)) % q;
                y[(tensorOffset+(rts<<1))*tupSize] = (y1 + ((ru2*y2) % q) + ((ru4*y3) % q) + ((ru1*y4) % q) + ((ru3*y5) % q)) % q;
                y[(tensorOffset+rts*3)*tupSize]    = (y1 + ((ru3*y2) % q) + ((ru1*y3) % q) + ((ru4*y4) % q) + ((ru2*y5) % q)) % q;
                y[(tensorOffset+(rts<<2))*tupSize] = (y1 + ((ru4*y2) % q) + ((ru3*y3) % q) + ((ru2*y4) % q) + ((ru1*y5) % q)) % q;
            }
        }
    }
    else if(p == 7)
    {
        hDim_t temp1 = rts*7;
        hInt_t ru1 = ru[rustride*tupSize];
        hInt_t ru2 = ru[(rustride<<1)*tupSize];
        hInt_t ru3 = ru[(rustride*3)*tupSize];
        hInt_t ru4 = ru[(rustride<<2)*tupSize];
        hInt_t ru5 = ru[(rustride*5)*tupSize];
        hInt_t ru6 = ru[(rustride*6)*tupSize];

        for(hDim_t blockOffset = 0; blockOffset < lts; blockOffset++)
        {
            hDim_t temp2 = blockOffset*temp1;
            for(hDim_t modOffset = 0; modOffset < rts; modOffset++)
            {
                tensorOffset = temp2 + modOffset;
                hInt_t y1, y2, y3, y4, y5, y6, y7;
                y1 = y[tensorOffset*tupSize];
                y2 = y[(tensorOffset+rts)*tupSize];
                y3 = y[(tensorOffset+(rts<<1))*tupSize];
                y4 = y[(tensorOffset+3*rts)*tupSize];
                y5 = y[(tensorOffset+(rts<<2))*tupSize];
                y6 = y[(tensorOffset+rts*5)*tupSize];
                y7 = y[(tensorOffset+rts*6)*tupSize];
                y[tensorOffset*tupSize]          = (y1 +     y2 +     y3 +     y4 +     y5 +     y6 +     y7) % q;
                y[(tensorOffset+rts)*tupSize]      = (y1 + ((ru1*y2) % q) + ((ru2*y3) % q) + ((ru3*y4) % q) + ((ru4*y5) % q) + ((ru5*y6) % q) + ((ru6*y7) % q)) % q;
                y[(tensorOffset+(rts<<1))*tupSize] = (y1 + ((ru2*y2) % q) + ((ru4*y3) % q) + ((ru6*y4) % q) + ((ru1*y5) % q) + ((ru3*y6) % q) + ((ru5*y7) % q)) % q;
                y[(tensorOffset+rts*3)*tupSize]    = (y1 + ((ru3*y2) % q) + ((ru6*y3) % q) + ((ru2*y4) % q) + ((ru5*y5) % q) + ((ru1*y6) % q) + ((ru4*y7) % q)) % q;
                y[(tensorOffset+(rts<<2))*tupSize] = (y1 + ((ru4*y2) % q) + ((ru1*y3) % q) + ((ru5*y4) % q) + ((ru2*y5) % q) + ((ru6*y6) % q) + ((ru3*y7) % q)) % q;
                y[(tensorOffset+rts*5)*tupSize]    = (y1 + ((ru5*y2) % q) + ((ru3*y3) % q) + ((ru1*y4) % q) + ((ru6*y5) % q) + ((ru4*y6) % q) + ((ru2*y7) % q)) % q;
                y[(tensorOffset+rts*6)*tupSize]    = (y1 + ((ru6*y2) % q) + ((ru5*y3) % q) + ((ru4*y4) % q) + ((ru3*y5) % q) + ((ru2*y6) % q) + ((ru1*y7) % q)) % q;
            }   
        }
    }
    else
    {
        hDim_t temp1 = rts*p;
        for(hDim_t blockOffset = 0; blockOffset < lts; blockOffset++)
        {
            hDim_t temp2 = blockOffset*temp1;
            for(hDim_t modOffset = 0; modOffset < rts; modOffset++)
            {
                tensorOffset = temp2 + modOffset;                
                for(hDim_t row = 0; row < p; row++)
                {
                    hInt_t acc = 0;
                    //p is small (<< 30 bits), so we can do p additions of mod-q values without overflow
                    for(hDim_t col = 0; col < p; col++)
                    {
                        acc += ((y[(tensorOffset+col*rts)*tupSize]*ru[((col*row) % p)*rustride*tupSize])%q);
                    }
                    tempSpace[row] = acc % q;
                }
                
                for(hDim_t row = 0; row < p; row++)
                {
                    y[(tensorOffset+rts*row)*tupSize] = tempSpace[row];
                }
            }
        }
    }
}

void crtpRq (hInt_t* y, hShort_t tupSize, hDim_t lts, hDim_t rts, 
             hDim_t p, hDim_t rustride, hInt_t* ru, hInt_t q)
{
    hDim_t tensorOffset;
    if(p == 2)
    {
        return;
    }
    else if(p == 3)
    {
        hDim_t temp1 = rts*2;
        hInt_t ru1 = ru[rustride*tupSize];
        hInt_t ru2 = ru[(rustride<<1)*tupSize];

        for(hDim_t blockOffset = 0; blockOffset < lts; blockOffset++)
        {
            hDim_t temp2 = blockOffset*temp1;
            for(hDim_t modOffset = 0; modOffset < rts; modOffset++)
            {
                tensorOffset = temp2 + modOffset;
                hInt_t y1, y2;
                y1 = y[tensorOffset*tupSize];
                y2 = y[(tensorOffset+rts)*tupSize];
                y[tensorOffset*tupSize]     = (y1 + ((ru1*y2)%q)) % q;
                y[(tensorOffset+rts)*tupSize] = (y1 + ((ru2*y2)%q)) % q;
            }   
        }
    }
    else if(p == 5)
    {
        hDim_t temp1 = rts*4;
        hInt_t ru1 = ru[rustride*tupSize];
        hInt_t ru2 = ru[(rustride<<1)*tupSize];
        hInt_t ru3 = ru[(rustride*3)*tupSize];
        hInt_t ru4 = ru[(rustride<<2)*tupSize];

        for(hDim_t blockOffset = 0; blockOffset < lts; blockOffset++)
        {
            hDim_t temp2 = blockOffset*temp1;
            for(hDim_t modOffset = 0; modOffset < rts; modOffset++)
            {
                tensorOffset = temp2 + modOffset;
                hInt_t y1, y2, y3, y4;
                y1 = y[tensorOffset*tupSize];
                y2 = y[(tensorOffset+rts)*tupSize];
                y3 = y[(tensorOffset+(rts<<1))*tupSize];
                y4 = y[(tensorOffset+3*rts)*tupSize];

                y[tensorOffset*tupSize]          = (y1 + ((ru1*y2) % q) + ((ru2*y3) % q) + ((ru3*y4) % q)) % q;
                y[(tensorOffset+rts)*tupSize]      = (y1 + ((ru2*y2) % q) + ((ru4*y3) % q) + ((ru1*y4) % q)) % q;
                y[(tensorOffset+(rts<<1))*tupSize] = (y1 + ((ru3*y2) % q) + ((ru1*y3) % q) + ((ru4*y4) % q)) % q;
                y[(tensorOffset+rts*3)*tupSize]    = (y1 + ((ru4*y2) % q) + ((ru3*y3) % q) + ((ru2*y4) % q)) % q;
            }   
        }
    }
    else if(p == 7)
    {
        hDim_t temp1 = rts*6;
        hInt_t ru1 = ru[rustride*tupSize];
        hInt_t ru2 = ru[(rustride<<1)*tupSize];
        hInt_t ru3 = ru[(rustride*3)*tupSize];
        hInt_t ru4 = ru[(rustride<<2)*tupSize];
        hInt_t ru5 = ru[(rustride*5)*tupSize];
        hInt_t ru6 = ru[(rustride*6)*tupSize];
        for(hDim_t blockOffset = 0; blockOffset < lts; blockOffset++)
        {
            hDim_t temp2 = blockOffset*temp1;
            for(hDim_t modOffset = 0; modOffset < rts; modOffset++)
            {
                tensorOffset = temp2 + modOffset;
                hInt_t y1, y2, y3, y4, y5, y6;
                y1 = y[tensorOffset*tupSize];
                y2 = y[(tensorOffset+rts)*tupSize];
                y3 = y[(tensorOffset+(rts<<1))*tupSize];
                y4 = y[(tensorOffset+3*rts)*tupSize];
                y5 = y[(tensorOffset+(rts<<2))*tupSize];
                y6 = y[(tensorOffset+rts*5)*tupSize];
                y[tensorOffset*tupSize]          = (y1 + ((ru1*y2) % q) + ((ru2*y3) % q) + ((ru3*y4) % q) + ((ru4*y5) % q) + ((ru5*y6) % q)) % q;
                y[(tensorOffset+rts)*tupSize]      = (y1 + ((ru2*y2) % q) + ((ru4*y3) % q) + ((ru6*y4) % q) + ((ru1*y5) % q) + ((ru3*y6) % q)) % q;
                y[(tensorOffset+(rts<<1))*tupSize] = (y1 + ((ru3*y2) % q) + ((ru6*y3) % q) + ((ru2*y4) % q) + ((ru5*y5) % q) + ((ru1*y6) % q)) % q;
                y[(tensorOffset+rts*3)*tupSize]    = (y1 + ((ru4*y2) % q) + ((ru1*y3) % q) + ((ru5*y4) % q) + ((ru2*y5) % q) + ((ru6*y6) % q)) % q;
                y[(tensorOffset+(rts<<2))*tupSize] = (y1 + ((ru5*y2) % q) + ((ru3*y3) % q) + ((ru1*y4) % q) + ((ru6*y5) % q) + ((ru4*y6) % q)) % q;
                y[(tensorOffset+rts*5)*tupSize]    = (y1 + ((ru6*y2) % q) + ((ru5*y3) % q) + ((ru4*y4) % q) + ((ru3*y5) % q) + ((ru2*y6) % q)) % q;
            }
        }
    }
    else
    {
        hInt_t* tempSpace = (hInt_t*)malloc((p-1)*sizeof(hInt_t));
        hDim_t temp1 = rts*(p-1);
        for(hDim_t blockOffset = 0; blockOffset < lts; blockOffset++)
        {
            hDim_t temp2 = blockOffset*temp1;
            for(hDim_t modOffset = 0; modOffset < rts; modOffset++)
            {
                tensorOffset = temp2 + modOffset;
                
                for(hDim_t row = 1; row < p; row++)
                {
                    hInt_t acc = 0;
                    for(hDim_t col = 0; col < p-1; col++)
                    {
                        acc += ((y[(tensorOffset+col*rts)*tupSize]*ru[((col*row) % p)*rustride*tupSize]) % q);
                    }
                    tempSpace[row-1] = acc % q;
                }
                
                for(hDim_t row = 0; row < p-1; row++)
                {
                    y[(tensorOffset+rts*row)*tupSize] = tempSpace[row];
                }
            }
        }
        free(tempSpace);
    }
}

//takes inverse rus
void crtpinvRq (hInt_t* y, hShort_t tupSize, hDim_t lts, hDim_t rts, 
                hDim_t p, hDim_t rustride, hInt_t* ruinv, hInt_t q)
{
    if(p ==2)
    {
        // need this case so that we can divide overall by mhat^(-1)
        return;
    }
    else
    {
        hDim_t tensorOffset,i;
        hInt_t* tempSpace = (hInt_t*)malloc((p-1)*sizeof(hInt_t));
        hDim_t temp1 = rts*(p-1);
        for(hDim_t blockOffset = 0; blockOffset < lts; blockOffset++)
        {
            hDim_t temp2 = blockOffset*temp1;
            for(hDim_t modOffset = 0; modOffset < rts; modOffset++)
            {
                tensorOffset = temp2 + modOffset;
                
                for(i = 0; i < p-1; i++)
                {
                    hInt_t sum = 0;
                    int j;
                    for(j = 0; j < p-1; j++)
                    {
                        int ruIdx = ((j+1)*i) % p;
                        sum += ((y[(tensorOffset+j*rts)*tupSize] * ruinv[ruIdx*rustride*tupSize]) % q);
                    }
                    tempSpace[i] = sum % q;
                }

                hInt_t shift = 0;
                for(i = 0; i < p-1; i++)
                {
                    // we were given the inverse rus, so we need to negate the indices
                    shift += ((y[(tensorOffset+i*rts)*tupSize] * ruinv[rustride*(p-(i+1))*tupSize]) % q);
                }

                for(i = 0; i < p-1; i++)
                {
                    y[(tensorOffset+i*rts)*tupSize] = (tempSpace[i] - shift) % q; 
                }
            }
        }
    }
}

void ppDFTRq (hInt_t* y, hShort_t tupSize, hDim_t lts, hDim_t rts, 
              PrimeExponent pe, hDim_t rustride, hInt_t* ru, hInt_t q, hInt_t* temp)
{
    hDim_t p = pe.prime;
    hShort_t e = pe.exponent;
    
    if(e == 0)
    {
        return;
    }
    
    hDim_t primeRuStride = rustride*ipow(p,e-1);
    
    hShort_t i;
    
    hDim_t ltsScale = ipow(p,e-1);
    hDim_t rtsScale = 1;
    hDim_t twidRuStride = rustride;
    for(i = 0; i < e; i++)
    {
        hDim_t rtsDim = rts*rtsScale;
        dftpRq (y, tupSize, lts*ltsScale, rtsDim, p, primeRuStride, ru, temp, q);
        dftptwidRq (y, tupSize, lts, rtsDim, pe, ltsScale*p, twidRuStride, ru, q);
        
        ltsScale /= p;
        rtsScale *= p;
        twidRuStride *= p;
        pe.exponent -= 1;
    }
}

void ppDFTInvRq (hInt_t* y, hShort_t tupSize, hDim_t lts, hDim_t rts, 
                 PrimeExponent pe, hDim_t rustride, hInt_t* ru, hInt_t q, hInt_t* temp)
{
    hDim_t p = pe.prime;
    hShort_t e = pe.exponent;
    
    if(e == 0)
    {
        return;
    }
    hDim_t primeRuStride = rustride*ipow(p,e-1);

    hShort_t i;
    
    hDim_t ltsScale = 1;
    hDim_t rtsScale = ipow(p,e-1);
    hDim_t twidRuStride = primeRuStride;
    pe.exponent = 1;
    for(i = 0; i < e; i++)
    {
        hDim_t rtsDim = rts*rtsScale;
        hDim_t ltsScaleP = ltsScale*p;
        dftptwidRq (y, tupSize, lts, rtsDim, pe, ltsScaleP, twidRuStride, ru, q);
        dftpRq (y, tupSize, lts*ltsScale, rtsDim, p, primeRuStride, ru, temp, q);
        
        ltsScale = ltsScaleP;
        rtsScale /= p;
        twidRuStride /= p;
        pe.exponent += 1;
    }
}

void ppcrtRq (void* y, hShort_t tupSize, hDim_t lts, hDim_t rts, 
              PrimeExponent pe, void* ru, hInt_t* qs)
{
    hDim_t p = pe.prime;
    hDim_t e = pe.exponent;
#ifdef DEBUG_MODE
    ASSERT(e != 0);
#endif
    hDim_t mprime = ipow(p,e-1);
    
#ifdef DEBUG_MODE
    printf("lts is %" PRId32 "\trts is %" PRId32 "\n", lts, rts);
    printf("rus for p=%" PRId32 ", e=%" PRId16 "\t[", pe.prime, pe.exponent);
    hDim_t i;
    for(i = 0; i < ipow(p,e); i++) {
        printf("%" PRId64 ",", ((hInt_t*)ru)[i*tupSize+tupIdx]);
    }
    printf("]\n");
#endif

    hInt_t* temp = 0;
    if(p >= DFTP_GENERIC_SIZE)
    {
        temp = (hInt_t*)malloc(p*sizeof(hInt_t));
    }
    for(int tupIdx = 0; tupIdx < tupSize; tupIdx++) {
        hInt_t* z = ((hInt_t*)y)+tupIdx;
        hInt_t q = qs[tupIdx];
        hInt_t* ruOffset = ((hInt_t*)ru)+tupIdx;
        crtpRq (z, tupSize, lts*mprime, rts, p, mprime, ruOffset, q);
        crtTwiddleRq (z, tupSize, lts, rts, pe, ruOffset, q);
        pe.exponent -= 1;
        ppDFTRq (z,  tupSize, lts, rts*(p-1), pe, p, ruOffset, q, temp);
        pe.exponent += 1;
    }

    if(p >= DFTP_GENERIC_SIZE)
    {
        free(temp);
    }
}

void ppcrtinvRq (void* y, hShort_t tupSize, hDim_t lts, hDim_t rts, 
                 PrimeExponent pe, void* ru, hInt_t* qs)
{
    hDim_t p = pe.prime;
    hDim_t e = pe.exponent;
#ifdef DEBUG_MODE
    ASSERT(e != 0);
#endif
    hDim_t mprime = ipow(p,e-1);
#ifdef DEBUG_MODE
    printf("lts is %" PRId32 "\trts is %" PRId32 "\n", lts, rts);
    printf("rus for p=%" PRId32 ", e=%" PRId16 "\t[", pe.prime, pe.exponent);
    hDim_t i;
    for(i = 0; i < ipow(p,e); i++) {
        printf("%" PRId64 ",", ((hInt_t*)ru)[i*tupSize+tupIdx]);
    }
    printf("]\n");
#endif

    hInt_t* temp = 0;
    if(p >= DFTP_GENERIC_SIZE)
    {
        temp = (hInt_t*)malloc(p*sizeof(hInt_t));
    }
    for(int tupIdx = 0; tupIdx < tupSize; tupIdx++) {
        hInt_t* z = ((hInt_t*)y)+tupIdx;
        hInt_t q = qs[tupIdx];
        hInt_t* ruOffset = ((hInt_t*)ru)+tupIdx;

        pe.exponent -= 1;
        ppDFTInvRq (z, tupSize, lts, rts*(p-1), pe, p, ruOffset, q, temp);
        pe.exponent += 1;
        crtTwiddleRq (z, tupSize, lts, rts, pe, ruOffset, q);
        crtpinvRq (z, tupSize, lts*mprime, rts, p, mprime, ruOffset, q);
    }

    if(p >= DFTP_GENERIC_SIZE)
    {
        free(temp);
    }
}

// EAC: Somebody who knows C/C++ should find a better way to handle pointers-to-pointers in a generic way
void tensorCRTRq (hShort_t tupSize, hInt_t* y, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE, hInt_t** ru, hInt_t* qs)
{
    hDim_t i;
#ifdef STATS
    struct timespec s1,s2,s3,s4,t1,t2,t3,t4;

    crtRqCtr++;

    clock_gettime(CLOCK_REALTIME, &s1);
    clock_gettime(CLOCK_MONOTONIC, &s2);
    clock_gettime(CLOCK_PROCESS_CPUTIME_ID, &s3);
    clock_gettime(CLOCK_THREAD_CPUTIME_ID, &s4);
#endif
#ifdef DEBUG_MODE
    printf("\n\nEntered tensorCRTRq\ttotm=%" PRId32 "\tnumFacts=%" PRId16 "\ttupSize=%" PRId16 "\n[", totm, sizeOfPE,tupSize);
    for(i = 0; i < tupSize; i++) {
        printf("%" PRId64 ",", qs[i]);
    }
    printf("]\n[");
    for(i = 0; i < totm*tupSize; i++) {
        printf("%" PRId64 ",", y[i]);
    }
    printf("]\n[");
    for(i = 0; i < sizeOfPE; i++) {
        printf("(%" PRId32 ",%" PRId16 "),", peArr[i].prime, peArr[i].exponent);
    }
    printf("]\n");
#endif

    void** rus = (void**)malloc(sizeOfPE*sizeof(void*));
    for(i = 0; i < sizeOfPE; i++)
    {
        rus[i] = (void*) (ru[i]);
    }

  tensorFuserCRT (y, tupSize, ppcrtRq, totm, peArr, sizeOfPE, rus, qs);

    for(i = 0; i < tupSize; i++) {
        hInt_t q = qs[i];
      for(hDim_t j = 0; j < totm; j++)
      {
          if(y[j*tupSize+i]<0)
          {
              y[j*tupSize+i]+=q;
          }
#ifdef DEBUG_MODE
          if(y[j*tupSize+i]<0)
          {
              printf("TENSOR CRT^T INV\n");
          }
#endif
        }
  }

  free(rus);
#ifdef STATS
    clock_gettime(CLOCK_REALTIME, &t1);
    clock_gettime(CLOCK_MONOTONIC, &t2);
    clock_gettime(CLOCK_PROCESS_CPUTIME_ID, &t3);
    clock_gettime(CLOCK_THREAD_CPUTIME_ID, &t4);

    crttime1 = tsAdd(crttime1, tsSubtract(t1,s1));
    crttime2 = tsAdd(crttime2, tsSubtract(t2,s2));
    crttime3 = tsAdd(crttime3, tsSubtract(t3,s3));
    crttime4 = tsAdd(crttime4, tsSubtract(t4,s4));
#endif
}

//takes inverse rus
void tensorCRTInvRq (hShort_t tupSize, hInt_t* y, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE, 
                    hInt_t** ruinv, hInt_t* mhatInv, hInt_t* qs)
{
  hDim_t i;
#ifdef STATS
    struct timespec s1,t1;
    crtInvRqCtr++;
    clock_gettime(CLOCK_REALTIME, &s1);
#endif
#ifdef DEBUG_MODE
    printf("\n\nEntered tensorCRTInvRq\ttotm=%" PRId32 "\tnumFacts=%" PRId16 "\ttupSize=%" PRId16 "\nqs[", totm, sizeOfPE, tupSize);
    for(i = 0; i < tupSize; i++) {
        printf("%" PRId64 ",", qs[i]);
    }
    printf("]\nmhatInv[");
    for(i = 0; i < tupSize; i++) {
        printf("%" PRId64 ",", mhatInv[i]);
    }
    printf("]\ny[");
    for(i = 0; i < totm*tupSize; i++) {
        printf("%" PRId64 ",", y[i]);
    }
    printf("]\npps[");
    for(i = 0; i < sizeOfPE; i++) {
        printf("(%" PRId32 ",%" PRId16 "),", peArr[i].prime, peArr[i].exponent);
    }
    printf("]\n");
#endif

  void** rus = (void**)malloc(sizeOfPE*sizeof(void*));
    for(i = 0; i < sizeOfPE; i++)
    {
        rus[i] = (void*) (ruinv[i]);
    }
  tensorFuserCRT (y, tupSize, ppcrtinvRq, totm, peArr, sizeOfPE, rus, qs);

    for (i = 0; i < tupSize; i++) {
        hInt_t q = qs[i];
      for (hDim_t j = 0; j < totm; j++)
      {
          y[j*tupSize+i] = (y[j*tupSize+i]*mhatInv[i])%q;
          if(y[j*tupSize+i] < 0)
          {
              y[j*tupSize+i] +=q;
          }
#ifdef DEBUG_MODE
          if(y[j*tupSize+i]<0)
          {
              printf("TENSOR CRT INV\n");
          }
#endif
        }
  }

  free(rus);
#ifdef STATS
    clock_gettime(CLOCK_REALTIME, &t1);
    crtInvRqTime = tsAdd(crtInvRqTime, tsSubtract(t1,s1));
#endif
}








































void crtTwiddleC (complex_t* y, hShort_t tupSize, hDim_t lts, hDim_t rts, PrimeExponent pe, complex_t* ru)
{
    hDim_t idx;
    hDim_t p = pe.prime;
    hShort_t e = pe.exponent;
    
#ifdef DEBUG_MODE
    ASSERT(e != 0);
#endif

    pe.exponent -= 1; // used for an argument to bitrev
    
    if(p == 2)
    {
        hDim_t mprime = 1<<(e-1);
        hDim_t blockDim = rts*mprime; // size of block in block diagonal tensor matrix

        for(hDim_t i0 = 1; i0 < mprime; i0++) // loops over i/(p-1) for i = 0..(m'-1), we can skip i0 = 0
        {
            hDim_t temp2 = i0*rts;
            complex_t twid = ru[bitrev(pe,i0)*tupSize];

            for(hDim_t blockIdx = 0; blockIdx < lts; blockIdx++)
            {
                hDim_t temp3 = blockIdx*blockDim + temp2;
                for(hDim_t modOffset = 0; modOffset < rts; modOffset++)
                {
                    idx = temp3 + modOffset;
                    CMPLX_IMUL(y[idx*tupSize],twid);
                }
            }
        }
    }
    else
    {
        hDim_t mprime = ipow(p,e-1);
        hDim_t blockDim = rts*(p-1)*mprime; // size of block in block diagonal tensor matrix

        for(hDim_t i0 = 1; i0 < mprime; i0++) // loops over i/(p-1) for i = 0..(m'-1), we can skip i0 = 0
        {
            hDim_t temp1 = i0*(p-1);
            for(hDim_t i1 = 0; i1 < (p-1); i1++) // loops over i%(p-1) for i = 0..(m'-1)
            {        
                hDim_t temp2 = (temp1+i1)*rts;
                complex_t twid = ru[bitrev(pe,i0)*(i1+1)*tupSize];

                for(hDim_t blockIdx = 0; blockIdx < lts; blockIdx++)
                {
                    hDim_t temp3 = blockIdx*blockDim + temp2;
                    for(hDim_t modOffset = 0; modOffset < rts; modOffset++)
                    {
                        idx = temp3 + modOffset;
                        CMPLX_IMUL(y[idx*tupSize],twid);
                    }
                }
            }
        }
    }
}
    
// dim is power of p
void dftptwidC (complex_t* y, hShort_t tupSize, hDim_t lts, hDim_t rts, PrimeExponent pe, hDim_t dim, hDim_t rustride, complex_t* ru)
{
    hDim_t idx;
    hDim_t p = pe.prime;
    pe.exponent -= 1; // used for an argument to bitrev
    
    if(p == 2)
    {
        hDim_t mprime = dim>>1; // divides evenly
        hDim_t temp1 = rts*dim; // for use in computing [modified] tensorOffset
        for(hDim_t i0 = 1; i0 < mprime; i0++) // loops over i/p for i = 0..(dim-1), but we skip i0=0
        {
            hDim_t temp3 = rts*(i0*p+1);
            complex_t twid = ru[bitrev(pe,i0)*rustride*tupSize];

            for(hDim_t blockOffset = 0; blockOffset < lts; blockOffset++)
            {
                hDim_t temp2 = blockOffset*temp1 + temp3;
                for(hDim_t modOffset = 0; modOffset < rts; modOffset++)
                {
                    idx = temp2 + modOffset;
                    CMPLX_IMUL(y[idx*tupSize],twid);
                }
            }
        }
    }
    else
    {
        hDim_t mprime = dim/p; // divides evenly
        hDim_t temp1 = rts*dim; // for use in computing [modified] tensorOffset
        for(hDim_t i0 = 1; i0 < mprime; i0++) // loops over i/p for i = 0..(dim-1), but we skip i0=0
        {
            for(hDim_t i1 = 1; i1 < p; i1++) // loops over i%p for i = 0..(dim-1), but we skip i1=0
            {
                hDim_t temp3 = rts*(i0*p+i1);
                complex_t twid = ru[bitrev(pe,i0)*i1*rustride*tupSize];

                for(hDim_t blockOffset = 0; blockOffset < lts; blockOffset++)
                {
                    hDim_t temp2 = blockOffset*temp1 + temp3;
                    for(hDim_t modOffset = 0; modOffset < rts; modOffset++)
                    {
                        idx = temp2 + modOffset;
                        CMPLX_IMUL(y[idx*tupSize],twid);
                    }
                }
            }
        }
    }
}

//implied length of ru is rustride*p
//implied length of tempSpace is p, if p is not a special case
void dftpC (complex_t* y, hShort_t tupSize, hDim_t lts, hDim_t rts, hDim_t p, hDim_t rustride, complex_t* ru, complex_t* tempSpace)
{
    hDim_t blockOffset, modOffset, tensorOffset;
    
    if(p == 2)
    {
        hDim_t temp1 = rts<<1;
        for(blockOffset = 0; blockOffset < lts; blockOffset++)
        {
            hDim_t temp2 = blockOffset*temp1;
            for(modOffset = 0; modOffset < rts; modOffset++)
            {
                tensorOffset = temp2 + modOffset;
                complex_t u = y[tensorOffset*tupSize];
                complex_t t = y[(tensorOffset+rts)*tupSize];
                y[tensorOffset*tupSize] = CMPLX_ADD(u,t);
                y[(tensorOffset+rts)*tupSize] = CMPLX_SUB(u,t);
            }
        }
    }
    else if(p == 3)
    {
        hDim_t temp1 = rts*3;
        complex_t ru1 = ru[rustride*tupSize];
        complex_t ru2 = ru[(rustride<<1)*tupSize];

        for(blockOffset = 0; blockOffset < lts; blockOffset++)
        {
            hDim_t temp2 = blockOffset*temp1;
            for(modOffset = 0; modOffset < rts; modOffset++)
            {
                tensorOffset = temp2 + modOffset;
                complex_t y1, y2, y3;
                y1 = y[tensorOffset*tupSize];
                y2 = y[(tensorOffset+rts)*tupSize];
                y3 = y[(tensorOffset+(rts<<1))*tupSize];
                y[tensorOffset*tupSize]            = CMPLX_ADD3(y1,               y2,                y3);
                y[(tensorOffset+rts)*tupSize]      = CMPLX_ADD3(y1, CMPLX_MUL(ru1,y2), CMPLX_MUL(ru2,y3));
                y[(tensorOffset+(rts<<1))*tupSize] = CMPLX_ADD3(y1, CMPLX_MUL(ru2,y2), CMPLX_MUL(ru1,y3));
            }   
        }
    }
    else if(p == 5)
    {
        hDim_t temp1 = rts*5;
        complex_t ru1 = ru[rustride*tupSize];
        complex_t ru2 = ru[(rustride<<1)*tupSize];
        complex_t ru3 = ru[(rustride*3)*tupSize];
        complex_t ru4 = ru[(rustride<<2)*tupSize];

        for(blockOffset = 0; blockOffset < lts; blockOffset++)
        {
            hDim_t temp2 = blockOffset*temp1;
            for(modOffset = 0; modOffset < rts; modOffset++)
            {
                tensorOffset = temp2 + modOffset;
                complex_t y1, y2, y3, y4, y5;
                y1 = y[tensorOffset*tupSize];
                y2 = y[(tensorOffset+rts)*tupSize];
                y3 = y[(tensorOffset+(rts<<1))*tupSize];
                y4 = y[(tensorOffset+3*rts)*tupSize];
                y5 = y[(tensorOffset+(rts<<2))*tupSize];
                y[tensorOffset*tupSize]          = CMPLX_ADD5(y1,               y2,                y3,                y4,                y5);
                y[(tensorOffset+rts)*tupSize]      = CMPLX_ADD5(y1, CMPLX_MUL(ru1,y2), CMPLX_MUL(ru2,y3), CMPLX_MUL(ru3,y4), CMPLX_MUL(ru4,y5));
                y[(tensorOffset+(rts<<1))*tupSize] = CMPLX_ADD5(y1, CMPLX_MUL(ru2,y2), CMPLX_MUL(ru4,y3), CMPLX_MUL(ru1,y4), CMPLX_MUL(ru3,y5));
                y[(tensorOffset+rts*3)*tupSize]    = CMPLX_ADD5(y1, CMPLX_MUL(ru3,y2), CMPLX_MUL(ru1,y3), CMPLX_MUL(ru4,y4), CMPLX_MUL(ru2,y5));
                y[(tensorOffset+(rts<<2))*tupSize] = CMPLX_ADD5(y1, CMPLX_MUL(ru4,y2), CMPLX_MUL(ru3,y3), CMPLX_MUL(ru2,y4), CMPLX_MUL(ru1,y5));
            }   
        }
    }
    else if(p == 7)
    {
        hDim_t temp1 = rts*7;
        complex_t ru1 = ru[rustride*tupSize];
        complex_t ru2 = ru[(rustride<<1)*tupSize];
        complex_t ru3 = ru[(rustride*3)*tupSize];
        complex_t ru4 = ru[(rustride<<2)*tupSize];
        complex_t ru5 = ru[(rustride*5)*tupSize];
        complex_t ru6 = ru[(rustride*6)*tupSize];

        for(blockOffset = 0; blockOffset < lts; blockOffset++)
        {
            hDim_t temp2 = blockOffset*temp1;
            for(modOffset = 0; modOffset < rts; modOffset++)
            {
                tensorOffset = temp2 + modOffset;
                complex_t y1, y2, y3, y4, y5, y6, y7;
                y1 = y[tensorOffset*tupSize];
                y2 = y[(tensorOffset+rts)*tupSize];
                y3 = y[(tensorOffset+(rts<<1))*tupSize];
                y4 = y[(tensorOffset+3*rts)*tupSize];
                y5 = y[(tensorOffset+(rts<<2))*tupSize];
                y6 = y[(tensorOffset+rts*5)*tupSize];
                y7 = y[(tensorOffset+rts*6)*tupSize];
                y[tensorOffset*tupSize]          = CMPLX_ADD7(y1,               y2,                y3,                y4,                y5,                y6,                y7);
                y[(tensorOffset+rts)*tupSize]      = CMPLX_ADD7(y1, CMPLX_MUL(ru1,y2), CMPLX_MUL(ru2,y3), CMPLX_MUL(ru3,y4), CMPLX_MUL(ru4,y5), CMPLX_MUL(ru5,y6), CMPLX_MUL(ru6,y7));
                y[(tensorOffset+(rts<<1))*tupSize] = CMPLX_ADD7(y1, CMPLX_MUL(ru2,y2), CMPLX_MUL(ru4,y3), CMPLX_MUL(ru6,y4), CMPLX_MUL(ru1,y5), CMPLX_MUL(ru3,y6), CMPLX_MUL(ru5,y7));
                y[(tensorOffset+rts*3)*tupSize]    = CMPLX_ADD7(y1, CMPLX_MUL(ru3,y2), CMPLX_MUL(ru6,y3), CMPLX_MUL(ru2,y4), CMPLX_MUL(ru5,y5), CMPLX_MUL(ru1,y6), CMPLX_MUL(ru4,y7));
                y[(tensorOffset+(rts<<2))*tupSize] = CMPLX_ADD7(y1, CMPLX_MUL(ru4,y2), CMPLX_MUL(ru1,y3), CMPLX_MUL(ru5,y4), CMPLX_MUL(ru2,y5), CMPLX_MUL(ru6,y6), CMPLX_MUL(ru3,y7));
                y[(tensorOffset+rts*5)*tupSize]    = CMPLX_ADD7(y1, CMPLX_MUL(ru5,y2), CMPLX_MUL(ru3,y3), CMPLX_MUL(ru1,y4), CMPLX_MUL(ru6,y5), CMPLX_MUL(ru4,y6), CMPLX_MUL(ru2,y7));
                y[(tensorOffset+rts*6)*tupSize]    = CMPLX_ADD7(y1, CMPLX_MUL(ru6,y2), CMPLX_MUL(ru5,y3), CMPLX_MUL(ru4,y4), CMPLX_MUL(ru3,y5), CMPLX_MUL(ru2,y6), CMPLX_MUL(ru1,y7));
            }   
        }
    }
    else
    {
        hDim_t temp1 = rts*p;
        for(blockOffset = 0; blockOffset < lts; blockOffset++)
        {
            hDim_t temp2 = blockOffset*temp1;
            for(modOffset = 0; modOffset < rts; modOffset++)
            {
                tensorOffset = temp2 + modOffset;
                hDim_t row, col;
                
                for(row = 0; row < p; row++)
                {
                    complex_t acc = ((complex_t){0,0});
                    for(col = 0; col < p; col++)
                    {
                        CMPLX_IADD(acc, CMPLX_MUL(y[(tensorOffset+col*rts)*tupSize], ru[((col*row) % p)*rustride*tupSize]));
                    }
                    tempSpace[row] = acc;
                }
                
                for(row = 0; row < p; row++)
                {
                    y[(tensorOffset+rts*row)*tupSize] = tempSpace[row];   
                }
            }
        }
    }
}

void crtpC (complex_t* y, hShort_t tupSize, hDim_t lts, hDim_t rts, hDim_t p, hDim_t rustride, complex_t* ru)
{
    hDim_t blockOffset, modOffset, tensorOffset;
    
    if(p == 2)
    {
        return;
    }
    else if(p == 3)
    {
        hDim_t temp1 = rts*2;
        complex_t ru1 = ru[rustride*tupSize];
        complex_t ru2 = ru[(rustride<<1)*tupSize];

        for(blockOffset = 0; blockOffset < lts; blockOffset++)
        {
            hDim_t temp2 = blockOffset*temp1;
            for(modOffset = 0; modOffset < rts; modOffset++)
            {
                tensorOffset = temp2 + modOffset;
                complex_t y1, y2;
                y1 = y[tensorOffset*tupSize];
                y2 = y[(tensorOffset+rts)*tupSize];
                y[tensorOffset*tupSize]     = CMPLX_ADD(y1, CMPLX_MUL(ru1,y2));
                y[(tensorOffset+rts)*tupSize] = CMPLX_ADD(y1, CMPLX_MUL(ru2,y2));
            }   
        }
    }
    else if(p == 5)
    {
        hDim_t temp1 = rts*4;
        complex_t ru1 = ru[rustride*tupSize];
        complex_t ru2 = ru[(rustride<<1)*tupSize];
        complex_t ru3 = ru[(rustride*3)*tupSize];
        complex_t ru4 = ru[(rustride<<2)*tupSize];

        for(blockOffset = 0; blockOffset < lts; blockOffset++)
        {
            hDim_t temp2 = blockOffset*temp1;
            for(modOffset = 0; modOffset < rts; modOffset++)
            {
                tensorOffset = temp2 + modOffset;
                complex_t y1, y2, y3, y4;
                y1 = y[tensorOffset*tupSize];
                y2 = y[(tensorOffset+rts)*tupSize];
                y3 = y[(tensorOffset+(rts<<1))*tupSize];
                y4 = y[(tensorOffset+3*rts)*tupSize];
                y[tensorOffset*tupSize]          = CMPLX_ADD4(y1, CMPLX_MUL(ru1,y2), CMPLX_MUL(ru2,y3), CMPLX_MUL(ru3,y4));
                y[(tensorOffset+rts)*tupSize]      = CMPLX_ADD4(y1, CMPLX_MUL(ru2,y2), CMPLX_MUL(ru4,y3), CMPLX_MUL(ru1,y4));
                y[(tensorOffset+(rts<<1))*tupSize] = CMPLX_ADD4(y1, CMPLX_MUL(ru3,y2), CMPLX_MUL(ru1,y3), CMPLX_MUL(ru4,y4));
                y[(tensorOffset+rts*3)*tupSize]    = CMPLX_ADD4(y1, CMPLX_MUL(ru4,y2), CMPLX_MUL(ru3,y3), CMPLX_MUL(ru2,y4));
            }   
        }
    }
    else if(p == 7)
    {
        hDim_t temp1 = rts*6;
        complex_t ru1 = ru[rustride*tupSize];
        complex_t ru2 = ru[(rustride<<1)*tupSize];
        complex_t ru3 = ru[(rustride*3)*tupSize];
        complex_t ru4 = ru[(rustride<<2)*tupSize];
        complex_t ru5 = ru[(rustride*5)*tupSize];
        complex_t ru6 = ru[(rustride*6)*tupSize];

        for(blockOffset = 0; blockOffset < lts; blockOffset++)
        {
            hDim_t temp2 = blockOffset*temp1;
            for(modOffset = 0; modOffset < rts; modOffset++)
            {
                tensorOffset = temp2 + modOffset;
                complex_t y1, y2, y3, y4, y5, y6;
                y1 = y[tensorOffset*tupSize];
                y2 = y[(tensorOffset+rts)*tupSize];
                y3 = y[(tensorOffset+(rts<<1))*tupSize];
                y4 = y[(tensorOffset+3*rts)*tupSize];
                y5 = y[(tensorOffset+(rts<<2))*tupSize];
                y6 = y[(tensorOffset+rts*5)*tupSize];
                y[tensorOffset*tupSize]          = CMPLX_ADD6(y1, CMPLX_MUL(ru1,y2), CMPLX_MUL(ru2,y3), CMPLX_MUL(ru3,y4), CMPLX_MUL(ru4,y5), CMPLX_MUL(ru5,y6));
                y[(tensorOffset+rts)*tupSize]      = CMPLX_ADD6(y1, CMPLX_MUL(ru2,y2), CMPLX_MUL(ru4,y3), CMPLX_MUL(ru6,y4), CMPLX_MUL(ru1,y5), CMPLX_MUL(ru3,y6));
                y[(tensorOffset+(rts<<1))*tupSize] = CMPLX_ADD6(y1, CMPLX_MUL(ru3,y2), CMPLX_MUL(ru6,y3), CMPLX_MUL(ru2,y4), CMPLX_MUL(ru5,y5), CMPLX_MUL(ru1,y6));
                y[(tensorOffset+rts*3)*tupSize]    = CMPLX_ADD6(y1, CMPLX_MUL(ru4,y2), CMPLX_MUL(ru1,y3), CMPLX_MUL(ru5,y4), CMPLX_MUL(ru2,y5), CMPLX_MUL(ru6,y6));
                y[(tensorOffset+(rts<<2))*tupSize] = CMPLX_ADD6(y1, CMPLX_MUL(ru5,y2), CMPLX_MUL(ru3,y3), CMPLX_MUL(ru1,y4), CMPLX_MUL(ru6,y5), CMPLX_MUL(ru4,y6));
                y[(tensorOffset+rts*5)*tupSize]    = CMPLX_ADD6(y1, CMPLX_MUL(ru6,y2), CMPLX_MUL(ru5,y3), CMPLX_MUL(ru4,y4), CMPLX_MUL(ru3,y5), CMPLX_MUL(ru2,y6));
            }   
        }
    }
    else
    {
        complex_t* tempSpace = (complex_t*)malloc((p-1)*sizeof(complex_t));
        hDim_t temp1 = rts*(p-1);
        for(blockOffset = 0; blockOffset < lts; blockOffset++)
        {
            hDim_t temp2 = blockOffset*temp1;
            for(modOffset = 0; modOffset < rts; modOffset++)
            {
                tensorOffset = temp2 + modOffset;
                hDim_t row, col;
                
                for(row = 1; row < p; row++)
                {
                    complex_t acc = ((complex_t){0,0});
                    for(col = 0; col < p-1; col++)
                    {
                        CMPLX_IADD(acc, CMPLX_MUL(y[(tensorOffset+col*rts)*tupSize], ru[((col*row) % p)*rustride*tupSize]));
                    }
                    tempSpace[row-1] = acc;
                }
                
                for(row = 0; row < p-1; row++)
                {
                    y[(tensorOffset+rts*row)*tupSize] = tempSpace[row];   
                }
            }
        }
        free(tempSpace);
    }
}

//takes inverse rus
void crtpinvC (complex_t* y, hShort_t tupSize, hDim_t lts, hDim_t rts, hDim_t p, hDim_t rustride, complex_t* ruinv)
{
    if(p ==2)
    {
        // need this case so that we can divide overall by mhat^(-1)
        return;
    }
    else
    {
        hDim_t tensorOffset,i;
        complex_t* tempSpace = (complex_t*)malloc(p*sizeof(complex_t));
        hDim_t temp1 = rts*(p-1);
        for(hDim_t blockOffset = 0; blockOffset < lts; blockOffset++)
        {
            hDim_t temp2 = blockOffset*temp1;
            for(hDim_t modOffset = 0; modOffset < rts; modOffset++)
            {
                tensorOffset = temp2 + modOffset;

                for(i = 0; i < p-1; i++)
                {
                    complex_t sum = ((complex_t){0,0});
                    int j;
                    for(j = 0; j < p-1; j++)
                    {
                        int ruIdx = (((j+1)*i) % p)*rustride;
                        CMPLX_IADD(sum, CMPLX_MUL(y[(tensorOffset+j*rts)*tupSize],ruinv[ruIdx*tupSize]));
                    }
                    tempSpace[i] = sum;
                }

                complex_t shift = ((complex_t){0,0});
                for(i = 0; i < p-1; i++)
                {
                    // we were given the inverse rus, so we need to negate the indices
                    int ruIdx = p-(i+1);
                    CMPLX_IADD(shift, CMPLX_MUL(y[(tensorOffset+i*rts)*tupSize], ruinv[rustride*ruIdx*tupSize]));
                }

                for(i = 0; i < p-1; i++)
                {
                    y[(tensorOffset+i*rts)*tupSize] = CMPLX_SUB(tempSpace[i], shift); 
                }
            }
        }
    }
}

void ppDFTC (complex_t* y, hShort_t tupSize, hDim_t lts, hDim_t rts, PrimeExponent pe, hDim_t rustride, complex_t* ru, complex_t* temp)
{
    hDim_t p = pe.prime;
    hShort_t e = pe.exponent;
    
    if(e == 0)
    {
        return;
    }
    
    hDim_t primeRuStride = rustride*ipow(p,e-1);    

    hShort_t i;
    
    hDim_t ltsScale = ipow(p,e-1);
    hDim_t rtsScale = 1;
    hDim_t twidRuStride = rustride;
    for(i = 0; i < e; i++)
    {
        hDim_t rtsDim = rts*rtsScale;
        dftpC (y, tupSize, lts*ltsScale, rtsDim, p, primeRuStride, ru, temp);
        dftptwidC (y, tupSize, lts, rtsDim, pe, ltsScale*p, twidRuStride, ru);
        
        ltsScale /= p;
        rtsScale *= p;
        twidRuStride *= p;
        pe.exponent -= 1;
    }
}

void ppDFTInvC (complex_t* y, hShort_t tupSize, hDim_t lts, hDim_t rts, PrimeExponent pe, hDim_t rustride, complex_t* ru, complex_t* temp)
{

    hDim_t p = pe.prime;
    hShort_t e = pe.exponent;
    
    if(e == 0)
    {
        return;
    }
    hDim_t primeRuStride = rustride*ipow(p,e-1);

    hShort_t i;
    
    hDim_t ltsScale = 1;
    hDim_t rtsScale = ipow(p,e-1);
    hDim_t twidRuStride = primeRuStride;
    pe.exponent = 1;
    for(i = 0; i < e; i++)
    {
        hDim_t rtsDim = rts*rtsScale;
        hDim_t ltsScaleP = ltsScale*p;
        dftptwidC (y, tupSize, lts, rtsDim, pe, ltsScaleP, twidRuStride, ru);
        dftpC (y, tupSize, lts*ltsScale, rtsDim, p, primeRuStride, ru, temp);
        
        ltsScale = ltsScaleP;
        rtsScale /= p;
        twidRuStride /= p;
        pe.exponent += 1;
    }
}

void ppcrtC (void* y, hShort_t tupSize, hDim_t lts, hDim_t rts, PrimeExponent pe, void* ru, hInt_t* qs)
{
    hDim_t p = pe.prime;
    hDim_t e = pe.exponent;
#ifdef DEBUG_MODE
    ASSERT(e != 0);
#endif
    hDim_t mprime = ipow(p,e-1);

#ifdef DEBUG_MODE
    printf("rus for p=%" PRId32 ", e=%" PRId16 "\t[", pe.prime, pe.exponent);
    hDim_t i;
    for(i = 0; i < ipow(p,e); i++) {
        printf("(%f,%f),", ((complex_t*)ru)[i].real, ((complex_t*)ru)[i].imag);
    }
    printf("]\n");
#endif

    complex_t* temp = 0;
    if(p >= DFTP_GENERIC_SIZE)
    {
        temp = (complex_t*)malloc(p*sizeof(complex_t));
    }

    for(int tupIdx = 0; tupIdx < tupSize; tupIdx++) {
        complex_t* z = ((complex_t*)y)+tupIdx;
        complex_t* ruOffset = ((complex_t*)ru)+tupIdx;
        crtpC (z, tupSize, lts*mprime, rts, p, mprime, ruOffset);
        crtTwiddleC (z, tupSize, lts, rts, pe, ruOffset);
        pe.exponent -= 1;
        ppDFTC (z, tupSize, lts, rts*(p-1), pe, p, ruOffset, temp);
        pe.exponent += 1;
    }

    if(p >= DFTP_GENERIC_SIZE)
    {
        free(temp);
    }
}

void ppcrtinvC (void* y, hShort_t tupSize, hDim_t lts, hDim_t rts, PrimeExponent pe, void* ru, hInt_t* qs)
{
    hDim_t p = pe.prime;
    hDim_t e = pe.exponent;
#ifdef DEBUG_MODE
    ASSERT(e != 0);
#endif
    hDim_t mprime = ipow(p,e-1);
    
    complex_t* temp = 0;
    if(p >= DFTP_GENERIC_SIZE)
    {
        temp = (complex_t*)malloc(p*sizeof(complex_t));
    }

    for(int tupIdx = 0; tupIdx < tupSize; tupIdx++) {
        complex_t* z = ((complex_t*)y)+tupIdx;
        complex_t* ruOffset = ((complex_t*)ru)+tupIdx;
        pe.exponent -= 1;
        ppDFTInvC (z, tupSize, lts, rts*(p-1), pe, p, ruOffset, temp);
        pe.exponent += 1;
        crtTwiddleC (z, tupSize, lts, rts, pe, ruOffset);
        crtpinvC (z, tupSize, lts*mprime, rts, p, mprime, ruOffset);
    }

    if(p >= DFTP_GENERIC_SIZE)
    {
        free(temp);
    }
}

void tensorCRTC (hShort_t tupSize, complex_t* y, hDim_t totm, PrimeExponent* peArr, hShort_t sizeOfPE, complex_t** ru)
{
#ifdef STATS
    struct timespec s1,t1;
    crtCCtr++;
    clock_gettime(CLOCK_REALTIME, &s1);
#endif
#ifdef DEBUG_MODE
    printf("\n\nEntered tensorCRTC\ttotm=%" PRId32 "\tnumFacts=%" PRId16 "\n[", totm, sizeOfPE);
    hDim_t j;
    for(j = 0; j < totm; j++) {
        printf("(%f,%f),", y[j].real, y[j].imag);
    }
    printf("]\n[");
    for(j = 0; j < sizeOfPE; j++) {
        printf("(%" PRId32 ",%" PRId16 "),", peArr[j].prime, peArr[j].exponent);
    }
    printf("]\n");
#endif
    void** rus = (void**)malloc(sizeOfPE*sizeof(void*));
    hShort_t i;
    for(i = 0; i < sizeOfPE; i++)
    {
        rus[i] = (void*) (ru[i]);
    }
  tensorFuserCRT (y, tupSize, ppcrtC, totm, peArr, sizeOfPE, rus, (hInt_t*)0);
  free(rus);
#ifdef STATS
    clock_gettime(CLOCK_REALTIME, &t1);
    crtCTime = tsAdd(crtCTime, tsSubtract(t1,s1));
#endif
}

//takes inverse rus
void tensorCRTInvC (hShort_t tupSize, complex_t* y, hDim_t totm, PrimeExponent* peArr, 
                    hShort_t sizeOfPE, complex_t** ruinv, complex_t* mhatInv)
{
#ifdef STATS
    struct timespec s1,t1;
    crtInvCCtr++;
    clock_gettime(CLOCK_REALTIME, &s1);
#endif
  hDim_t i;
  
  void** rus = (void**)malloc(sizeOfPE*sizeof(void*));
    for(i = 0; i < sizeOfPE; i++)
    {
        rus[i] = (void*) (ruinv[i]);
    }
  
  tensorFuserCRT (y, tupSize, ppcrtinvC, totm, peArr, sizeOfPE, rus, (hInt_t*)0);

    for (i = 0; i < tupSize; i++) {
      for (hDim_t j = 0; j < totm; j++)
      {
          CMPLX_IMUL(y[j*tupSize+i], mhatInv[i]);
      }
    }
  
  free(rus);
#ifdef STATS
    clock_gettime(CLOCK_REALTIME, &t1);
    crtInvCTime = tsAdd(crtInvCTime, tsSubtract(t1,s1));
#endif
}