%I #6 Sep 01 2021 22:21:41
%S 1,1,-4,120,-1344,3024000,-1576143360,101708006400,-2522591034163200,
%T 6686974460694528000,-1287307431968882688000,
%U 160078872315904478576640000,-53718579665963356985229312000,574898901006059006921736192000000,-241364461951740682229320388129587200000
%N a(n) = Bernoulli(2*n) * (2*n+1)! if 2*n+1 is a prime, otherwise a(n) = Bernoulli(2*n) * (2*n)!.
%H <a href="/index/Be#Bernoulli">Index entries for sequences related to Bernoulli numbers</a>
%F a(n) is the numerator of Bernoulli(2*n) * (2*n)! (for denominators see A128059).
%F a(n) is the numerator of (2*n)!^2 * [x^(2*n)] x * coth(x/2) / 2.
%F a(n) is the numerator of b(2*n) where b(n) = -Sum_{k=1..n} binomial(n,k)^2 * k! * b(n-k) / (k+1), b(0) = 1.
%e Bernoulli(2*n) * (2*n)! = [ 1, 1/3, -4/5, 120/7, -1344, 3024000/11, -1576143360/13, 101708006400, -2522591034163200/17, 6686974460694528000/19, ... ].
%t a[n_] := If[PrimeQ[2 n + 1], BernoulliB[2 n] (2 n + 1)!, BernoulliB[2 n] (2 n)!]; Table[a[n], {n, 0, 14}]
%t Table[Numerator[BernoulliB[2 n] (2 n)!], {n, 0, 14}]
%t Table[Numerator[(2 n)!^2 SeriesCoefficient[x Coth[x/2]/2, {x, 0, 2 n}]], {n, 0, 14}]
%t b[0] = 1; b[n_] := b[n] = -Sum[Binomial[n, k]^2 k! b[n - k]/(k + 1), {k, 1, n}]; a[n_] := Numerator[b[2 n]]; Table[a[n], {n, 0, 14}]
%o (PARI) a(n) = numerator(bernfrac(2*n)*(2*n)!); \\ _Michel Marcus_, Sep 01 2021
%Y Cf. A000367, A001332, A002431, A002445, A036278, A128059.
%K sign,frac
%O 0,3
%A _Ilya Gutkovskiy_, Sep 01 2021