login
a(n) = product{ p prime such that p divides n + 1 and p - 1 does not divide n }.
4

%I #16 Feb 24 2020 15:38:04

%S 1,1,1,1,1,3,1,1,1,5,1,3,1,7,5,1,1,3,1,5,7,11,1,3,1,13,1,7,1,15,1,1,

%T 11,17,35,3,1,19,13,5,1,21,1,11,1,23,1,3,1,5,17,13,1,3,55,7,19,29,1,

%U 15,1,31,7,1,13,33,1,17,23,35,1,3,1,37,5,19,77,39

%N a(n) = product{ p prime such that p divides n + 1 and p - 1 does not divide n }.

%H Charles R Greathouse IV, <a href="/A226040/b226040.txt">Table of n, a(n) for n = 0..10000</a>

%H Peter Luschny, <a href="http://www.luschny.de/math/euler/GeneralizedBernoulliNumbers.html">Generalized Bernoulli numbers</a>.

%F a(n) = A225481(n) / A141056(n).

%e a(41) = 21 = 3*7 = product({2,3,7} setminus {2}).

%p s:= (p, n) -> ((n+1) mod p = 0) and (n mod (p-1) <> 0);

%p A226040 := n -> mul(z, z = select(p->s(p,n), select('isprime', [$2..n])));

%p seq(A226040(n), n=0..77);

%t a[n_] := Times @@ Select[ FactorInteger[n+1][[All, 1]], !Divisible[n, #-1] &]; a[0] = 1; Table[a[n], {n, 0, 77}] (* _Jean-François Alcover_, Jun 27 2013, after Maple *)

%o (Sage)

%o def A226040(n):

%o F = filter(lambda p: ((n+1) % p == 0) and (n % (p-1)), primes(n))

%o return mul(F)

%o [A226040(n) for n in (0..77)]

%o (PARI) a(n)=my(f=factor(n+1)[,1],s=1);prod(i=1,#f,if(n%(f[i]-1),f[i],1)) \\ _Charles R Greathouse IV_, Jun 27 2013

%Y Cf. A160014, A226038, A226039, A225481, A141056.

%K nonn

%O 0,6

%A _Peter Luschny_, May 26 2013