sq[n_,q_] := sq[n,q] =
	1+Apply[Plus, Module[{r}, Table[ 2 q^(r^2), {r,1,n} ]
	]]

sq[n,q,0] = 1

sq[n_,q_,1] := sq[n,q]

sq[n_,q_,r_] := sq[n,q,r] =
	Normal[ Series[ sq[n,q,r-1] * sq[n,q], {q, 0, n^2}]]