Loading ...
Sorry, an error occurred while loading the content.

Re: [PrimeNumbers] Digest Number 2715

Expand Messages
  • Kermit Rose
    ... 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
    • Maximilian Hasler
      ... 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
      • Maximilian Hasler
        ... ... 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
        • Phil Carmody
          ... 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
          • Maximilian Hasler
            ... 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.