mersenneforum.org  

Go Back   mersenneforum.org > Extra Stuff > Blogorrhea > enzocreti

Reply
 
Thread Tools
Old 2020-02-13, 08:13   #1
enzocreti
 
Mar 2018

17·31 Posts
Default 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)
enzocreti is offline   Reply With Quote
Old 2020-02-14, 08:12   #2
LaurV
Romulan Interpreter
 
LaurV's Avatar
 
Jun 2011
Thailand

22AC16 Posts
Default

This sounds like a honest question.

There are:

\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

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

Last fiddled with by LaurV on 2020-02-14 at 08:52 Reason: forgot the lfoors and rfloors
LaurV is offline   Reply With Quote
Old 2020-02-14, 09:06   #3
ATH
Einyen
 
ATH's Avatar
 
Dec 2003
Denmark

2,969 Posts
Default

Here is the OEIS sequence for it:

http://oeis.org/A070428

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

Last fiddled with by ATH on 2020-02-14 at 09:07
ATH is offline   Reply With Quote
Old 2020-02-14, 09:12   #4
LaurV
Romulan Interpreter
 
LaurV's Avatar
 
Jun 2011
Thailand

22×7×317 Posts
Default

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 >
(it seems like my boss is not pushing me too hard today... )

Last fiddled with by LaurV on 2020-02-14 at 09:13
LaurV is offline   Reply With Quote
Old 2020-02-14, 09:49   #5
axn
 
axn's Avatar
 
Jun 2003

472010 Posts
Default

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
EDIT:- Apparently I reinvented Möbius function

Last fiddled with by axn on 2020-02-14 at 10:15
axn is offline   Reply With Quote
Old 2020-02-14, 12:15   #6
sweety439
 
sweety439's Avatar
 
Nov 2016

2×3×5×79 Posts
Default

Quote:
Originally Posted by enzocreti View Post
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
sweety439 is offline   Reply With Quote
Old 2020-02-14, 14:51   #7
Dr Sardonicus
 
Dr Sardonicus's Avatar
 
Feb 2017
Nowhere

358310 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
Old 2020-02-15, 02:14   #8
LaurV
Romulan Interpreter
 
LaurV's Avatar
 
Jun 2011
Thailand

22×7×317 Posts
Default

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

Last fiddled with by LaurV on 2020-02-15 at 03:37
LaurV is offline   Reply With Quote
Old 2020-02-15, 04:09   #9
axn
 
axn's Avatar
 
Jun 2003

472010 Posts
Default

Quote:
Originally Posted by LaurV View Post
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.
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.
?
axn is offline   Reply With Quote
Old 2020-02-15, 07:38   #10
LaurV
Romulan Interpreter
 
LaurV's Avatar
 
Jun 2011
Thailand

100010101011002 Posts
Default

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)

Last fiddled with by LaurV on 2020-02-15 at 07:51
LaurV is offline   Reply With Quote
Old 2020-02-15, 13:23   #11
Dr Sardonicus
 
Dr Sardonicus's Avatar
 
Feb 2017
Nowhere

67778 Posts
Default

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
<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
Dr Sardonicus is offline   Reply With Quote
Reply

Thread Tools


Similar Threads
Thread Thread Starter Forum Replies Last Post
New confirmed pi(10^27), pi(10^28) prime counting function records kwalisch Computer Science & Computational Number Theory 34 2020-10-29 10:04
Prime counting function records D. B. Staple Computer Science & Computational Number Theory 48 2020-10-28 07:54
Prime counting function Steve One Miscellaneous Math 20 2018-03-03 22:44
Legendre's prime counting function pbewig Information & Answers 0 2011-07-14 00:47
Counting on Fingers Orgasmic Troll Lounge 3 2005-12-31 22:22

All times are UTC. The time now is 02:47.

Fri Oct 30 02:47:27 UTC 2020 up 49 days, 23:58, 1 user, load averages: 2.07, 2.08, 2.15

Powered by vBulletin® Version 3.8.11
Copyright ©2000 - 2020, Jelsoft Enterprises Ltd.

This forum has received and complied with 0 (zero) government requests for information.

Permission is granted to copy, distribute and/or modify this document under the terms of the GNU Free Documentation License, Version 1.2 or any later version published by the Free Software Foundation.
A copy of the license is included in the FAQ.