%I #23 Sep 08 2022 08:46:12
%S 3,13,21,57,43,157,73,241,183,343,157,813,211,601,601,993,343,1561,
%T 421,1807,1057,1333,601,3661,993,1807,1641,3193,931,5257,1057,4033,
%U 2353,2971,2353,8373,1483,3661,3193,8191,1807,9313,1981,7141,6163,5257,2353
%N a(n) = 1 + sigma(n) + sigma(n)^2.
%H Robert Price, <a href="/A258774/b258774.txt">Table of n, a(n) for n = 1..10000</a>
%H OEIS Wiki, <a href="https://oeis.org/wiki/Cyclotomic Polynomials at x=n, n! and sigma(n)">Cyclotomic Polynomials at x=n, n! and sigma(n)</a>
%F a(n) = 1 + A000203(n) + A000203(n)^2.
%F a(n) = 1 + A000203(n) + A072861(n). - _Omar E. Pol_, Jun 19 2015
%p with(numtheory): A258774:=n->1+sigma(n)+sigma(n)^2: seq(A258774(n), n=1..100); # _Wesley Ivan Hurt_, Jul 09 2015
%t Table[1 + DivisorSigma[1, n] + DivisorSigma[1, n]^2, {n, 10000}]
%t Table[Cyclotomic[3, DivisorSigma[1, n]], {n, 10000}]
%o (Magma) [1+SumOfDivisors(n)+ SumOfDivisors(n)^2: n in [1..50]]; // _Vincenzo Librandi_, Jun 10 2015
%o (PARI) a(n)=my(s=sigma(n)); s^2+s+1 \\ _Charles R Greathouse IV_, Jun 10 2015
%o (Python)
%o from sympy import divisor_sigma
%o def A258774(n):
%o ....return (lambda x: x*(x+1)+1)(divisor_sigma(n)) # _Chai Wah Wu_, Jun 10 2015
%Y Cf. A000203 (sum of divisors of n).
%Y Cf. A258775 (indices of primes in this sequence), A258776 (corresponding primes).
%K easy,nonn
%O 1,1
%A _Robert Price_, Jun 09 2015