A351131 - OEIS

%I #14 Apr 01 2022 23:47:21

%S 0,1,3,6,10,66,78,105,231,325,465,561,595,861,1378,2211,2278,2485,

%T 3081,3570,3655,4278,4465,5253,6441,6670,8515,8778,9453,9870,10011,

%U 10153,12561,13530,15051,18145,21115,21945,22578,23005,25878,27966,28441,40470,45753

%N Triangular numbers (A000217) whose second arithmetic derivative (A068346) is also a triangular number.

%e 6 = A000217(3), 6'' = 5' = 1 = A000217(1), so 6 is a term.

%e 66 = A000217(11), 66'' = 61' = 1 = A000217(1), so 66 is a term.

%e 325 = A000217(25), 325'' = 155' = 36 = A000217(8), so 325 is a term.

%t d[0] = d[1] = 0; d[n_] := n*Plus @@ ((Last[#]/First[#]) & /@ FactorInteger[n]); Select[Table[n*(n + 1)/2, {n, 0, 300}], IntegerQ[Sqrt[8*d[d[#]] + 1]] &] (* _Amiram Eldar_, Feb 07 2022 *)

%o (Magma) tr:=func<m|IsSquare(8*m+1)>;  f:=func<n |n le 1 select 0 else n*(&+[Factorisation(n)[i][2] / Factorisation(n)[i][1]: i in [1..#Factorisation(n)]])>; [n:n in [d*(d+1) div 2:d in [0..310]]| tr(Floor(f(Floor(f(n)))))];

%o (PARI) der(n) = my(f=factor(n)); vecsum([n/f[1]*f[2]|f<-factor(n+!n)~]); \\ A003415

%o isok(m) = ispolygonal(m, 3) && ispolygonal(der(der(m)), 3); \\ _Michel Marcus_, Feb 16 2022

%o (Python)

%o from itertools import count, islice

%o from sympy import factorint, integer_nthroot, isprime, nextprime

%o def istri(n): return integer_nthroot(8*n+1, 2)[1]

%o def ad(n):

%o     return 0 if n < 2 else sum(n*e//p for p, e in factorint(n).items())

%o def agen(): # generator of terms

%o     for i in count(0):

%o         t = i*(i+1)//2

%o         if istri(ad(ad(t))):

%o             yield t

%o print(list(islice(agen(), 45))) # _Michael S. Branicky_, Feb 16 2022

%Y Cf. A000217, A003415, A068346.

%K nonn

%O 1,3

%A _Marius A. Burtea_, Feb 07 2022