Pierpont primes

From Rosetta Code
Revision as of 11:57, 22 August 2019 by Thundergnat (talk | contribs) (→‎{{header|Perl 6}}: further simplifications)
Pierpont primes is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

A Pierpont prime is a prime number of the form: 2u3v + 1 for some non-negative integers u and v .


A Pierpont prime of the second kind is a prime number of the form: 2u3v - 1 for some non-negative integers u and v .


The term "Pierpont primes" is generally understood to mean the first definition, but will be called "Pierpont primes of the first kind" on this page to distinguish them.


Task
  • Write a routine (function, procedure, whatever) to find Pierpont primes of the first & second kinds.
  • Use the routine to find and display here, on this page, the first 50 Pierpont primes of the first kind.
  • Use the routine to find and display here, on this page, the first 50 Pierpont primes of the second kind
  • If your language supports large integers, find and display here, on this page, the 250th Pierpont prime of the first kind and the 250th Pierpont prime of the second kind.


See also


Factor

<lang factor>USING: fry grouping io kernel locals make math math.functions math.primes prettyprint sequences sorting ;

pierpont ( ulim vlim quot -- seq )
   '[
       _ <iota> _ <iota> [
           [ 2 ] [ 3 ] bi* [ swap ^ ] 2bi@ * 1 @
           dup prime? [ , ] [ drop ] if
       ] cartesian-each
   ] { } make natural-sort ; inline
.fifty ( seq -- ) 50 head 10 group simple-table. nl ;

[let

   [ + ] [ - ] [ [ 120 80 ] dip pierpont ] bi@
   :> ( first second )
   
   "First 50 Pierpont primes of the first kind:" print
   first .fifty
   "First 50 Pierpont primes of the second kind:" print
   second .fifty
   "250th Pierpont prime of the first kind: " write
   249 first nth . nl
   "250th Pierpont prime of the second kind: " write
   249 second nth .

]</lang>

Output:
First 50 Pierpont primes of the first kind:
2      3      5       7       13      17      19      37      73      97
109    163    193     257     433     487     577     769     1153    1297
1459   2593   2917    3457    3889    10369   12289   17497   18433   39367
52489  65537  139969  147457  209953  331777  472393  629857  746497  786433
839809 995329 1179649 1492993 1769473 1990657 2654209 5038849 5308417 8503057

First 50 Pierpont primes of the second kind:
2      3      5      7       11      17       23       31       47       53
71     107    127    191     383     431      647      863      971      1151
2591   4373   6143   6911    8191    8747     13121    15551    23327    27647
62207  73727  131071 139967  165887  294911   314927   442367   472391   497663
524287 786431 995327 1062881 2519423 10616831 17915903 18874367 25509167 30233087

250th Pierpont prime of the first kind: 62518864539857068333550694039553

250th Pierpont prime of the second kind: 4111131172000956525894875083702271

Go

Brute force approach

Despite being inherently inefficient, this still works very quickly (less than 0.4 seconds on my machine).

A GMP wrapper, rather than Go's math/big package, has been used for extra speed (about 3.5 times quicker).

However, in order to be sure that the first 250 Pierpont primes will be generated, it is necessary to choose the loop sizes to produce somewhat more than this and then sort them into order. <lang go>package main

import (

   "fmt"
   big "github.com/ncw/gmp"
   "sort"

)

var (

   one   = new(big.Int).SetUint64(1)
   two   = new(big.Int).SetUint64(2)
   three = new(big.Int).SetUint64(3)

)

