先を越された…
double c[8][8],s;int x,y,w[64];di(){for(x=64;x--;c[y=x/8][x%8]=cos((2*(x%8)+1)*
y*M_PI/16));}dc(int*b){for(x=64;x--;w[x]=s*(x/8?1:1/sqrt(2))*(x%8?1:1/sqrt(2))/
4)for(s=0,y=64;y--;)s+=b[y]*c[x%8][y%8]*c[x/8][y/8];}