%I #21 Apr 19 2020 23:40:19
%S 1,-1,1,0,-1,0,1,0,-1,0,5,0,-691,0,35,0,-3617,0,43867,0,-1222277,0,
%T 854513,0,-236364091,0,8553103,0,-23749461029,0,8615841276005,0,
%U -84802531453387,0,90219075042845,0,-26315271553053477373,0,38089920879940267
%N a(n) = BS(n) * W(n) where BS = Sum_{k=0..n} ((-1)^k*k!/(k+1)) S(n, k) and S(n, k) the Stirling subset numbers A048993(n, k). W(n) = Product_{ p primes <= n+1 such that p divides n+1 or p-1 divides n } = A225481(n).
%C a(n)/A225481(n) is a representation of the Bernoulli numbers. This is case m = 1 of the scaled generalized Bernoulli numbers defined as Sum_{k=0..n} ((-1)^k*k!/(k+1)) S_{m}(n,k) where S_{m}(n,k) are generalized Stirling subset numbers. A225481(n) can be seen as an analog of the Clausen numbers A141056(n). Reduced to lowest terms a(n)/A225481(n) becomes A027641(n)/A027642(n).
%H Peter Luschny, <a href="http://www.luschny.de/math/euler/StirlingFrobeniusNumbers.html">Stirling-Frobenius numbers</a>
%H Peter Luschny, <a href="http://www.luschny.de/math/euler/GeneralizedBernoulliNumbers.html">Generalized Bernoulli numbers</a>.
%e The numerators of 1/1, -1/2, 1/6, 0/2, -1/30, 0/6, 1/42, 0/2, -1/30, 0/10, 5/66, 0/6, -691/2730, 0/14, 35/30, 0/2, -3617/510, 0/6, 43867/798, ... (the denominators are A225481(n)).
%t BS[n_] := Sum[((-1)^k*k!/(k + 1)) StirlingS2[n, k], {k, 0, n}];
%t W[n_] := Product[If[Divisible[n + 1, p] || Divisible[n, p - 1], p, 1], {p, Prime /@ Range[PrimePi[n + 1]]}];
%t a[n_] := BS[n] W[n];
%t Table[a[n], {n, 0, 38}] (* _Jean-François Alcover_, Jul 08 2019 *)
%o (Sage)
%o @CachedFunction
%o def EulerianNumber(n, k, m) : # -- The Eulerian numbers --
%o if n == 0: return 1 if k == 0 else 0
%o return (m*(n-k)+m-1)*EulerianNumber(n-1, k-1, m) + \
%o (m*k+1)*EulerianNumber(n-1, k, m)
%o @CachedFunction
%o def SF_BS(n, m): # -- The scaled Stirling-Frobenius Bernoulli numbers --
%o return add(add(EulerianNumber(n, j, m)*binomial(j, n - k) \
%o for j in (0..n))/((-m)^k*(k+1)) for k in (0..n))
%o def A226156(n): # -- The numerators of SF_BS(n, 1) relative to A225481 --
%o C = mul(filter(lambda p: ((n+1)%p == 0) or (n%(p-1) == 0), primes(n+2)))
%o return C*SF_BS(n, 1)
%o [A226156(n) for n in (0..25)]
%Y Cf. A225481, A226157.
%K sign,frac
%O 0,11
%A _Peter Luschny_, May 30 2013