login
Number of cubefree integers not exceeding 10^n.
5

%I #29 Aug 08 2024 11:12:14

%S 1,9,85,833,8319,83190,831910,8319081,83190727,831907372,8319073719,

%T 83190737244,831907372522,8319073725828,83190737258105,

%U 831907372580692,8319073725807178,83190737258070643,831907372580707771

%N Number of cubefree integers not exceeding 10^n.

%C An alternate definition specifying "less than 10^n" would yield the same sequence except for the first 3 terms: 0, 8, 84, 833, 8319, etc. (since powers of 10 beyond 1000 are not cubefree anyhow).

%C The limit of a(n)/10^n is the inverse of Apery's constant, 1/zeta(3), whose digits are given by A088453.

%H Chai Wah Wu, <a href="/A160112/b160112.txt">Table of n, a(n) for n = 0..31</a> (terms 0..29 from Gerard P. Michon)

%H Gerard P. Michon, <a href="http://www.numericana.com/answer/counting.htm#cubefree">On the number of cubefree integers not exceeding N</a>.

%F a(n) = Sum_{i=1..floor(10^(n/3))} A008683(i)*floor(10^n/i^3).

%e a(0)=1 because 1 <= 10^0 is not a multiple of the cube of a prime.

%e a(1)=9 because the 9 numbers 1,2,3,4,5,6,7,9,10 are cubefree; 8 is not.

%e a(2)=85 because there are 85 cubefree integers equal to 100 or less.

%e a(3)=833 because there are 833 cubefree integers below 1000 (which is not cubefree itself).

%p with(numtheory): A160112:=n->add(mobius(i)*floor(10^n/(i^3)), i=1..10^n): (1,9,85,seq(A160112(n), n=3..5)); # _Wesley Ivan Hurt_, Aug 01 2015

%t Table[ Sum[ MoebiusMu[x]*Floor[10^n/(x^3)], {x, 10^(n/3)}], {n, 0, 18}] (* _Robert G. Wilson v_, May 27 2009 *)

%o (Python)

%o from sympy import mobius, integer_nthroot

%o def A160112(n): return sum(mobius(k)*(10**n//k**3) for k in range(1, integer_nthroot(10**n,3)[0]+1)) # _Chai Wah Wu_, Aug 06 2024

%o (Python)

%o from bitarray import bitarray

%o from sympy import integer_nthroot

%o def A160112(n): # faster program

%o q = 10**n

%o m = integer_nthroot(q,3)[0]+1

%o a, b = bitarray(m), bitarray(m)

%o a[1], p, i, c = 1, 2, 4, q-sum(q//k**3 for k in range(2,m))

%o while i < m:

%o j = 2

%o while i < m:

%o if j==p:

%o c -= (b[i]^1 if a[i] else -1)*(q//i**3)

%o j, a[i], b[i] = 0, 1, 1

%o else:

%o t1, t2 = a[i], b[i]

%o if (t1&t2)^1:

%o a[i], b[i] = (t1^1)&t2, ((t1^1)&t2)^1

%o c += (t2 if t1 else 2)*(q//i**3) if (t1^1)&t2 else (t2-2 if t1 else 0)*(q//i**3)

%o i += p

%o j += 1

%o p += 1

%o while a[p]|b[p]:

%o p += 1

%o i = p<<1

%o return c # _Chai Wah Wu_, Aug 06 2024

%Y A004709 (cubefree numbers). A088453 (limit of the string of digits). A160113 (binary counterpart for cubefree integers). A071172 & A053462 (decimal counterpart for squarefree integers). A143658 (binary counterpart for squarefree integers).

%K easy,nice,nonn

%O 0,2

%A _Gerard P. Michon_, May 02 2009, May 06 2009