A010846 Number of numbers <= n whose set of prime factors is a subset of the set of prime factors of n.
1, 2, 2, 3, 2, 5, 2, 4, 3, 6, 2, 8, 2, 6, 5, 5, 2, 10, 2, 8, 5, 7, 2, 11, 3, 7, 4, 8, 2, 18, 2, 6, 6, 8, 5, 14, 2, 8, 6, 11, 2, 19, 2, 9, 8, 8, 2, 15, 3, 12, 6, 9, 2, 16, 5, 11, 6, 8, 2, 26, 2, 8, 8, 7, 5, 22, 2, 10, 6, 20, 2, 18, 2, 9, 9, 10, 5, 23, 2, 14, 5, 9, 2, 28, 5, 9, 7, 11, 2, 32, 5, 10
Offset: 1
Examples
From _Wolfdieter Lang_, Jun 30 2014: (Start)
a(1) = 1 because the empty set is a subset of any set.
a(6) = 5 from the five numbers: 1 with the empty set, 2 with the set {2}, 3 with {3}, 4 with {2} and 6 with {2,3}, which are all subsets of {2,3}. 5 is out because {5} is not a subset of {2,3}. (End)
From _David A. Corneth_, Feb 10 2015: (Start)
Let p# be the product of primes up to p, A002110. Then,
a(13#) = 1161
a(17#) = 4843
a(19#) = 19985
a(23#) = 83074
a(29#) = 349670
a(31#) = 1456458
a(37#) = 6107257
a(41#) = 25547835
(End)
Links
- Michael De Vlieger, Table of n, a(n) for n = 1..10000 (first 5000 terms from T. D. Noe)
- Dorian Goldfeld, Modular forms, elliptic curves, and the ABC conjecture
Crossrefs
Cf. A162306 (numbers for each n).
Programs
-
Maple
A:= proc(n) local F, S, s,j,p; F:= numtheory:-factorset(n); S:= {1}; for p in F do S:= {seq(seq(s*p^j, j=0..floor(log[p](n/s))),s=S)} od; nops(S) end proc; seq(A(n),n=1..1000); # Robert Israel, Jun 27 2014 -
Mathematica
pf[n_] := If[n==1, {}, Transpose[FactorInteger[n]][[1]]]; SubsetQ[lst1_, lst2_] := Intersection[lst1,lst2]==lst1; Table[pfn=pf[n]; Length[Select[Range[n], SubsetQ[pf[ # ],pfn] &]], {n,100}] (* T. D. Noe, Jun 30 2009 *) Table[Total[MoebiusMu[#] Floor[n/#] &@ Select[Range@ n, CoprimeQ[#, n] &]], {n, 92}] (* Michael De Vlieger, May 08 2016 *) -
PARI
a(n,f=factor(n)[,1])=if(#f>1,my(v=f[1..#f-1],p=f[#f],s); while(n>0, s+=a(n,v); n\=p); s, if(#f&&n>0,log(n+.5)\log(f[1])+1,n>0)) \\ Charles R Greathouse IV, Jun 27 2013
-
PARI
a(n) = sum(k=1,n,if(gcd(n,k)-1,0,moebius(k)*(n\k))) \\ Benoit Cloitre, May 07 2016
-
PARI
a(n,f=factor(n)[,1])=if(#f<2, return(if(#f, valuation(n,f[1])+1, 0))); my(v=f[1..#f-1],p=f[#f],s); while(n, s+=a(n,v); n\=p); s \\ Charles R Greathouse IV, Nov 03 2021
-
Python
def A010846(n): return sum((m:=n**k)//k-(m-1)//k for k in range(1,n+1)) # Chai Wah Wu, Aug 15 2024
-
Python
from math import gcd from sympy import mobius def A010846(n): return sum(mobius(k)*(n//k) for k in range(1,n+1) if gcd(n,k)==1) # Chai Wah Wu, Apr 23 2025
Formula
a(n) = |{k<=n, k|n^(tau(k)-1)}|. - Vladeta Jovovic, Sep 13 2006
a(n) = Sum_{j = 1..n} Product_{primes p | j} delta(n mod p,0) where delta is the Kronecker delta. - Robert Israel, Feb 09 2015
a(n) = Sum_{1<=k<=n,(n,k)=1} mu(k)*floor(n/k). - Benoit Cloitre, May 07 2016
a(n) = Sum_{k=1..n} floor(n^k/k)-floor((n^k -1)/k). - Anthony Browne, May 28 2016
Extensions
Definition made more precise at the suggestion of Wolfdieter Lang
Comments