OFFSET
1,2
LINKS
Robert Israel, Table of n, a(n) for n = 1..10000
FORMULA
a(n) = A000005(n) + A004788(n-1). - Vladeta Jovovic, Apr 07 2009 (Corrected by Altug Alkan, Nov 25 2015)
EXAMPLE
For n = 8 we have the following proper divisors for k <= n: {1}, {1}, {1}, {1, 2}, {1}, {1, 2, 3}, {1}, {1, 2, 4}. Only k = 6 has a proper divisor that is not a divisor of 8, viz. 3. Hence a(8) = 7.
MAPLE
N:= 1000: # to get a(1) to a(N)
A:= Vector(N, numtheory:-tau):
for p in select(isprime, [2, seq(i, i=3..N, 2)]) do
for d from 0 to floor(log[p](N))-1 do
R:= [seq(seq(p^d*(i+p*j), j=1..floor((N/p^d - i)/p)), i=1..p-1)];
A[R]:= map(`+`, A[R], 1);
od
od:
convert(A, list); # Robert Israel, Nov 24 2015
MATHEMATICA
f[n_] := Block[{d = Most@ Divisors@ n}, Select[Range@ n, Union[Most@ Divisors@ #, d] == d &]]; Array[Length@ f@ # &, {75}] (* Michael De Vlieger, Nov 25 2015 *)
PROG
(Magma) [ #[ k: k in [1..n] | forall(t){ d: d in Divisors(k) | d eq k or d in Divisors(n) } ]: n in [1..75] ];
(PARI) a004788(n) = {sfp = Set(); for (k=1, n-1, sfp = setunion(sfp, Set(factor(binomial(n, k))[, 1]))); return (length(sfp)); }
a(n) = numdiv(n) + a004788(n-1); \\ Altug Alkan, Nov 25 2015
CROSSREFS
KEYWORD
nonn
AUTHOR
Jaroslav Krizek, Apr 01 2009
EXTENSIONS
Edited and extended by Klaus Brockhaus, Apr 06 2009
STATUS
approved