Sorry, an error occurred while loading the content.

## Re: [PrimeNumbers] Digest Number 2715

Expand Messages
• ... 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 10:47 AM
• 0 Attachment
primenumbers@yahoogroups.com wrote:
> ________________________________________________________________________
> ________________________________________________________________________
> 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.
>
>
>

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 2 of 5 , Jun 9 1:16 PM
• 0 Attachment
> 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 3 of 5 , Jun 9 1:36 PM
• 0 Attachment
> 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 4 of 5 , Jun 10 1:56 PM
• 0 Attachment
--- On Tue, 6/9/09, Maximilian Hasler <maximilian.hasler@...> wrote:
> From: Maximilian Hasler <maximilian.hasler@...>
> 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.

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 5 of 5 , Jun 10 5:39 PM
• 0 Attachment
--- 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.