\\ product of Eisenstein integers
prd(u,v) = {
	(real(u)*real(v) - imag(u)*imag(v)) + (imag(u)*real(v)+real(u)*imag(v)-imag(u)*imag(v))*I;
}

\\ powers of 1+w
z = 1
p1w = vector(6, x, [z,z=prd(z,1+I)][1])

g(d) = {
	if (d==0,
			0,
			my (u=(d-1)\2, v=(d-1)%2);
			prd(p1w[1+u], if (v, 2+I, 1))
	);
}

a(n) = {
	my (d=Vecrev(digits(n, 13)), p=1, z=0);
	for (k=1, #d,
		z+=prd(g(d[k]), p);
		p=prd(p, 4+I);
	);
	real(z)
}

{
	for (n=0, 13^3-1,
		print (n " "a(n));
	);
}

quit