OFFSET
1,2
COMMENTS
a(n) seems to be divisible by any positive integer as n approaches infinity. Example: There seem to be 6 terms without a 2 (i.e., 1, 3, 5, 9, 15, 105), 14 terms without a 3 (as large as 280), 22 terms without a factor of 4 (as large as 6930), and 29 terms without a 5 (as large as 3276). True for a(n) < 10^8.
For any term a(n) > 1, it seems that there exists at least one term a(x) such that a(x) * prime(y) = a(n) where prime(y) <= (7/2) * A053669(a(x)). True for a(n) < 10^8.
The ratio of prime(n)/a(n) for indices {1..8, 11, 15, 17} is almost 2 (with an error of at most 1 on the numerator). The exact ratios are 2/1, 3/2, 5/3, 7/4, 11/5, 13/6, 17/8, 19/9, 31/15, 47/24, 59/30. - Matthew Russell Downey, Jul 25 2023
LINKS
Matthew Russell Downey, On Bifurcations And Beauty (2023): pages 40-43.
FORMULA
Conjecture: a(n) ~ exp(4*(n-1)/21). - Matthew Russell Downey, Jul 25 2023
EXAMPLE
MATHEMATICA
f[p_, e_] := (e*(p - 1)/p + 1)*p^e; pil[n_] := Times @@ (f @@@ FactorInteger[n]); seq[nmax_] := Module[{p, pm = 0, s = {}}, Do[If[(p = pil[n]) > pm, pm = p; AppendTo[s, n]], {n, 1, nmax}]; s]; seq[1200] (* Amiram Eldar, Feb 13 2023 *)
PROG
(Python)
from sympy import factorint
from math import prod
def A018804(m): return prod(p**(e-1)*((p-1)*e+p) for p, e in factorint(m).items())
record = 0
for m in range(1, 2000):
value = A018804(m)
if value > record:
record = value
print(m, end=", ")
(PARI) f(n) = sumdiv(n, d, n*eulerphi(d)/d); \\ A018804
lista(nn) = my(r=0, list=List()); for (n=1, nn, my(m=f(n)); if (m > r, listput(list, n); r = m); ); Vec(list); \\ Michel Marcus, Feb 26 2023
CROSSREFS
KEYWORD
nonn
AUTHOR
Matthew Russell Downey, Feb 08 2023
STATUS
approved