+static double c[8] = {1.0,
+0.980785280403230, // cos(Pi*1/16), sin(Pi*7/16)
+0.923879532511287, // cos(Pi*2/16), sin(Pi*6/16)
+0.831469612302545, // cos(Pi*3/16), sin(Pi*5/16)
+0.707106781186548, // cos(Pi*4/16), sin(Pi*4/16), 1/sqrt(2)
+0.555570233019602, // cos(Pi*5/16), sin(Pi*3/16)
+0.382683432365090, // cos(Pi*6/16), sin(Pi*2/16)
+0.195090322016128 // cos(Pi*7/16), sin(Pi*1/16)
+};
+
+static double cc[8];
+int ccquant = -1;
+
+static void preparequant(int quant)
+{
+ if(ccquant == quant)
+ return;
+ cc[0] = c[0]/(quant*2*4);
+ cc[1] = c[1]/(quant*2*4);
+ cc[2] = c[2]/(quant*2*4);
+ cc[3] = c[3]/(quant*2*4);
+ cc[4] = c[4]/(quant*2*4);
+ cc[5] = c[5]/(quant*2*4);
+ cc[6] = c[6]/(quant*2*4);
+ cc[7] = c[7]/(quant*2*4);
+ ccquant = quant;
+}
+
+inline static void innerdct(double*a,double*b, double*c)
+{
+ // c1*c7*2 = c6
+ // c2*c6*2 = c4
+ // c3*c5*2 = c2
+ // c4*c4*2 = 1
+
+ //{ 1, 3, 5, 7, -7, -5, -3, -1},
+ //{ 3, -7, -1, -5, 5, 1, 7, -3},
+ //{ 5, -1, 7, 3, -3, -7, 1, -5},
+ //{ 7, -5, 3, -1, 1, -3, 5, -7}
+ double b0,b1,b2,b3,b4,b5;
+ b2 = (a[0]+a[7]);
+ b3 = (a[1]+a[6]);
+ b4 = (a[2]+a[5]);
+ b5 = (a[3]+a[4]);
+
+ b0 = (b2+b5)*c[4];
+ b1 = (b3+b4)*c[4];
+ b[0*8] = b0 + b1;
+ b[4*8] = b0 - b1;
+ b[2*8] = (b2-b5)*c[2] + (b3-b4)*c[6];
+ b[6*8] = (b2-b5)*c[6] + (b4-b3)*c[2];
+
+ b0 = (a[0]-a[7]);
+ b1 = (a[1]-a[6]);
+ b2 = (a[2]-a[5]);
+ b3 = (a[3]-a[4]);
+
+ b[1*8] = b0*c[1] + b1*c[3] + b2*c[5] + b3*c[7];
+ b[3*8] = b0*c[3] - b1*c[7] - b2*c[1] - b3*c[5];
+ b[5*8] = b0*c[5] - b1*c[1] + b2*c[7] + b3*c[3];
+ b[7*8] = b0*c[7] - b1*c[5] + b2*c[3] - b3*c[1];
+}
+
+static void dct2(double*src, int*dest)
+{
+ double tmp[64], tmp2[64];
+ double*p;
+ int u,x,v,t;
+
+ for(v=0;v<8;v++)
+ {
+ double* a=&src[v*8];
+ double* b=&tmp[v];
+ innerdct(a,b,c);
+ }
+ for(v=0;v<8;v++)
+ {
+ double* a=&tmp[v*8];
+ double* b=&tmp2[v];
+ innerdct(a,b,cc);
+ }
+ for(t=0;t<64;t++) {
+ dest[zigzagtable[t]] = (int)(tmp2[t]);
+ }
+}
+
+