Evaluating the formula from:
http://www.johndcook.com/blog/2015/08/20/when-is-a-triangle-a-square/
from sympy import var, sqrt, init_printing
init_printing()
var("n")
e = ((17+12*sqrt(2))**n + (17-12*sqrt(2))**n - 2 ) / 32
e
for n_ in range(1, 20):
print n_, e.subs(n, n_).expand()
1 1 2 36 3 1225 4 41616 5 1413721 6 48024900 7 1631432881 8 55420693056 9 1882672131025 10 63955431761796 11 2172602007770041 12 73804512832419600 13 2507180834294496361 14 85170343853180456676 15 2893284510173841030625 16 98286503002057414584576 17 3338847817559778254844961 18 113422539294030403250144100 19 3853027488179473932250054441