func pierpont(ulim, vlim int, first bool) []*big.Int {

   p := new(big.Int)
   p2 := new(big.Int).Set(one)
   p3 := new(big.Int).Set(one)   
   var pp []*big.Int
   for v := 0; v < vlim; v++ {
       for u := 0; u < ulim; u++ {
           p.Mul(p2, p3)
           if first {
               p.Add(p, one)
           } else {
               p.Sub(p, one)
           }
           if p.ProbablyPrime(10) {
               q := new(big.Int)
               q.Set(p)
               pp = append(pp, q)
           }
           p2.Mul(p2, two)
       }
       p3.Mul(p3, three)
       p2.Set(one)
   }
   sort.Slice(pp, func(i, j int) bool {
       return pp[i].Cmp(pp[j]) < 0
   })
   return pp

} func main() {

   fmt.Println("First 50 Pierpont primes of the first kind:")
   pp := pierpont(120, 80, true)
   for i := 0; i < 50; i++ {
       fmt.Printf("%8d ", pp[i])
       if (i-9)%10 == 0 {
           fmt.Println()
       }
   }
   fmt.Println("\nFirst 50 Pierpont primes of the second kind:")
   pp2 := pierpont(120, 80, false)
   for i := 0; i < 50; i++ {
       fmt.Printf("%8d ", pp2[i])
       if (i-9)%10 == 0 {
           fmt.Println()
       }
   }
   fmt.Println("\n250th Pierpont prime of the first kind:", pp[249])
   fmt.Println("\n250th Pierpont prime of the second kind:", pp2[249])

}</lang>

Output:
First 50 Pierpont primes of the first kind:
       2        3        5        7       13       17       19       37       73       97 
     109      163      193      257      433      487      577      769     1153     1297 
    1459     2593     2917     3457     3889    10369    12289    17497    18433    39367 
   52489    65537   139969   147457   209953   331777   472393   629857   746497   786433 
  839809   995329  1179649  1492993  1769473  1990657  2654209  5038849  5308417  8503057 

First 50 Pierpont primes of the second kind:
       2        3        5        7       11       17       23       31       47       53 
      71      107      127      191      383      431      647      863      971     1151 
    2591     4373     6143     6911     8191     8747    13121    15551    23327    27647 
   62207    73727   131071   139967   165887   294911   314927   442367   472391   497663 
  524287   786431   995327  1062881  2519423 10616831 17915903 18874367 25509167 30233087 

250th Pierpont prime of the first kind: 62518864539857068333550694039553

250th Pierpont prime of the second kind: 4111131172000956525894875083702271

3-Smooth approach

The strategy here is to generate successive 3-smooth numbers, add (or subtract) one, check if prime and, if so, append to a slice of 'big' integers until the required number is reached.

The Pierpoint primes of the first and second kind are generated at the same time so there is no real need for parallel processing.

This approach is considerably faster than the first version. Although not shown below, it produces the first 250 Pierpont primes (of both kinds) in under 0.2 seconds and the first 1000 in about 7.4 seconds. However, the first 2000 takes around 100 seconds.

These timings are for my Celeron @1.6GHz and should therefore be much faster on a more modern machine.

<lang go>package main

import (

   "fmt"
   big "github.com/ncw/gmp"
   "time"

)

var (

   one   = new(big.Int).SetUint64(1)
   two   = new(big.Int).SetUint64(2)
   three = new(big.Int).SetUint64(3)

)

func min(i, j int) int {

   if i < j {
       return i
   }
   return j

}

func pierpont(n int, first bool) (p [2][]*big.Int) {

   p[0] = make([]*big.Int, n)
   p[1] = make([]*big.Int, n)
   for i := 0; i < n; i++ {
       p[0][i] = new(big.Int)
       p[1][i] = new(big.Int)
   }
   p[0][0].Set(two)
   count, count1, count2 := 0, 1, 0
   var s []*big.Int
   s = append(s, new(big.Int).Set(one))
   i2, i3, k := 0, 0, 1
   n2, n3, t := new(big.Int), new(big.Int), new(big.Int)
   for count < n {
       n2.Mul(s[i2], two)
       n3.Mul(s[i3], three)
       if n2.Cmp(n3) < 0 {
           t.Set(n2)
           i2++
       } else {
           t.Set(n3)
           i3++
       }
       if t.Cmp(s[k-1]) > 0 {
           s = append(s, new(big.Int).Set(t))
           k++
           t.Add(t, one)
           if count1 < n && t.ProbablyPrime(10) {
               p[0][count1].Set(t)
               count1++
           }
           t.Sub(t, two)
           if count2 < n && t.ProbablyPrime(10) {
               p[1][count2].Set(t)
               count2++
           }
           count = min(count1, count2)
       }
   }
   return p

}

