%I #10 Mar 31 2020 12:41:57
%S 28,86,628,2058,9310,25298,73220,126168,357238,423828,882418,1132550,
%T 1954860,2371648,2600598,3968188,4627280,6585390,7501858,10156328,
%U 14088548,24754940,26936208,32941678,47503218,61839490,72120200
%N If p and q are (odd) twin primes and q > p then p*q^2+(p+q)+1 is divisible by 3; a(n) = (p*q^2+(p+q)+1)/3.
%F a(n) = 2*A151990(n). - _R. J. Mathar_, Sep 18 2009
%p A001359 := proc(n) if n = 1 then 3; else for p from procname(n-1)+2 by 2 do if isprime(p) and isprime(p+2) then RETURN(p) ; fi; od: fi; end: A164689 := proc(n) p := A001359(n) ; (p+1)*(p^2+3*p+3)/3 ; end: seq(A164689(n),n=1..80) ; # _R. J. Mathar_, Sep 18 2009
%t (* b = A001359 *)
%t b[n_] := b[n] = If[n == 1, 3, Module[{p = NextPrime[b[n - 1]]}, While[ !PrimeQ[p + 2], p = NextPrime[p]]; p]];
%t a[n_] := With[{p = b[n]}, (p + 1)(p^2 + 3 p + 3)/3];
%t Array[a, 27] (* _Jean-François Alcover_, Mar 31 2020 *)
%Y Cf. A151990.
%K nonn
%O 1,1
%A Tanin (Mirza Sabbir Hossain Beg) (mirzasabbirhossainbeg(AT)yahoo.com), Aug 22 2009
%E Incorrect leading term deleted by _N. J. A. Sloane_, Sep 14 2009
%E More terms from _R. J. Mathar_, Sep 18 2009