View Single Post
 2020-02-14, 14:51 #7 Dr Sardonicus     Feb 2017 Nowhere 3×5×239 Posts 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