OFFSET
1,2
COMMENTS
Note that 1 is regarded as free of squares of primes and is also a square number and a unitary divisor.
LINKS
Antti Karttunen, Table of n, a(n) for n = 1..10000
Steven R. Finch, Unitarism and Infinitarism, February 25, 2004. [Cached copy, with permission of the author]
FORMULA
Multiplicative with a(p) = 2 and a(p^e) = 1 for e > 1. a(n) = 2^A056169(n). - Vladeta Jovovic, Nov 01 2001
From Vaclav Kotesovec, Feb 11 2023: (Start)
Dirichlet g.f.: zeta(s) * Product_{primes p} (1 + 1/p^s - 1/p^(2*s)).
Dirichlet g.f.: zeta(s)^2 * Product_{primes p} (1 - 2/p^(2*s) + 1/p^(3*s)), (with a product that converges for s=1).
Let f(s) = Product_{primes p} (1 - 2/p^(2*s) + 1/p^(3*s)), then Sum_{k=1..n} a(k) ~ n * (f(1) * (log(n) + 2*gamma - 1) + f'(1)), where f(1) = Product_{primes p} (1 - 2/p^2 + 1/p^3) = A065464 = 0.42824950567709444021876..., f'(1) = f(1) * Sum_{primes p} (4*p-3) * log(p) / (p^3 - 2*p + 1) = 0.808661108949590913395... and gamma is the Euler-Mascheroni constant A001620. (End)
a(n) = Sum_{d|n, gcd(d,n/d)=1} mu(d)^2. - Wesley Ivan Hurt, May 25 2023
a(n) = Sum_{d|n} A343443(d)*mu(n/d). - Ridouane Oudra, Dec 18 2023
EXAMPLE
n = 252 = 2*2*3*3*7 has 18 divisors, 8 unitary and 8 squarefree divisors of which 2 are unitary and squarefree, divisors {1,7};
n = 2520 = 2*2*2*3*3*5*7 has 48 divisors, 16 unitary and 16 squarefree divisors of which {1,5,7,35} are both, thus a(2520) = 4.
a(2520) = a(2^3*3^2*5*7) = a(2^3)*a(3^2)*a(5)*a(7) = 1*1*2*2 = 4.
MATHEMATICA
Array[DivisorSigma[0, #] &@ Denominator[#/Apply[Times, FactorInteger[#][[All, 1]]]^2] &, 105] (* or *)
Table[DivisorSum[n, 1 &, And[SquareFreeQ@ #, CoprimeQ[#, n/#]] &], {n, 105}] (* Michael De Vlieger, Jul 19 2017 *)
f[p_, e_] := If[e==1, 2, 1]; a[1] = 1; a[n_] := Times @@ (f @@@ FactorInteger[n]); Array[a, 100] (* Amiram Eldar, May 14 2019 *)
PROG
(PARI)
A057521(n) = { my(f=factor(n)); prod(i=1, #f~, if(f[i, 2]>1, f[i, 1]^f[i, 2], 1)); } \\ Charles R Greathouse IV, Aug 13 2013
\\ Or:
A055229(n) = { my(c=core(n)); gcd(c, n/c); }; \\ This function from Charles R Greathouse IV, Nov 20 2012
(PARI) for(n=1, 100, print1(direuler(p=2, n, (1 + X - X^2)/(1-X))[n], ", ")) \\ Vaclav Kotesovec, Feb 11 2023
(Scheme) (define (A056671 n) (if (= 1 n) n (* (if (= 1 (A067029 n)) 2 1) (A056671 (A028234 n))))) ;; (After the given multiplicative formula) - Antti Karttunen, Jul 19 2017
(Python)
from sympy import factorint, prod
def a(n): return 1 if n==1 else prod([2 if e==1 else 1 for p, e in factorint(n).items()])
print([a(n) for n in range(1, 51)]) # Indranil Ghosh, Jul 19 2017
CROSSREFS
KEYWORD
mult,nonn
AUTHOR
Labos Elemer, Aug 10 2000
STATUS
approved