Browse Groups

• ... The problem is not with either function. The problem is that, in any computer, real numbers are represented by dyadic rationals. Try
Message 1 of 5 , Jun 9, 2009
View Source
> ________________________________________________________________________
> ________________________________________________________________________
> 2a. Bug in Pari gp ?
> Date: Tue Jun 9, 2009 12:11 am ((PDT))
>
> ? for(n=1,30,print([n,ceil(sqrtn(n^2,2))]))
> [1, 1]
>
> [9, 10]
> [10, 11]
> [11, 11]
>
> [14, 15]
> [15, 16]
> [16, 16]
>
> [18, 19]
> [19, 20]
> [20, 21]
> [21, 21]
>
> [27, 28]
> [28, 29]
> [29, 30]
> [30, 31]
>

> ? for(n=1,30,print([n,floor(sqrtn(n^2,2))]))
> [1, 1]
>
> [10, 10]
> [11, 10]
> [12, 11]
> [13, 13]
>
> [22, 21]
> [23, 23]
> [24, 23]
> [25, 25]
>
>
> Is it a problem with the sqrtn(x,n) (n-th root of x)? Or wih floor and ceil functions?
>
> I wanted to search if some numbers were perfect powers using the trick:
> If(floor(sqrt(x,n))==ceil(sqrt(x,n)),.....
>
> and trying to use n (the n-th root ) as a variable in a for loop to look for many possible n-roots if this was possible technically.
>
>
>

The problem is not with either function. The problem is that, in any
computer, real numbers are represented by dyadic rationals.

Try

for(n=1,30,print([n,ceil( -.001 + sqrtn(n^2,2))]))

for(n=1,30,print([n,floor(.001 + sqrtn(n^2,2))]))

> What is the way to solve my problem to detect if some kind of numbers are perfect powers?
>
>
>

Of course the easiest solution is that suggested by Pete Kosinar.

Use ispower()

Kermit
• ... Well, but if they are integers, it is conceivable that the results be exact in some range largely exceeding 30^2. Btw, using sqrt() instead of sqrtn(,2)
Message 1 of 5 , Jun 9, 2009
View Source
> The problem is not with either function. The problem is that, in any
> computer, real numbers are represented by dyadic rationals.

Well, but if they are integers, it is conceivable that the results be exact in some range largely exceeding 30^2.

Btw, using sqrt() instead of sqrtn(,2) gives much less "errors".

> Try
> for(n=1,30,print([n,ceil( -.001 + sqrtn(n^2,2))]))
> for(n=1,30,print([n,floor(.001 + sqrtn(n^2,2))]))

no need to be that shy, you can as well use
ceil( sqrtn( n^2-.5, 2))
or
floor( sqrtn(n^2+1, 2))
to get the mathematically correct result over a larger range.

In place of the latter, you can also use
sqrtint(n^2+n)
since the next higher square is (n+1)^2 = nÂ²+2n+1

>> What is the way to solve my problem to detect
>> if some kind of numbers are perfect powers?

An elementary way without using the ispower() function is to check if the exponents in the prime factorization have a gcd > 1 :
If n = m^k with m = product p_i ^ e_i,
then all exponents f_i in n = product p_i ^ f_i are multiples of k.

for( n=1,99, gcd( factor(n)[,2] )>1 & print1(n","))
4,8,9,16,25,27,32,36,49,64,81,
=> http://www.research.att.com/~njas/sequences/A001597

Of course the ispower function does not need to calculate the full factorization but can stop at p_i if gcd [f_1, ..., f_i] becomes = 1.

Maximilian
• ... ... so the following timings may be somehow surprising: ? for( n=1,99999, gcd( factor(n)[,2] ) 1 ) time = 136 ms. ? for( n=1,99999, ispower(n) ) time =
Message 1 of 5 , Jun 9, 2009
View Source
> Of course the ispower function does not need to calculate the full factorization but can stop at p_i if gcd [f_1, ..., f_i] becomes = 1.
>
> Maximilian

... so the following timings may be somehow surprising:

? for( n=1,99999, gcd( factor(n)[,2] )>1 )
time = 136 ms.
? for( n=1,99999, ispower(n) )
time = 2,445 ms.

? \v
GP/PARI CALCULATOR Version 2.4.2 (development)
i686 running linux (ix86 kernel) 32-bit version
compiled: Feb 11 2008, gcc-4.1.2 20061115 (prerelease) (Debian 4.1.1-21)

???

Maximilian
• ... Very much so. However the former method is roughly equal in speed at n~10^10, and gets slower at a depressing rate. The latter method slows down slowly.
Message 1 of 5 , Jun 10, 2009
View Source
--- On Tue, 6/9/09, Maximilian Hasler <maximilian.hasler@...> wrote:
> From: Maximilian Hasler <maximilian.hasler@...>
> Subject: [PrimeNumbers] Re: Digest Number 2715
> Date: Tuesday, June 9, 2009, 11:36 PM
> > Of course the ispower function
> does not need to calculate the full factorization but can
> stop at p_i if gcd [f_1, ..., f_i] becomes = 1.
> >
> > Maximilian
>
> ... so the following timings may be somehow surprising:
>
> ? for( n=1,99999, gcd( factor(n)[,2] )>1 )
> time = 136 ms.
> ? for( n=1,99999, ispower(n) )
> time = 2,445 ms.

Very much so.

However the former method is roughly equal in speed at n~10^10, and gets slower at a depressing rate. The latter method slows down slowly.

> ? \v
> GP/PARI CALCULATOR Version 2.4.2 (development)
> i686 running linux (ix86 kernel) 32-bit version
> compiled: Feb 11 2008, gcc-4.1.2 20061115 (prerelease)
> (Debian 4.1.1-21)

ix86 kernel? yikes - grab a GMP one!

Phil
• ... That s not the problem, but I forgot to post follow-up with Karim s immediate reply: actually the problem is fixed in current SVN. Grab the latest from the
Message 1 of 5 , Jun 10, 2009
View Source
--- In primenumbers@yahoogroups.com, Phil Carmody <thefatphil@...> wrote:
>> ...somehow surprising...
> Very much so.
>
> However the former method is roughly equal in speed at n~10^10, and gets slower at a depressing rate. The latter method slows down slowly.
>
> ix86 kernel? yikes - grab a GMP one!

That's not the problem, but I forgot to post follow-up with Karim's immediate reply:
actually the problem is fixed in current SVN. Grab the latest from the PARI headquarters, for me it compiled out of the box (once I added bison on my debian system).
Performance of ispower() is great now over any range, and soon drops below the gcd(factor) version:

GP/PARI CALCULATOR Version 2.4.3 (development svn-11745)
i686 running linux (ix86/GMP-4.2.4 kernel) 32-bit version
...
(22:46) gp > for( n=1,99999, ispower(n) )
time = 172 ms.
(22:46) gp > for( n=1,99999, gcd( factor(n)[,2] )>1 )
time = 160 ms.
(22:46) gp > for( n=10^8,10^8+100000, gcd( factor(n)[,2] )>1 )
time = 757 ms.
(22:47) gp > for( n=10^8,10^8+100000, ispower(n) )
time = 188 ms.

Maximilian
Your message has been successfully submitted and would be delivered to recipients shortly.
• Changes have not been saved
Press OK to abandon changes or Cancel to continue editing
• Your browser is not supported
Kindly note that Groups does not support 7.0 or earlier versions of Internet Explorer. We recommend upgrading to the latest Internet Explorer, Google Chrome, or Firefox. If you are using IE 9 or later, make sure you turn off Compatibility View.