- primenumbers@yahoogroups.com wrote:
> ________________________________________________________________________

The problem is not with either function. The problem is that, in any

> ________________________________________________________________________

> 2a. Bug in Pari gp ?

> Posted by: "sopadeajo2001" sopadeajo2001@... sopadeajo2001

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

>

>

>

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 > The problem is not with either function. The problem is that, in any

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

> computer, real numbers are represented by dyadic rationals.

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

> Try

no need to be that shy, you can as well use

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

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

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

An elementary way without using the ispower() function is to check if the exponents in the prime factorization have a gcd > 1 :

>> if some kind of numbers are perfect powers?

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

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

>

> Maximilian

? 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- --- On Tue, 6/9/09, Maximilian Hasler <maximilian.hasler@...> wrote:
> From: Maximilian Hasler <maximilian.hasler@...>

Very much so.

> Subject: [PrimeNumbers] Re: Digest Number 2715

> To: primenumbers@yahoogroups.com

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

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

ix86 kernel? yikes - grab a GMP one!

> 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)

Phil - --- In primenumbers@yahoogroups.com, Phil Carmody <thefatphil@...> wrote:
>> ...somehow surprising...

That's not the problem, but I forgot to post follow-up with Karim's immediate reply:

> 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!

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