Year-end appeal: Please make a donation to the OEIS Foundation to support ongoing development and maintenance of the OEIS. We are now in our 61st year, we have over 378,000 sequences, and we’ve reached 11,000 citations (which often say “discovered thanks to the OEIS”).
%I #64 Nov 17 2024 18:55:50
%S 12,1222,1437286,3441373,1032893366969
%N Fixed points of A341885.
%C Numbers n such that A341885(n) = n.
%C Includes 2*p*q if p and q are primes such that p^2-4*p*q+q^2+p+q+6 = 0. This includes 12 for p=2, q=3, 1222 for p=13,q=47, 1437286 for p=439, q=1637, and 76498942675946443126 for p=3201392659, q=11947760057.
%C Another term: 6538810199342921107066977217325653068509 = 13 * 4401624135264074597*114272683103433355069. - _Chai Wah Wu_, Feb 25 2021
%F A341885(a(n)) = a(n).
%e a(2) = 1222 is a term because 1222 = 2*13*47 and A341885(1222) = 2*3/2 + 13*14/2 + 47*48/2 = 1222.
%p f:= proc(n) local F,t;
%p F:= ifactors(n)[2];
%p add(t[1]*(t[1]+1)/2*t[2],t=F)
%p end proc:
%p select(t -> f(t)=t, [$1..4000000]);
%t Block[{a = {}}, Monitor[Do[If[# == i, AppendTo[a, i]] &@ Total[PolygonalNumber@ Flatten[ConstantArray[#1, #2] & @@@ FactorInteger[i]]], {i, 2, 4*10^6}], i]; a] (* _Michael De Vlieger_, Feb 22 2021 *)
%o (Python)
%o from sympy import factorint
%o A340834_list = [n for n in range(2,10**4) if n == sum(k*m*(m+1)//2 for m,k in factorint(n).items())] # _Chai Wah Wu_, Feb 25 2021
%Y Cf. A341885.
%K nonn,more
%O 1,1
%A _J. M. Bergot_ and _Robert Israel_, Feb 22 2021
%E a(5) from _Martin Ehrenstein_, Mar 07 2021