Is there a power-counting function?

[enzocreti]

How many powers are there less than a given number x?
Is there a function that "counts" the powers less than x?


For example powers less than 10 are three (4,8,9)

[LaurV]

This sounds like a honest question.

There are:



(seriously, this would take milliseconds, for really huge numbers....)


Edit: actually, not. I was a bit stupid, as this would count the powers of powers again, for example the powers of 4 or 9 or 36 will be counted again, they were already counted as powers of 2, 3, 6 respectively. But this is easy to avoid.

[ATH]

Here is the OEIS sequence for it:

http://oeis.org/A070428

http://oeis.org/A070428/b070428.txt

[LaurV]

gp > cnt=0;for(i=1,100000,if(ispower(i),print(cnt++": "i)))
.....

<lot of lines removed>
....

359: 96721
360: 97336
361: 97344
362: 97969
363: 98596
364: 99225
365: 99856
366: 100000
gp > getnum(n)=s=0;for(i=2,log(n)/log(2),z=floor(n^(1/i)); s+=z-getnum(z)-1);s
%2 = (n)->s=0;for(i=2,log(n)/log(2),z=floor(n^(1/i));s+=z-getnum(z)-1);s
gp > getnum(10)
%3 = 3
gp > getnum(16)
%4 = 4
gp > getnum(10000)
%5 = 124
gp > ##
  ***   last result computed in 0 ms.
gp > getnum(100000)
%6 = 366
gp > ##
  ***   last result computed in 0 ms.
gp > getnum(10^50)
time = 15 ms.
%7 = 10000000046415898134516526
gp > ##
  ***   last result computed in 15 ms.
gp >

(it seems like my boss is not pushing me too hard today... )

[axn]

A faster version using inclusion/exclusion principle rather than recursion:

howmanyfactors(n)=#factor(n)[,1];
getnum(n)=my(s=0, z); for(i=2,log(n)/log(2), if(!issquarefree(i), next()); z=floor(n^(1/i))-1; if(howmanyfactors(i)%2==0, s-=z, s+=z)); s

EDIT:- Apparently I reinvented Möbius function

[sweety439]

Originally Posted by enzocreti

How many powers are there less than a given number x?
Is there a function that "counts" the powers less than x?


For example powers less than 10 are three (4,8,9)

powerpi(x)=floor(x^(1/2))+floor(x^(1/3))+floor(x^(1/5))+floor(x^(1/7))+floor(x^(1/11))+...

e.g. powerpi(9453)=97+21+6+3+2+2=131

[Dr Sardonicus]

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.

? 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:

? 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.

[LaurV]

Hihi, none of those are faster than the recursive version I posted, it takes longer than 15 ms to factor n in 10^50 range, as in my example. Well, there could be a "workaround" to always chose n extremely easy to factor (like 10^50 itself, which is trivial), because the result of the function won't change in a large range (around 10^50 there is no other "power" for quite a while! the precedent one is at 10^50-19999999999999999999999999)**, but for a general random n, mine is faster by far***
-------------
**challenge: the fastest "nextpower()" and "precpower()" functions?
***you may however, need to redefine realprecision to match the range in this case

[axn]

Originally Posted by LaurV

Hihi, none of those are faster than the recursive version I posted, it takes longer than 15 ms to factor n in 10^50 range, as in my example.

What are you talking about? No one is factoring n in the 10^50 range. The only numbers being factored are in the log2(n) range.

? \p 50
   realprecision = 57 significant digits (50 digits displayed)

getnumlaurv(n)=my(s=0, z); for(i=2,log(n)/log(2),z=floor(n^(1/i)); s+=z-getnumlaurv(z)-1); s
howmanyfactors(n)=#factor(n)[,1];
getnumaxn(n)=my(s=0, z); for(i=2,log(n)/log(2), if(!issquarefree(i), next()); z=floor(n^(1/i))-1; if(howmanyfactors(i)%2==0, s-=z, s+=z)); s
getnumdrs(n)=my(s=0, t); for(i=2,n,if(issquarefree(i),t=-(floor(n^(1/i))-1)*moebius(i);s+=t);if(t==0,break)); s++
? #
   timer = 1 (on)
? r=10^50+random(2^64)
%5 = 100000000000000000000000000000013282407956253574712
? 
? for(i=1, 1000, getnumlaurv(r))
time = 8,918 ms.
? for(i=1, 1000, getnumaxn(r))
time = 437 ms.
? for(i=1, 1000, getnumdrs(r))
time = 426 ms.
?

[LaurV]

Whoops...
fast reading, you confused me with those "n" and "i" mixture
Of course, you are right.
(I would however, test to see if pari is not caching the result, hihi, don't want to give up so easy, which is possible for your function, but not for mine, due to recursion - you know, like when reading files in windows, if you read the same file many times, subsequent accesses are much faster, hihi)
(just kidding)

[Dr Sardonicus]

My clunky Pari-GP script (including a print() statement for each term in the sum, and an inevitable superfluous zero term in the sum), my ancient version of Pari-GP, and doddering old computer still did the count for 10^50 in 2 milliseconds. Since I count 1 as a power, my answer is 1 greater than yours.

If you take out the print() statement, replace floor(n^(1/i)) with sqrtnint(n,i), and (if you want) take out the s++ statement at the end if you want to ignore the power 1, and run it on your machine, you'll probably do 10^50 in less than a millisecond, and also be able to do arbitrary n. With the default precision, my script craps out at around 10^77 due to insufficient precision.

n=10^50;
? 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 9999999999999999999999999 9999999999999999999999999
3 10000000046415888336127786 46415888336127787
5 10000000046415898336127784 9999999998
6 10000000046415898120684316 -215443468
7 10000000046415898134579269 13894953
10 10000000046415898134479271 -99998
11 10000000046415898134514381 35110
13 10000000046415898134521397 7016
14 10000000046415898134517671 -3726
15 10000000046415898134515518 -2153
17 10000000046415898134516390 872
19 10000000046415898134516817 427
21 10000000046415898134516578 -239
22 10000000046415898134516392 -186

<cut a whole bunch of printed lines>

163 10000000046415898134516526 1
165 10000000046415898134516527 1
166 10000000046415898134516526 -1
167 10000000046415898134516526 0
time = 2 ms.
%2 = 10000000046415898134516527