OFFSET
1,32
LINKS
Michael De Vlieger, Table of n, a(n) for n = 1..10000
FORMULA
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]
From Amiram Eldar, Sep 26 2022: (Start)
Asymptotic mean: Limit_{m->oo} (1/m) * Sum_{k=1..m} a(k) = zeta(2)*(1 - 1/zeta(3)) = 0.27650128922802056073... . (End)
EXAMPLE
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.
MATHEMATICA
Table[DivisorSum[n, 1 &, And[IntegerQ@ Sqrt@ #, ! CoprimeQ[#, n/#]] &], {n, 105}] (* Michael De Vlieger, Jul 28 2017 *)
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 *)
PROG
(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
(PARI) a(n) = sumdiv(n, d, if(gcd(d, n/d)!=1, issquare(d))); \\ Michel Marcus, Jul 29 2017
(Python)
from math import prod
from sympy import factorint
def A056626(n):
f = factorint(n).values()
return prod((e>>1)+1 for e in f)-(1<<sum(e&1^1 for e in f)) # Chai Wah Wu, Aug 04 2024
CROSSREFS
KEYWORD
nonn
AUTHOR
Labos Elemer, Aug 08 2000
EXTENSIONS
a(32) and a(96) corrected by Michael De Vlieger, Jul 29 2017
STATUS
approved