\\ 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, k, [z, z=prd(z, 1+I)][1])

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

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);
	);
	imag(z)
}

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

quit