func main() {

   start := time.Now()
   p := pierpont(2000, true)
   fmt.Println("First 50 Pierpont primes of the first kind:")
   for i := 0; i < 50; i++ {
       fmt.Printf("%8d ", p[0][i])
       if (i-9)%10 == 0 {
           fmt.Println()
       }
   }
   fmt.Println("\nFirst 50 Pierpont primes of the second kind:")
   for i := 0; i < 50; i++ {
       fmt.Printf("%8d ", p[1][i])
       if (i-9)%10 == 0 {
           fmt.Println()
       }
   }
   fmt.Println("\n250th Pierpont prime of the first kind:", p[0][249])
   fmt.Println("\n250th Pierpont prime of the second kind:", p[1][249])
   
   fmt.Println("\n1000th Pierpont prime of the first kind:", p[0][999])
   fmt.Println("\n1000th Pierpont prime of the second kind:", p[1][999])
   
   fmt.Println("\n2000th Pierpont prime of the first kind:", p[0][1999])
   fmt.Println("\n2000th Pierpont prime of the second kind:", p[1][1999])
   
   elapsed := time.Now().Sub(start)
   fmt.Printf("\nTook %s\n", elapsed)

}</lang>

Output:
First 50 Pierpont primes of the first kind:
       2        3        5        7       13       17       19       37       73       97 
     109      163      193      257      433      487      577      769     1153     1297 
    1459     2593     2917     3457     3889    10369    12289    17497    18433    39367 
   52489    65537   139969   147457   209953   331777   472393   629857   746497   786433 
  839809   995329  1179649  1492993  1769473  1990657  2654209  5038849  5308417  8503057 

First 50 Pierpont primes of the second kind:
       2        3        5        7       11       17       23       31       47       53 
      71      107      127      191      383      431      647      863      971     1151 
    2591     4373     6143     6911     8191     8747    13121    15551    23327    27647 
   62207    73727   131071   139967   165887   294911   314927   442367   472391   497663 
  524287   786431   995327  1062881  2519423 10616831 17915903 18874367 25509167 30233087 

250th Pierpont prime of the first kind: 62518864539857068333550694039553

250th Pierpont prime of the second kind: 4111131172000956525894875083702271

1000th Pierpont prime of the first kind: 69269314716439690250482558089997110961545818230232043107188537422260188701607997086273960899938499201024414931399264696270849

1000th Pierpont prime of the second kind: 1308088756227965581249669045506775407896673213729433892383353027814827286537163695213418982500477392209371001259166465228280492460735463423

2000th Pierpont prime of the first kind: 23647056334818750458979408107288138983957799805326855934519920502493109431728722178351835778368596067773810122477389192659352731519830867553659739507195398662712180250483714053474639899675114018023738461139103130959712720686117399642823861502738433

2000th Pierpont prime of the second kind: 1702224134662426018061116932011222570937093650174807121918750428723338890211147039320296240754205680537318845776107057915956535566573559841027244444877454493022783449689509569107393738917120492483994302725479684822283929715327187974256253064796234576415398735760543848603844607

Took 1m40.781726122s

Julia

The generator method is very fast but does not guarantee the primes are generated in order. Therefore we generate two times the primes needed and then sort and return the lower half. <lang julia>using Primes

function pierponts(N, firstkind = true)

   ret, incdec = BigInt[], firstkind ? 1 : -1
   for k2 in 0:10000, k3 in 0:k2, switch in false:true
       i, j = switch ? (k3, k2) : (k2, k3)
       n = BigInt(2)^i * BigInt(3)^j + incdec
       if isprime(n) && !(n in ret)
           push!(ret, n)
           if length(ret) == N * 2
               return sort(ret)[1:N]
           end
       end
   end
   throw("Failed to find $(N * 2) primes")

end

println("The first 50 Pierpont primes (first kind) are: ", pierponts(50))

println("\nThe first 50 Pierpont primes (second kind) are: ", pierponts(50, false))

println("\nThe 250th Pierpont prime (first kind) is: ", pierponts(250)[250])

println("\nThe 250th Pierpont prime (second kind) is: ", pierponts(250, false)[250])

println("\nThe 1000th Pierpont prime (first kind) is: ", pierponts(1000)[1000])

println("\nThe 1000th Pierpont prime (second kind) is: ", pierponts(1000, false)[1000])

println("\nThe 2000th Pierpont prime (first kind) is: ", pierponts(2000)[2000])

