# encoding: utf-8

"""
by ulrik.sverdrup@gmail.com
lic: GPL3
"""

from __future__ import division
import scipy as sc
import scipy.linalg

sz = 8
base = sc.zeros([sz, sz, sz, sz])
v = sc.array(xrange(sz))

alpha = sc.ones([sz])* sc.sqrt(2/8)
alpha[0] = sc.sqrt(1/8)

# notice that sc.outer is the wrong one!
outer = sc.multiply.outer

cos = sc.cos( outer(v, sc.pi/8*(v + 1/2)) )
supercos = outer(cos,cos)
superalpha = outer(alpha, alpha)

for x in v:
	for y in v:
		base[x,y] = superalpha[x,y] * supercos[x,:,y,:]

superalpha = outer(alpha,alpha)
ncos = supercos
print superalpha
print supercos

print base



m = sc.array([
[52 , 55 , 61 , 66 , 70 , 61 , 64 , 73 ],
[63 , 59 , 55 , 90 , 109 , 85 , 69 , 72],
[62 , 59 , 68 , 113 , 144 , 104 , 66 , 73],
[63 , 58 , 71 , 122 , 154 , 106 , 70 , 69],
[67 , 61 , 68 , 104 , 126 , 88 , 68 , 70],
[79 , 65 , 60 , 70 , 77 , 68 , 58 , 75],
[85 , 71 , 64 , 59 , 55 , 61 , 65 , 83],
[87 , 79 , 69 , 68 , 65 , 76 , 78 , 94]])

m = m - 128
print m

def dct(m):
	res = sc.zeros([sz,sz])
	for x in v:
		for y in v:
			res[x,y] = sc.sum(sc.multiply(m,base[x,y]))
	return res

print sc.round_(dct(m))