%I #23 Feb 10 2019 23:02:35
%S 1,5,39,369,3898,44239,528083,6544745,83496720,1090091650,14501708246,
%T 195954553755,2682953977174,37150480629539,519455719162283,
%U 7325383709872345,104080732316126716,1488685017986884528,21420051312840487968
%N a(n) = Sum_{k=1..n} binomial(n+k-1, n)^2 / n.
%C Note that Sum_{k=1..n} binomial(n+k-1, n) / n = Catalan(n) = A000108(n).
%C p divides a((p-1)/2) for prime p = {5, 13, 17, 29, 37, 41, 53, ...) = A002144 Pythagorean primes: primes of form 4n + 1. - _Alexander Adamchuk_, Dec 27 2013
%F G.f.: has an anti-derivative of a hypergeometric function, see Maple program. - _Mark van Hoeij_, May 05 2013
%F Recurrence: 2*n^2*(2*n + 1)*(21*n^2 - 62*n + 46)*a(n) = (1365*n^5 - 6067*n^4 + 9948*n^3 - 7478*n^2 + 2640*n - 360)*a(n-1) - 4*(n-2)*(2*n - 3)^2*(21*n^2 - 20*n + 5)*a(n-2). - _Vaclav Kotesovec_, Mar 02 2014
%F a(n) ~ 16^n / (3*Pi*n^2). - _Vaclav Kotesovec_, Mar 02 2014
%p ogf := (4-x)^(1/2)*x^(-3/2)*Int((x+5/4)*hypergeom([1/2, 1/2],[1],16*x)/((4-x)^(3/2)*x^(1/2)),x) - 5/(8*x);
%p series(eval(ogf, Int = proc(a,x) int(series(a,x=0,30),x) end), x=0, 30); # _Mark van Hoeij_, May 05 2013
%t Table[ Sum[ Binomial[ n+k-1, n ]^2, {k,1,n} ] / n, {n,1,30} ]
%Y Cf. A000108 (Catalan numbers), A002144 (Pythagorean primes).
%K nonn
%O 1,2
%A _Alexander Adamchuk_, May 15 2007