println("\nThe 2000th Pierpont prime (second kind) is: ", pierponts(2000, false)[2000])

</lang>

Output:
The first 50 Pierpont primes (first kind) are: BigInt[2, 3, 5, 7, 13, 17, 19, 37, 73, 97, 109, 163, 193, 257, 433, 487, 577, 769, 1153, 1297, 1459, 2593, 2917, 3457, 3889, 10369, 12289, 17497, 18433, 39367, 52489, 65537, 139969, 147457, 209953, 331777, 472393, 629857, 746497, 786433, 839809, 995329, 1179649, 1492993, 1769473, 1990657, 2654209, 5038849, 5308417, 8503057]

The first 50 Pierpont primes (second kind) are: BigInt[2, 3, 5, 7, 11, 17, 23, 31, 47, 53, 71, 107, 127, 191, 383, 431, 647, 863, 971, 1151, 2591, 4373, 6143, 6911, 8191, 8747, 13121, 15551, 23327, 27647, 62207, 73727, 131071, 139967, 165887, 294911, 314927, 442367, 472391, 497663, 524287, 786431, 995327, 1062881, 2519423, 10616831, 17915903, 18874367, 25509167, 30233087]

The 250th Pierpont prime (first kind) is: 62518864539857068333550694039553

The 250th Pierpont prime (second kind) is: 4111131172000956525894875083702271

The 1000th Pierpont prime (first kind) is: 69269314716439690250482558089997110961545818230232043107188537422260188701607997086273960899938499201024414931399264696270849

The 1000th Pierpont prime (second kind) is: 1308088756227965581249669045506775407896673213729433892383353027814827286537163695213418982500477392209371001259166465228280492460735463423

The 2000th Pierpont prime (first kind) is: 23647056334818750458979408107288138983957799805326855934519920502493109431728722178351835778368596067773810122477389192659352731519830867553659739507195398662712180250483714053474639899675114018023738461139103130959712720686117399642823861502738433

The 2000th Pierpont prime (second kind) is: 1702224134662426018061116932011222570937093650174807121918750428723338890211147039320296240754205680537318845776107057915956535566573559841027244444877454493022783449689509569107393738917120492483994302725479684822283929715327187974256253064796234576415398735760543848603844607

Perl 6

Works with: Rakudo version 2019.07.1

This finesse version never produces more Pierpont numbers than it needs to fulfill the requested number of primes. It uses a series of parallel iterators with additional iterators added as required. No need to speculatively generate an overabundance. No need to rely on magic numbers. No need to sort them. It produces exactly what is needed, in order, on demand.

<lang perl6>use ntheory:from<Perl5> <is_prime>;

sub pierpont ($kind is copy = 1) {

   fail "Unknown type: $kind. Must be one of 1 (default) or 2" if $kind !== 1|2;
   $kind = -1 if $kind == 2;
   my $po3 = 3;
   my $add-one = 3;
   my @iterators = [1,2,4,8 … *].iterator, [3,9,27 … *].iterator;
   my @head = @iterators».pull-one;
   gather {
       loop {
           my $key = @head.pairs.min( *.value ).key;
           my $min = @head[$key];
           @head[$key] = @iterators[$key].pull-one;
           take $min + $kind if "{$min + $kind}".&is_prime;
           if $min >= $add-one {
               @iterators.push: ([|((2,4,8).map: * * $po3) … *]).iterator;
               $add-one = @head[+@iterators - 1] = @iterators[+@iterators - 1].pull-one;
               $po3 *= 3;
           }
       }
   }

}

say "First 50 Pierpont primes of the first kind:\n" ~ pierpont[^50].rotor(10)».fmt('%8d').join: "\n";

say "\nFirst 50 Pierpont primes of the second kind:\n" ~ pierpont(2)[^50].rotor(10)».fmt('%8d').join: "\n";

say "\n250th Pierpont prime of the first kind: " ~ pierpont[249];

say "\n250th Pierpont prime of the second kind: " ~ pierpont(2)[249];

  1. And the question absolutely nobody was asking:

say "\n1000th Pierpont prime of the first kind:\n" ~ pierpont[999];

say "\n1000th Pierpont prime of the second kind:\n" ~ pierpont(2)[999];</lang>

