View Single Post
Old 2020-02-14, 14:51   #7
Dr Sardonicus
 
Dr Sardonicus's Avatar
 
Feb 2017
Nowhere

3×5×239 Posts
Default

Note: The OP apparently forgot 1.

I worked the answer out for n = 10000.

In the first routine, I hand-calculated the limit on primes. The general bound is floor(log(n)/log(2)). For n = 10000 this is 13.

Code:
? n=10000
%5 = 10000

? j=0;for(i=1,n,if(ispower(i,2)||ispower(i,3)||ispower(i,5)||ispower(i,7)||ispower(i,11)||ispower(i,13),j++));print(j)
125
Of course, this takes a while. The following is much faster:

Code:
? s=0;for(i=2,n,if(issquarefree(i),t=-(floor(n^(1/i))-1)*moebius(i);s+=t;print(i" "s" "t));if(t==0,break));s++
2 99 99
3 119 20
5 124 5
6 121 -3
7 123 2
10 122 -1
11 123 1
13 124 1
14 124 0
%6 = 125
The print() statement can of course be eliminated. I merely put it in to show the "inclusion-exclusion" at work.

If n is very large, it is conceivable that there would be insufficient sig figs to calculate the roots to the nearest integer. My ancient version of Pari-GP has sqrtint(n), but doesn't have sqrtnint(n,i). Using sqrtnint(n,i) - 1 in place of floor(n^(1/i)) - 1 may avoid this problem.

Last fiddled with by Dr Sardonicus on 2020-02-14 at 14:57 Reason: xignif ostpy
Dr Sardonicus is offline   Reply With Quote