Is there a power-counting function?
 

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)

This sounds like a honest question. :razz:

There are:

[TEX]\lfloor\sqrt x-1\rfloor+\lfloor\sqrt[3]x-1\rfloor+\lfloor\sqrt[4]x-1\rfloor+...+\lfloor\sqrt[\lfloor\log_2x\rfloor]x-1\rfloor[/TEX]

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

Here is the OEIS sequence for it:

[url]http://oeis.org/A070428[/url]

[url]http://oeis.org/A070428/b070428.txt[/url]

[CODE]
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 >
[/CODE](it seems like my boss is not pushing me too hard today... :razz:)

A faster version using inclusion/exclusion principle rather than recursion:
[CODE]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[/CODE]

EDIT:- Apparently I reinvented Möbius function

[QUOTE=enzocreti;537486]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)[/QUOTE]

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

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[/code]

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[/code]

The print() statement can of course be eliminated. I merely put it in to show the "inclusion-exclusion" at work.

If n is [i]very[/i] 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.

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*** :razz:
-------------
**challenge: the fastest "nextpower()" and "precpower()" functions?
***you may however, need to redefine realprecision to match the range in this case

[QUOTE=LaurV;537614]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.[/QUOTE]
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.
[CODE]? \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.
?
[/CODE]

Whoops... :redface:
fast reading, you confused me with those "n" and "i" mixture :razz:
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, :razz: 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)

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.

[code]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[/code]
<cut a whole bunch of printed lines>
[code]163 10000000046415898134516526 1
165 10000000046415898134516527 1
166 10000000046415898134516526 -1
167 10000000046415898134516526 0
time = 2 ms.
%2 = 10000000046415898134516527[/code]