Output:
First 50 Pierpont primes of the first kind:
       2        3        5        7       13       17       19       37       73       97
     109      163      193      257      433      487      577      769     1153     1297
    1459     2593     2917     3457     3889    10369    12289    17497    18433    39367
   52489    65537   139969   147457   209953   331777   472393   629857   746497   786433
  839809   995329  1179649  1492993  1769473  1990657  2654209  5038849  5308417  8503057

First 50 Pierpont primes of the second kind:
       2        3        5        7       11       17       23       31       47       53
      71      107      127      191      383      431      647      863      971     1151
    2591     4373     6143     6911     8191     8747    13121    15551    23327    27647
   62207    73727   131071   139967   165887   294911   314927   442367   472391   497663
  524287   786431   995327  1062881  2519423 10616831 17915903 18874367 25509167 30233087

250th Pierpont prime of the first kind: 62518864539857068333550694039553

250th Pierpont prime of the second kind: 4111131172000956525894875083702271

1000th Pierpont prime of the first kind:
69269314716439690250482558089997110961545818230232043107188537422260188701607997086273960899938499201024414931399264696270849

1000th Pierpont prime of the second kind:
1308088756227965581249669045506775407896673213729433892383353027814827286537163695213418982500477392209371001259166465228280492460735463423

REXX

The REXX language has a "big num" capability to handle almost any amount of decimal digits,   but
it lacks a robust   isPrime   function.   Without that, verifying very large primes is problematic. <lang rexx>/*REXX program finds and displays Pierpont primes of the first and second kinds. */ parse arg n . /*obtain optional argument from the CL.*/ if n== | n=="," then n= 50 /*Not specified? Then use the default.*/ numeric digits n /*ensure enough decimal digs (bit int).*/ big= copies(9, digits() ) /*BIG: used as a max number (a limit).*/

      do t=1  to -1  by -2;  usum= 0;   vsum= 0;    s= 0
      #= 0                                      /*number of Pierpont primes  (so far). */
      w= 0                                      /*the max width of a Pierpont prime.   */
      $=;    do j=0  until #>=n                 /*$: the list of the Pierpont primes.  */
             if usum<=s  then usum= get(2, 3);    if vsum<=s  then vsum= get(3, 2)
             s= min(vsum, usum);  if \isPrime(s)  then iterate  /*get min;  is prime?  */
             #= # + 1;            $= $ s                        /*bump counter; append.*/
             w= max(w, length(s) )                              /*find max prime width.*/
             end   /*j*/
      say
      if t==1  then @= '1st'                                    /*choose word for type.*/
               else @= '2nd'                                    /*   "     "   "    "  */
      say center(n   " Pierpont primes of the "   @   ' kind', max(10*(w+1) -1, 79), "═")
      call show $                                               /*display the primes.  */
      end   /*type*/

exit /*stick a fork in it, we're all done. */ /*──────────────────────────────────────────────────────────────────────────────────────*/ show: do j=1 by 10 to words($); _=

          do k=j  for 10;                 _= _ right( word($, k), w)
          end   /*k*/
        if _\==  then say substr( strip(_, 'T'), 2)
        end     /*j*/;                                        return

/*──────────────────────────────────────────────────────────────────────────────────────*/ isPrime: procedure; parse arg x; if x<2 then return 0 /*not a prime.*/

        if wordpos(x, '2 3 5 7')\==0     then return 1                  /*it's a prime.*/
        if x//2==0  then return 0;       if x//3==0  then return 0      /*not  a prime.*/
          do j=5  by 6  until j*j>x
          if x//j==0  then return 0;     if x//(j+2)==0  then return 0  /*not  a prime.*/
          end   /*j*/;                                        return 1  /*it's a prime.*/

/*──────────────────────────────────────────────────────────────────────────────────────*/ get: parse arg c1,c2; m=big; do ju=0; pu= c1**ju; if pu+t>s then return min(m, pu+t)

                               do jv=0;  pv= c2**jv;  if pv  >s  then iterate ju
                               _= pu*pv + t;          if _   >s  then m= min(_, m)
                               end   /*jv*/
                             end     /*ju*/           /*see the    RETURN    (above).  */</lang>
