/******************************************************** * Some code. Copyright (C) 2003 by Pascal Massimino.   * * All Rights Reserved.      (http://skal.planet-d.net) * * For Educational/Academic use ONLY. See 'LICENSE.TXT'.* ********************************************************//* *  skl_dct16_LLM.cpp * *  "Fast and precise" LLM implementation of FDCT/IDCT, where  *  rotations are decomposed using: *    tmp = (x+y).cos t *    x' = tmp + y.(sin t - cos t) *    y' = tmp - x.(sin t + cos t) * *  See details at http://skl.planet-d.net/coding/dct.html *  and at the end of this file... * * Reference (e.g.): *  Loeffler C., Ligtenberg A., and Moschytz C.S.:  *    Practical Fast 1D DCT Algorithm with Eleven Multiplications,  *  Proc. ICASSP 1989, 988-991. * *  IEEE-1180-like error specs for FDCT: * Peak error:   1.0000 * Peak MSE:     0.0340 * Overall MSE:  0.0200 * Peak ME:      0.0191 * Overall ME:   -0.0033 * *  error specs for IDCT: * Peak error:   1.0000  * Peak MSE:     0.0065 * Overall MSE:  0.0051 * Peak ME:      0.0015 * Overall ME:   0.0000 * ********************************************************/#include "skl.h"#include <math.h>extern "C" {typedef SKL_INT16 TYPE;//////////////////////////////////////////////////////////#define LOAD_BUTF(m1, m2, a, b, tmp, S) \  (m1) = (S)[(a)] + (S)[(b)];           \  (m2) = (S)[(a)] - (S)[(b)]#define BUTF(a, b, tmp) \  (tmp) = (a)+(b);      \  (b) = (a)-(b);        \  (a) = (tmp)#define ROTATE(m1,m2,c,k1,k2,tmp,Fix,Rnd) \  (tmp) = ( (m1) + (m2) )*(c);            \  (m1) *= k1;                             \  (m2) *= k2;                             \  (tmp) += (Rnd);                         \  (m1) = ((m1)+(tmp))>>(Fix);             \  (m2) = ((m2)+(tmp))>>(Fix);#define ROTATE2(m1,m2,c,k1,k2,tmp) \  (tmp) = ( (m1) + (m2) )*(c);     \  (m1) *= k1;                      \  (m2) *= k2;                      \  (m1) = (m1)+(tmp);               \  (m2) = (m2)+(tmp);#define ROTATE0(m1,m2,c,k1,k2,tmp) \  (m1) = ( (m2) )*(c);             \  (m2) = (m2)*k2+(m1);#define SHIFTL(x,n)   ((x)<<(n))#define SHIFTR(x, n)  ((x)>>(n))#define HALF(n)       (1<<((n)-1))#define IPASS 3#define FPASS 2#define FIX  16#if 1#define ROT6_C     35468#define ROT6_SmC   50159#define ROT6_SpC  121095#define ROT17_C    77062#define ROT17_SmC  25571#define ROT17_SpC 128553#define ROT37_C    58981#define ROT37_SmC  98391#define ROT37_SpC  19571#define ROT13_C   167963#define ROT13_SmC 134553#define ROT13_SpC 201373#else#define FX(x) ( (int)floor((x)*(1<<FIX) + .5 ) )static const double c1 = cos(1.*M_PI/16);static const double c2 = cos(2.*M_PI/16);static const double c3 = cos(3.*M_PI/16);static const double c4 = cos(4.*M_PI/16);static const double c5 = cos(5.*M_PI/16);static const double c6 = cos(6.*M_PI/16);static const double c7 = cos(7.*M_PI/16);static const int ROT6_C   = FX(c2-c6);  // 0.541static const int ROT6_SmC = FX(2*c6);   // 0.765static const int ROT6_SpC = FX(2*c2);   // 1.847static const int ROT17_C   = FX(c1+c7);  // 1.175static const int ROT17_SmC = FX(2*c7);   // 0.390static const int ROT17_SpC = FX(2*c1);   // 1.961static const int ROT37_C   = FX((c3-c7)/c4);  // 0.899static const int ROT37_SmC = FX(2*(c5+c7));   // 1.501static const int ROT37_SpC = FX(2*(c1-c3));   // 0.298static const int ROT13_C   = FX((c1+c3)/c4);  // 2.562static const int ROT13_SmC = FX(2*(c3+c7));   // 2.053static const int ROT13_SpC = FX(2*(c1+c5));   // 3.072#endif//////////////////////////////////////////////////////////void Skl_Dct16_C( TYPE *In ){  TYPE *pIn;  int i;  pIn = In;  for(i=8; i>0; --i)  {    int mm0, mm1, mm2, mm3, mm4, mm5, mm6, mm7, Spill;      // odd    LOAD_BUTF(mm1,mm6, 1, 6, mm0, pIn);    LOAD_BUTF(mm2,mm5, 2, 5, mm0, pIn);    LOAD_BUTF(mm3,mm4, 3, 4, mm0, pIn);    LOAD_BUTF(mm0,mm7, 0, 7, Spill, pIn);    BUTF(mm1, mm2, Spill);    BUTF(mm0, mm3, Spill);    ROTATE(mm3, mm2, ROT6_C, ROT6_SmC, -ROT6_SpC, Spill, FIX-FPASS, HALF(FIX-FPASS));    pIn[2] = mm3;    pIn[6] = mm2;    BUTF(mm0, mm1, Spill);    pIn[0] = SHIFTL(mm0, FPASS);    pIn[4] = SHIFTL(mm1, FPASS);      // even    mm3 = mm5 + mm7;    mm2 = mm4 + mm6;    ROTATE(mm2, mm3,  ROT17_C, -ROT17_SpC, -ROT17_SmC, mm0, FIX-FPASS, HALF(FIX-FPASS));    ROTATE(mm4, mm7, -ROT37_C,  ROT37_SpC,  ROT37_SmC, mm0, FIX-FPASS, HALF(FIX-FPASS));    mm7 += mm3;    mm4 += mm2;    pIn[1] = mm7;    pIn[7] = mm4;    ROTATE(mm5, mm6, -ROT13_C,  ROT13_SmC,  ROT13_SpC, mm0, FIX-FPASS, HALF(FIX-FPASS));    mm5 += mm3;    mm6 += mm2;    pIn[3] = mm6;    pIn[5] = mm5;    pIn  += 8;  }  pIn = In;  for(i=8; i>0; --i)  {    int mm0, mm1, mm2, mm3, mm4, mm5, mm6, mm7, Spill;      // odd    LOAD_BUTF(mm1,mm6, 1*8, 6*8, mm0, pIn);    LOAD_BUTF(mm2,mm5, 2*8, 5*8, mm0, pIn);    BUTF(mm1, mm2, mm0);    LOAD_BUTF(mm3,mm4, 3*8, 4*8, mm0, pIn);    LOAD_BUTF(mm0,mm7, 0*8, 7*8, Spill, pIn);    BUTF(mm0, mm3, Spill);    ROTATE(mm3, mm2, ROT6_C, ROT6_SmC, -ROT6_SpC, Spill, 0,  HALF(FIX+FPASS+3));    pIn[2*8] = (TYPE)SHIFTR(mm3,FIX+FPASS+3);    pIn[6*8] = (TYPE)SHIFTR(mm2,FIX+FPASS+3);    mm0 += HALF(FPASS+3) - 1;    BUTF(mm0, mm1, Spill);    pIn[0*8] = (TYPE)SHIFTR(mm0, FPASS+3);    pIn[4*8] = (TYPE)SHIFTR(mm1, FPASS+3);      // even    mm3 = mm5 + mm7;    mm2 = mm4 + mm6;    ROTATE(mm2, mm3,  ROT17_C, -ROT17_SpC, -ROT17_SmC, mm0, 0, HALF(FIX+FPASS+3));    ROTATE2(mm4, mm7, -ROT37_C,  ROT37_SpC,  ROT37_SmC, mm0);    mm7 += mm3;    mm4 += mm2;    pIn[7*8] = (TYPE)SHIFTR(mm4,FIX+FPASS+3);    pIn[1*8] = (TYPE)SHIFTR(mm7,FIX+FPASS+3);    ROTATE2(mm5, mm6, -ROT13_C,  ROT13_SmC,  ROT13_SpC, mm0);    mm5 += mm3;    mm6 += mm2;    pIn[5*8] = (TYPE)SHIFTR(mm5,FIX+FPASS+3);    pIn[3*8] = (TYPE)SHIFTR(mm6,FIX+FPASS+3);    pIn++;  }}//////////////////////////////////////////////////////////void Skl_IDct16_C( TYPE *In ){  TYPE *pIn;  int i;  pIn = In;  for (i=8; i>0; --i)  {    int mm0, mm1, mm2, mm3, mm4, mm5, mm6, mm7, Spill;      // odd        mm4 = (int)pIn[7];    mm5 = (int)pIn[5];    mm6 = (int)pIn[3];    mm7 = (int)pIn[1];    mm2 = mm4 + mm6;    mm3 = mm5 + mm7;    ROTATE2(mm2, mm3,  ROT17_C, -ROT17_SpC, -ROT17_SmC, mm1);    ROTATE2(mm4, mm7, -ROT37_C,  ROT37_SpC,  ROT37_SmC, mm1);    ROTATE2(mm5, mm6, -ROT13_C,  ROT13_SmC,  ROT13_SpC, mm1);    mm4 += mm2;    mm5 += mm3;    mm6 += mm2;    mm7 += mm3;      // even    mm3 = (int)pIn[2];    mm2 = (int)pIn[6];    ROTATE2(mm3, mm2, ROT6_C, ROT6_SmC, -ROT6_SpC, mm1);    LOAD_BUTF(mm0, mm1, 0, 4, Spill, pIn);    mm0 = SHIFTL(mm0, FIX) + HALF(FIX-IPASS);    mm1 = SHIFTL(mm1, FIX) + HALF(FIX-IPASS);    BUTF(mm0, mm3, Spill);    BUTF(mm1, mm2, Spill);    BUTF(mm0, mm7, Spill);    pIn[0] = SHIFTR(mm0, FIX-IPASS);    pIn[7] = SHIFTR(mm7, FIX-IPASS);    BUTF(mm1, mm6, mm0);    pIn[1] = SHIFTR(mm1, FIX-IPASS);    pIn[6] = SHIFTR(mm6, FIX-IPASS);    BUTF(mm2, mm5, mm0);    pIn[2] = SHIFTR(mm2, FIX-IPASS);    pIn[5] = SHIFTR(mm5, FIX-IPASS);    BUTF(mm3, mm4, mm0);    pIn[3] = SHIFTR(mm3, FIX-IPASS);    pIn[4] = SHIFTR(mm4, FIX-IPASS);    pIn += 8;  }  pIn = In;  for (i=8; i>0; --i)  {    int mm0, mm1, mm2, mm3, mm4, mm5, mm6, mm7, Spill;      // odd        mm4 = (int)pIn[7*8];    mm5 = (int)pIn[5*8];    mm6 = (int)pIn[3*8];    mm7 = (int)pIn[1*8];    mm2 = mm4 + mm6;    mm3 = mm5 + mm7;    ROTATE2(mm2, mm3,  ROT17_C, -ROT17_SpC, -ROT17_SmC, mm1);    ROTATE2(mm4, mm7, -ROT37_C,  ROT37_SpC,  ROT37_SmC, mm1);    ROTATE2(mm5, mm6, -ROT13_C,  ROT13_SmC,  ROT13_SpC, mm1);    mm4 += mm2;    mm5 += mm3;    mm6 += mm2;    mm7 += mm3;          // even    mm3 = (int)pIn[2*8];    mm2 = (int)pIn[6*8];    ROTATE2(mm3, mm2, ROT6_C, ROT6_SmC, -ROT6_SpC, mm1);    LOAD_BUTF(mm0, mm1, 0*8, 4*8, Spill, pIn);    mm0 = SHIFTL(mm0, FIX) + HALF(FIX+IPASS+3);    mm1 = SHIFTL(mm1, FIX) + HALF(FIX+IPASS+3);    BUTF(mm0, mm3, Spill);    BUTF(mm1, mm2, Spill);    BUTF(mm0, mm7, Spill);    pIn[8*0] = (TYPE) SHIFTR(mm0, FIX+IPASS+3);    pIn[8*7] = (TYPE) SHIFTR(mm7, FIX+IPASS+3);    BUTF(mm1, mm6, mm0);    pIn[8*1] = (TYPE) SHIFTR(mm1, FIX+IPASS+3);    pIn[8*6] = (TYPE) SHIFTR(mm6, FIX+IPASS+3);    BUTF(mm2, mm5, mm0);    pIn[8*2] = (TYPE) SHIFTR(mm2, FIX+IPASS+3);    pIn[8*5] = (TYPE) SHIFTR(mm5, FIX+IPASS+3);    BUTF(mm3, mm4, mm0);    pIn[8*3] = (TYPE) SHIFTR(mm3, FIX+IPASS+3);    pIn[8*4] = (TYPE) SHIFTR(mm4, FIX+IPASS+3);    pIn++;  }}void Skl_IDct16_Sparse_C( TYPE *In ){