%I #28 Dec 15 2022 06:30:08
%S 1,5,10,25,26,50,50,125,100,130,122,250,170,250,260,625,290,500,362,
%T 650,500,610,530,1250,676,850,1000,1250,842,1300,962,3125,1220,1450,
%U 1300,2500,1370,1810,1700,3250,1682,2500,1850,3050,2600,2650,2210,6250,2500
%N a(n) = A289310(n)^2 + A289311(n)^2.
%C This sequence is totally multiplicative.
%C a(n) > n^2 for any n > 1.
%C If n is a square, then a(n) is a square.
%C If a(n) and a(m) are squares, then a(n*m) is a square.
%C a(n) is also a square for nonsquares n = 42, 168, 246, 287, 378, 672, 984, 1050, 1148, 1434, 1512, 1673, 2058, 2214, 2583, 2688, ...
%H Rémy Sigrist, <a href="/A289320/b289320.txt">Table of n, a(n) for n = 1..10000</a>
%F Totally multiplicative, with a(p^k) = (1 + p^2)^k for any prime p and k > 0.
%F Sum_{k=1..n} a(k) ~ c * n^3, where c = 2/(Pi^2 * Product_{p prime} (1 - 1/p^2 - 1/p^3 - 1/p^4)) = 0.4778963213... . - _Amiram Eldar_, Nov 13 2022
%F Sum_{n>=1} 1/a(n) = 15/Pi^2 (A082020). - _Amiram Eldar_, Dec 15 2022
%t f[p_, e_] := (p^2 + 1)^e; a[1] = 1; a[n_] := Times @@ f @@@ FactorInteger[n]; Array[a, 50] (* _Amiram Eldar_, Nov 13 2022 *)
%o (PARI) a(n) = my (f=factor(n)); return (prod(i=1, #f~, (1 + f[i,1]^2) ^ f[i,2]))
%o (Python)
%o from sympy import factorint
%o from operator import mul
%o from functools import reduce
%o def a(n): return 1 if n==1 else reduce(mul, [(1 + p**2)**k for p, k in factorint(n).items()])
%o print([a(n) for n in range(1, 101)]) # _Indranil Ghosh_, Aug 03 2017
%Y Cf. A066872, A082020, A289310, A289311.
%K nonn,mult
%O 1,2
%A _Rémy Sigrist_, Jul 02 2017