%I #22 Aug 28 2023 09:41:33
%S 1,2,4,30,24,1210,18396,235998,4793456,76168850,1282320348,
%T 25100418046,481341997032,10452086347274,237925595533164,
%U 5524220670435982,136705837928870368,3444192369181374754,89772662325079950436,2431910317560215089758,67517711482300160612104
%N Number of ways of writing n^2 as a sum of n squares.
%H Alois P. Heinz, <a href="/A232173/b232173.txt">Table of n, a(n) for n = 0..200</a> (first 101 terms from Paul D. Hanna)
%F a(n) equals the coefficient of x^(n^2) in the n-th power of Jacobi theta_3(x) where theta_3(x) = 1 + 2*Sum_{n>=1} x^(n^2).
%e There are a(4) = 24 solutions (w,x,y,z) of 4^2 = w^2 + x^2 + y^2 + z^2:
%e (2,2,2,2), (-2,-2,-2,-2), 6 permutations of (2,2,-2,-2),
%e 4 permutations of (2,2,2,-2), 4 permutations of (2,-2,-2,-2),
%e 4 permutations of (4,0,0,0), and 4 permutations of (-4,0,0,0).
%e To illustrate a(n) = the coefficient of x^(n^2) in theta_3(x)^n, where
%e theta_3(x) = 1 + 2*x + 2*x^4 + 2*x^9 + 2*x^16 + 2*x^25 + 2*x^36 + 2*x^49 +...,
%e form a table of coefficients of x^k in theta_3(x)^n, n>=0, like so:
%e n\k:0..1...2...3...4...5...6...7...8...9..10..11..12..13..14..15..16....
%e 0:[(1),0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...];
%e 1: [1,(2), 0, 0, 2, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 2,...];
%e 2: [1, 4, 4, 0, (4), 8, 0, 0, 4, 4, 8, 0, 0, 8, 0, 0, 4,...];
%e 3: [1, 6, 12, 8, 6, 24, 24, 0, 12,(30),24, 24, 8, 24, 48, 0, 6,...];
%e 4: [1, 8, 24, 32, 24, 48, 96, 64, 24,104,144, 96, 96,112,192,192,(24),...];
%e 5: [1,10, 40, 80, 90,112,240,320,200,250,560,560,400,560,800,960,730,...];
%e then the coefficients in parenthesis form the initial terms of this sequence.
%p b:= proc(n, t) option remember; `if`(n=0, 1, `if`(n<0 or t<1, 0,
%p b(n, t-1) +2*add(b(n-j^2, t-1), j=1..isqrt(n))))
%p end:
%p a:= n-> b(n^2, n):
%p seq(a(n), n=0..20); # _Alois P. Heinz_, Mar 10 2023
%t b[n_, t_] := b[n, t] = If[n == 0, 1, If[n < 0 || t < 1, 0, b[n, t - 1] + 2*Sum[b[n - j^2, t - 1], {j, 1, Floor@Sqrt[n]}]]];
%t a[n_] := b[n^2, n];
%t Table[a[n], {n, 0, 20}] (* _Jean-François Alcover_, Aug 28 2023, after _Alois P. Heinz_ *)
%o (PARI) {a(n)=local(THETA3=1+2*sum(m=1, n+1, x^(m^2))+x*O(x^(n^2))); polcoeff(THETA3^n, n^2)}
%o for(n=0, 30, print1(a(n), ", "))
%Y Cf. A066535.
%Y Main diagonal of A302996.
%K nonn
%O 0,2
%A _Paul D. Hanna_, Nov 19 2013