 mersenneforum.org Is there a power-counting function?
 Register FAQ Search Today's Posts Mark Forums Read 2020-02-13, 08:13 #1 enzocreti   Mar 2018 2×3×7×11 Posts 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)   2020-02-14, 08:12 #2 LaurV Romulan Interpreter   Jun 2011 Thailand 3×31×89 Posts 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. Last fiddled with by LaurV on 2020-02-14 at 08:52 Reason: forgot the lfoors and rfloors   2020-02-14, 09:06 #3 ATH Einyen   Dec 2003 Denmark 3·13·71 Posts 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   2020-02-14, 09:12 #4 LaurV Romulan Interpreter   Jun 2011 Thailand 3×31×89 Posts Code: gp > cnt=0;for(i=1,100000,if(ispower(i),print(cnt++": "i))) ..... .... 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   2020-02-14, 09:49 #5 axn   Jun 2003 22·3·373 Posts 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   2020-02-14, 12:15   #6
sweety439

Nov 2016

110111011102 Posts Quote:
 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   2020-02-14, 14:51 #7 Dr Sardonicus   Feb 2017 Nowhere 25·3·29 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   2020-02-15, 02:14 #8 LaurV Romulan Interpreter   Jun 2011 Thailand 205516 Posts 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   2020-02-15, 04:09   #9
axn

Jun 2003

22×3×373 Posts Quote:
 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.
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.
?   2020-02-15, 07:38 #10 LaurV Romulan Interpreter   Jun 2011 Thailand 3×31×89 Posts 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   2020-02-15, 13:23 #11 Dr Sardonicus   Feb 2017 Nowhere 53408 Posts 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: 163 10000000046415898134516526 1 165 10000000046415898134516527 1 166 10000000046415898134516526 -1 167 10000000046415898134516526 0 time = 2 ms. %2 = 10000000046415898134516527  Thread Tools Show Printable Version Email this Page Similar Threads Thread Thread Starter Forum Replies Last Post kwalisch Computer Science & Computational Number Theory 22 2018-09-09 17:25 Steve One Miscellaneous Math 20 2018-03-03 22:44 D. B. Staple Computer Science & Computational Number Theory 43 2016-12-25 22:15 pbewig Information & Answers 0 2011-07-14 00:47 Orgasmic Troll Lounge 3 2005-12-31 22:22

All times are UTC. The time now is 23:55.

Sun Feb 23 23:55:29 UTC 2020 up 23 days, 18:27, 2 users, load averages: 2.26, 2.46, 2.52