%I #28 Mar 20 2023 09:02:54
%S 1,7,8,18,14,31,20,36,31,56,32,54,38,90,44,72,57,98,72,90,62,127,68,
%T 144,74,140,80,126,108,180,112,144,98,217,104,162,110,248,144,180,133,
%U 224,128,252,160,270,140,216,180,266,152,288,158,378,164,252,183,308
%N Sum of divisors of 3*n + 1.
%C Cubic AGM theta functions: a(q) (see A004016), b(q) (A005928), c(q) (A005882).
%H Amiram Eldar, <a href="/A144614/b144614.txt">Table of n, a(n) for n = 0..10000</a>
%F Expansion of q^(-1/3) * a(q) * c(q) / 3 in powers of q where a(), c() are cubic AGM theta functions. - _Michael Somos_, May 30 2012
%F From _Michael Somos_, Jun 09 2012: (Start)
%F Expansion of q^(-1/3) * (eta(q)^3 + 9 * eta(q^9)^3) * eta(q^3)^2 / eta(q) in powers of q.
%F a(n) = A000203(3*n + 1).
%F Convolution of A004016 and A033687. (End)
%F Sum_{k=1..n} a(k) = (2*Pi^2/9) * n^2 + O(n*log(n)). - _Amiram Eldar_, Dec 16 2022
%e G.f. = 1 + 7*x + 8*x^2 + 18*x^3 + 14*x^4 + 31*x^5 + 20*x^6 + 36*x^7 + 31*x^8 + 56*x^9 +...
%e G.f. = q + 7*q^4 + 8*q^7 + 18*q^10 + 14*q^13 + 31*q^16 + 20*q^19 + 36*q^22 + 31*q^25 + ...
%t a[ n_] := If[ n < 0, 0, DivisorSigma[1, 3 n + 1]]; (* _Michael Somos_, May 26 2014 *)
%t DivisorSigma[1,3*Range[0,60]+1] (* _Harvey P. Dale_, Mar 20 2023 *)
%o (PARI) {a(n) = if( n<0, 0, sigma( 3*n + 1))}; /* _Michael Somos_, May 30 2012 */
%o (PARI) {a(n) = my(A); if( n<0, 0, A = x * O(x^n); polcoeff( (eta(x + A)^3 + 9 * x * eta(x^9 + A)^3) * eta(x^3 + A)^2 / eta(x + A), n))}; /* _Michael Somos_, Jun 09 2012 */
%o (Sage) ModularForms( Gamma0(9), 2, prec=70).1; # _Michael Somos, May 26 2014 */
%o (Magma) Basis( ModularForms( Gamma0(9), 2), 173)[2]; /* _Michael Somos_, Jun 10 2015 */
%Y Cf. A000203, A004016, A005882, A005928, A016777, A033687.
%K nonn
%O 0,2
%A _N. J. A. Sloane_, Jan 15 2009