output   when using the default input:
═════════════════════50  Pierpont primes of the  1st  kind═════════════════════
      2       3       5       7      13      17      19      37      73      97
    109     163     193     257     433     487     577     769    1153    1297
   1459    2593    2917    3457    3889   10369   12289   17497   18433   39367
  52489   65537  139969  147457  209953  331777  472393  629857  746497  786433
 839809  995329 1179649 1492993 1769473 1990657 2654209 5038849 5308417 8503057

══════════════════════════50  Pierpont primes of the  2nd  kind══════════════════════════
       2        3        5        7       11       17       23       31       47       53
      71      107      127      191      383      431      647      863      971     1151
    2591     4373     6143     6911     8191     8747    13121    15551    23327    27647
   62207    73727   131071   139967   165887   294911   314927   442367   472391   497663
  524287   786431   995327  1062881  2519423 10616831 17915903 18874367 25509167 30233087

zkl

Translation of: Julia
Library: GMP

GNU Multiple Precision Arithmetic Library

Using GMP's probabilistic primes makes it is easy and fast to test for primeness. <lang zkl>var [const] BI=Import("zklBigNum"); // libGMP

fcn pierPonts(N, firstkind=True){

  var pps1=List(), pps2=List();	// cache
  pps,incDec := ( if(firstkind) pps1 else pps2 ), firstkind and 1 or -1;
  foreach k2,k3,swap in ([0..10_000], [0..k2], T(False,True)){
     if(swap) i,j := k3,k2; else i,j := k2,k3;
     n:=BI(2).pow(i) * BI(3).pow(j) + incDec;
     if(n.probablyPrime() and not pps.holds(n)){ // prime testing 1st a bit faster

pps.merge(n); // keep list sorted if(pps.len() == N*2) return(pps[0,N]);

     }
  }
  throw(Exception.NotFoundError("Failed to find %d primes".fmt(N * 2)));

}</lang> <lang zkl>println("The first 50 Pierpont primes (first kind): ", pierPonts(50).concat(", ")); println("\nThe first 50 Pierpont primes (second kind): ", pierPonts(50, False).concat(", "));

foreach n in (T(250, 1_000, 2_000)){ // slow, this is why we cached

  println("\n%4dth Pierpont prime, first kind: ".fmt(n), pierPonts(n)[-1]);
  println( "                      second kind: ", pierPonts(n,False)[-1])

}</lang>

Output:
The first 50 Pierpont primes (first kind): 2, 3, 5, 7, 13, 17, 19, 37, 73, 97, 109, 163, 193, 257, 433, 487, 577, 769, 1153, 1297, 1459, 2593, 2917, 3457, 3889, 10369, 12289, 17497, 18433, 39367, 52489, 65537, 139969, 147457, 209953, 331777, 472393, 629857, 746497, 786433, 839809, 995329, 1179649, 1492993, 1769473, 1990657, 2654209, 5038849, 5308417, 8503057

The first 50 Pierpont primes (second kind): 2, 3, 5, 7, 11, 17, 23, 31, 47, 53, 71, 107, 127, 191, 383, 431, 647, 863, 971, 1151, 2591, 4373, 6143, 6911, 8191, 8747, 13121, 15551, 23327, 27647, 62207, 73727, 131071, 139967, 165887, 294911, 314927, 442367, 472391, 497663, 524287, 786431, 995327, 1062881, 2519423, 10616831, 17915903, 18874367, 25509167, 30233087

 250th Pierpont prime, first kind: 62518864539857068333550694039553
                      second kind: 4111131172000956525894875083702271

1000th Pierpont prime, first kind: 69269314716439690250482558089997110961545818230232043107188537422260188701607997086273960899938499201024414931399264696270849
                      second kind: 1308088756227965581249669045506775407896673213729433892383353027814827286537163695213418982500477392209371001259166465228280492460735463423

2000th Pierpont prime, first kind: 23647056334818750458979408107288138983957799805326855934519920502493109431728722178351835778368596067773810122477389192659352731519830867553659739507195398662712180250483714053474639899675114018023738461139103130959712720686117399642823861502738433
                      second kind: 1702224134662426018061116932011222570937093650174807121918750428723338890211147039320296240754205680537318845776107057915956535566573559841027244444877454493022783449689509569107393738917120492483994302725479684822283929715327187974256253064796234576415398735760543848603844607