OFFSET
1,2
COMMENTS
a(n) is not a divisor of n for n = 21, 24, 42, 69, 84, 105, 115, 120, 138, 168, 171, ..., (A330736).
The sequence for the numerators only has terms for 1, 3, 5, 7, 11, 13, 17, 19, 23, 25, 29, 31, 33, ..., .
LINKS
Antti Karttunen, Table of n, a(n) for n = 1..10000
Antti Karttunen, Data supplement: n, a(n) computed for n = 1..65537
FORMULA
a(n) = n iff n is a power of a prime (A000961).
a(n) < n iff n is a member of A024619.
a(n) >= A034699(n). - Robert Israel, Aug 11 2019
gcd(a(n), n) = A330691(n). - Antti Karttunen, Dec 29 2019
MAPLE
H:= 1: B[1]:= 1:
for n from 2 to 200 do H:= H + 1/n; B[n]:= denom(H) od:
f:= proc(n) local F, t0, t;
t0:= max(seq(t[1]^t[2], t=ifactors(n)[2]));
for t from t0 do if B[t] mod n = 0 then return t fi od
end proc:
f(1):= 1:
map(f, [$1..100]); # Robert Israel, Aug 11 2019
MATHEMATICA
s = 0; k = 1; t[_] := 0; While[k < 101, s = s + 1/k; lst = Select[ Range@ 100, Mod[Denominator@ s, #] == 0 &]; If[t[#] == 0, t[#] = k] & /@ lst; k++]; t@# & /@ Range@75
PROG
(PARI) f(n) = denominator(sum(k=2, n, 1/k)); \\ A002805
a(n) = my(k=1); while(f(k) % n, k++); k; \\ Michel Marcus, Aug 11 2019
(PARI) A309639list(up_to) = { my(s=0, v002805=vector(up_to), v309639=vector(up_to)); v002805[1] = 1; for(k=2, up_to, s += 1/k; v002805[k] = denominator(s)); for(n=1, up_to, for(j=1, up_to, if(!(v002805[j]%n), v309639[n] = j; break))); (v309639); }; \\ Antti Karttunen, Dec 29 2019
CROSSREFS
KEYWORD
nonn
AUTHOR
Robert G. Wilson v, Aug 11 2019
STATUS
approved