A351130 - OEIS

%I #22 Apr 01 2022 23:47:11

%S 0,1,3,21,351,43956,187578,246753,570846,1200475,4890628,15671601,

%T 83663580,442903203,3776109156,35358717628,1060996913571,

%U 2443123072855,65801068940503,598914888327003,1364298098094561

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

%C Triangular numbers in A229511.

%e 21 = A000217(6), 21' = 10 = A000217(4), so 21 is a term.

%e 351 = A000217(26), 351' = 378 = A000217(27), so 351 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, 10^5}], IntegerQ[Sqrt[8*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..90000]]| tr(Floor(f(n)))];

%o (PARI) lista(nn) = my(t); for (n=0, nn, if (ispolygonal(der(t=n*(n+1)/2), 3), print1(t, ", "))); \\ _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(t)):

%o             yield t

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

%Y Intersection of A000217 and A229511.

%Y Cf. A003415.

%K nonn,more

%O 1,3

%A _Marius A. Burtea_, Feb 07 2022

%E a(20)-a(21) from _Michael S. Branicky_, Feb 17 2022