This is a front-end for the Online Encyclopedia of Integer Sequences, made by Christian Perfect. The idea is to provide OEIS entries in non-ancient HTML, and then to think about how they're presented visually. The source code is on GitHub.
%I A056626 #37 Aug 04 2024 01:25:43
%S A056626 0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,1,0,0,0,0,2,0,0,
%T A056626 0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,1,0,1,0,0,0,0,0,0,0,2,0,0,0,0,
%U A056626 0,0,0,2,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,1,0,0,0,0,0,0,0,2,0,0,0,0,0,0,0,1,0
%N A056626 Number of non-unitary square divisors of n.
%H A056626 Michael De Vlieger, <a href="/A056626/b056626.txt">Table of n, a(n) for n = 1..10000</a>
%F A056626 a(n) = A046951(n) - 2^r(n), where r(n) is the number of distinct prime factors of the largest unitary square divisor of n. [Corrected by _Amiram Eldar_, Aug 03 2024]
%F A056626 a(n) = A046951(n) - 2^(A162641(n)). - _David A. Corneth_, Jul 28 2017
%F A056626 From _Amiram Eldar_, Sep 26 2022: (Start)
%F A056626 a(n) = A046951(n) - A056624(n).
%F A056626 Asymptotic mean: Limit_{m->oo} (1/m) * Sum_{k=1..m} a(k) = zeta(2)*(1 - 1/zeta(3)) = 0.27650128922802056073... . (End)
%e A056626 n = p^u prime power has u+1 square divisors of which 2 (i.e., 1 and n) are unitary but u-1 are not unitary, so a(p^u) = u - 1. E.g., n = 4^4 = 256, has 5 square divisors {1, 4, 16, 64, 256} of which {4, 16, 64} are not unitary, so a(256)=3.
%t A056626 Table[DivisorSum[n, 1 &, And[IntegerQ@ Sqrt@ #, ! CoprimeQ[#, n/#]] &], {n, 105}] (* _Michael De Vlieger_, Jul 28 2017 *)
%t A056626 f1[p_, e_] := 1 + Floor[e/2]; f2[p_, e_] := 2^(1 - Mod[e, 2]); a[1] = 0; a[n_] := Times @@ f1 @@@ (fct = FactorInteger[n])- Times @@ f2 @@@ fct; Array[a, 100] (* _Amiram Eldar_, Sep 26 2022 *)
%o A056626 (PARI) a(n) = {my(f = factor(n), r=0, m = 0); prod(i=1,#f~,f[i,2]>>1 + 1) - 2^(omega(f) - omega(core(f)))} \\ _David A. Corneth_, Jul 28 2017
%o A056626 (PARI) a(n) = sumdiv(n, d, if(gcd(d, n/d)!=1, issquare(d))); \\ _Michel Marcus_, Jul 29 2017
%o A056626 (Python)
%o A056626 from math import prod
%o A056626 from sympy import factorint
%o A056626 def A056626(n):
%o A056626 f = factorint(n).values()
%o A056626 return prod((e>>1)+1 for e in f)-(1<<sum(e&1^1 for e in f)) # _Chai Wah Wu_, Aug 04 2024
%Y A056626 Cf. A000188, A008833, A034444, A046951, A055229, A056624, A162641.
%K A056626 nonn
%O A056626 1,32
%A A056626 _Labos Elemer_, Aug 08 2000
%E A056626 a(32) and a(96) corrected by _Michael De Vlieger_, Jul 29 2017