f1 = 1/(2 (x1^2 + y1^2))
f2 = 2/((x1 + Sqrt[3] x2)^2 + (y1 + Sqrt[3] y2)^2)
n = 100000
f = Log[Exp[n*f1] + Exp[n*f2]]/n
g = 0.5*Log[f]
c = -D[D[g, x1], x1] - D[D[g, y1], y1] - D[g, x2]^2 - D[g, y2]^2