Cuban primes

From Rosetta Code
Revision as of 22:25, 1 February 2019 by rosettacode>Horst.h (added Pascal)
Cuban 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.

The name   cuban   has nothing to do with Cuba,   but has to do with the fact that cubes   (3rd powers)   play a role in its definition.


Some definitions of cuban primes
  •   primes which are the difference of two consecutive cubes.
  •   primes of the form:   (n+1)3 - n3.
  •   primes of the form:   n3 - (n-1)3.
  •   primes   p   such that   n2(p+n)   is a cube for some   n>0.
  •   primes   p   such that   4p = 1 + 3n2.


Cuban primes were named in 1923 by Allan Joseph Champneys Cunningham.


Task requirements
  •   show the first   200   cuban primes   (in a multi─line horizontal format).
  •   show the   100,000th   cuban prime.
  •   show all cuban primes with commas.
  •   show all output here.


Also see
  •   MathWorld entry:   cuban prime.
  •   The OEIS entry:   A2407.     The   100,000th   cuban prime can be verified in the   2nd   example   on this OEIS web page.



Go

<lang go>package main

import ( 

   "fmt"
   "math/big"

)

func commatize(n uint64) string { 

   s := fmt.Sprintf("%d", n)
   le := len(s)
   for i := le - 3; i >= 1; i -= 3 {
       s = s[0:i] + "," + s[i:]
   }
   return s

}

func main() { 

   var z big.Int
   var cube1, cube2, cube100k, diff uint64
   cubans := make([]string, 200)
   cube1 = 1
   count := 0
   for i := 1; ; i++ {
       j := i + 1
       cube2 = uint64(j * j * j)
       diff = cube2 - cube1
       z.SetUint64(diff)
       if z.ProbablyPrime(0) { // 100% accurate for z < 2 ^ 64
           if count < 200 {
               cubans[count] = commatize(diff)
           }
           count++
           if count == 100000 {
               cube100k = diff
               break
           }
       }
       cube1 = cube2
   }
   fmt.Println("The first 200 cuban primes are:-")
   for i := 0; i < 20; i++ {
       j := i * 10 
       fmt.Printf("%9s\n", cubans[j : j+10]) // 10 per line say
   }
   fmt.Println("\nThe 100,000th cuban prime is", commatize(cube100k))

}</lang>

Output:
The first 200 cuban primes are:-
[        7        19        37        61       127       271       331       397       547       631]
[      919     1,657     1,801     1,951     2,269     2,437     2,791     3,169     3,571     4,219]
[    4,447     5,167     5,419     6,211     7,057     7,351     8,269     9,241    10,267    11,719]
[   12,097    13,267    13,669    16,651    19,441    19,927    22,447    23,497    24,571    25,117]
[   26,227    27,361    33,391    35,317    42,841    45,757    47,251    49,537    50,311    55,897]
[   59,221    60,919    65,269    70,687    73,477    74,419    75,367    81,181    82,171    87,211]
[   88,237    89,269    92,401    96,661   102,121   103,231   104,347   110,017   112,327   114,661]
[  115,837   126,691   129,169   131,671   135,469   140,617   144,541   145,861   151,201   155,269]
[  163,567   169,219   170,647   176,419   180,811   189,757   200,467   202,021   213,067   231,019]
[  234,361   241,117   246,247   251,431   260,191   263,737   267,307   276,337   279,991   283,669]
[  285,517   292,969   296,731   298,621   310,087   329,677   333,667   337,681   347,821   351,919]
[  360,187   368,551   372,769   374,887   377,011   383,419   387,721   398,581   407,377   423,001]
[  436,627   452,797   459,817   476,407   478,801   493,291   522,919   527,941   553,411   574,219]
[  584,767   590,077   592,741   595,411   603,457   608,851   611,557   619,711   627,919   650,071]
[  658,477   666,937   689,761   692,641   698,419   707,131   733,591   742,519   760,537   769,627]
[  772,669   784,897   791,047   812,761   825,301   837,937   847,477   863,497   879,667   886,177]
[  895,987   909,151   915,769   925,741   929,077   932,419   939,121   952,597   972,991   976,411]
[  986,707   990,151   997,057 1,021,417 1,024,921 1,035,469 1,074,607 1,085,407 1,110,817 1,114,471]
[1,125,469 1,155,061 1,177,507 1,181,269 1,215,397 1,253,887 1,281,187 1,285,111 1,324,681 1,328,671]
[1,372,957 1,409,731 1,422,097 1,426,231 1,442,827 1,451,161 1,480,519 1,484,737 1,527,247 1,570,357]

The 100,000th cuban prime is 1,792,617,147,127

Pascal

Library: primTrial
Works with: Free Pascal

uses trial division to check primility.Slow in such number ranges.
OutNthCubPrime(10000) takes only 0,950 s.
100: 283,669; 1000: 65,524,807; 10000: 11,712,188,419; 100000: 1,792,617,147,127

<lang pascal>program CubanPrimes; {$IFDEF FPC} 

 {$MODE DELPHI}
 {$OPTIMIZATION ON,Regvar,PEEPHOLE,CSE,ASMCSE}
 {$CODEALIGN proc=32}

{$ENDIF} uses 

 primTrial;

const 

 COLUMNCOUNT = 10*10;

procedure FormOut(Cuban:Uint64;ColSize:Uint32); var 

 s : String;
 pI,pJ :pChar;
 i,j : NativeInt;

Begin 

 str(Cuban,s);
 i := length(s);
 If i>3 then
 Begin
   //extend s by the count of comma to be inserted
   j := i+ (i-1) div 3;
   setlength(s,j);
   pI := @s[i];
   pJ := @s[j];
   while i > 3 do
   Begin
      // copy 3 digits
      pJ^ := pI^;dec(pJ);dec(pI);
      pJ^ := pI^;dec(pJ);dec(pI);
      pJ^ := pI^;dec(pJ);dec(pI);
      // insert comma
      pJ^ := ',';dec(pJ);
      dec(i,3);
   end;
   //the digits in front are in the right place
 end;
 write(s:ColSize);

end;

procedure OutFirstCntCubPrimes(Cnt : Int32;ColCnt : Int32); var 

 cbDelta1,
 cbDelta2  : Uint64;
 ClCnt,ColSize : NativeInt;

Begin 

 If Cnt <= 0 then
   EXIT;
 IF ColCnt <= 0 then
   ColCnt := 1;
 ColSize := COLUMNCOUNT DIV ColCnt;
 dec(ColCnt);

 ClCnt := ColCnt;
 cbDelta1 := 0;
 cbDelta2 := 1;

 repeat
   if isPrime(cbDelta2) then
   Begin
     FormOut(cbDelta2,ColSize);
     dec(Cnt);

     dec(ClCnt);
     If ClCnt < 0 then
     Begin
       Writeln;
       ClCnt := ColCnt;
     end;
   end;
   inc(cbDelta1,6);// 0,6,12,18...
   inc(cbDelta2,cbDelta1);//1,7,19,35...
 until Cnt<= 0;

 writeln;

end;

procedure OutNthCubPrime(n : Int32); var 

 cbDelta1,
 cbDelta2  : Uint64;

Begin 

 If n <= 0 then
   EXIT;
 cbDelta1 := 0;
 cbDelta2 := 1;

 repeat
   inc(cbDelta1,6);
   inc(cbDelta2,cbDelta1);
   if isPrime(cbDelta2) then
     dec(n);
 until n<=0;

 FormOut(cbDelta2,20);
 writeln;

end;

Begin 

 OutFirstCntCubPrimes(200,10);
 OutNthCubPrime(100000);

end.</lang>

Output:
         7        19        37        61       127       271       331       397       547       631
       919     1,657     1,801     1,951     2,269     2,437     2,791     3,169     3,571     4,219
     4,447     5,167     5,419     6,211     7,057     7,351     8,269     9,241    10,267    11,719
    12,097    13,267    13,669    16,651    19,441    19,927    22,447    23,497    24,571    25,117
    26,227    27,361    33,391    35,317    42,841    45,757    47,251    49,537    50,311    55,897
    59,221    60,919    65,269    70,687    73,477    74,419    75,367    81,181    82,171    87,211
    88,237    89,269    92,401    96,661   102,121   103,231   104,347   110,017   112,327   114,661
   115,837   126,691   129,169   131,671   135,469   140,617   144,541   145,861   151,201   155,269
   163,567   169,219   170,647   176,419   180,811   189,757   200,467   202,021   213,067   231,019
   234,361   241,117   246,247   251,431   260,191   263,737   267,307   276,337   279,991   283,669
   285,517   292,969   296,731   298,621   310,087   329,677   333,667   337,681   347,821   351,919
   360,187   368,551   372,769   374,887   377,011   383,419   387,721   398,581   407,377   423,001
   436,627   452,797   459,817   476,407   478,801   493,291   522,919   527,941   553,411   574,219
   584,767   590,077   592,741   595,411   603,457   608,851   611,557   619,711   627,919   650,071
   658,477   666,937   689,761   692,641   698,419   707,131   733,591   742,519   760,537   769,627
   772,669   784,897   791,047   812,761   825,301   837,937   847,477   863,497   879,667   886,177
   895,987   909,151   915,769   925,741   929,077   932,419   939,121   952,597   972,991   976,411
   986,707   990,151   997,057 1,021,417 1,024,921 1,035,469 1,074,607 1,085,407 1,110,817 1,114,471
 1,125,469 1,155,061 1,177,507 1,181,269 1,215,397 1,253,887 1,281,187 1,285,111 1,324,681 1,328,671
 1,372,957 1,409,731 1,422,097 1,426,231 1,442,827 1,451,161 1,480,519 1,484,737 1,527,247 1,570,357

   1,792,617,147,127  //user    2m1.950s

Perl 6

Works with: Rakudo version 2018.12

Not the most efficient, but concise, and good enough for this task. <lang perl6>sub comma { $^i.flip.comb(3).join(',').flip }

my @cubans = lazy (1..Inf).hyper(:8degree).map({ ($_+1)³ - .³ }).grep: *.is-prime;

put @cubans[^200]».&comma».fmt("%9s").rotor(10).join: "\n";

put ;

put @cubans[99_999]., # zero indexed</lang>

Output:
        7        19        37        61       127       271       331       397       547       631
      919     1,657     1,801     1,951     2,269     2,437     2,791     3,169     3,571     4,219
    4,447     5,167     5,419     6,211     7,057     7,351     8,269     9,241    10,267    11,719
   12,097    13,267    13,669    16,651    19,441    19,927    22,447    23,497    24,571    25,117
   26,227    27,361    33,391    35,317    42,841    45,757    47,251    49,537    50,311    55,897
   59,221    60,919    65,269    70,687    73,477    74,419    75,367    81,181    82,171    87,211
   88,237    89,269    92,401    96,661   102,121   103,231   104,347   110,017   112,327   114,661
  115,837   126,691   129,169   131,671   135,469   140,617   144,541   145,861   151,201   155,269
  163,567   169,219   170,647   176,419   180,811   189,757   200,467   202,021   213,067   231,019
  234,361   241,117   246,247   251,431   260,191   263,737   267,307   276,337   279,991   283,669
  285,517   292,969   296,731   298,621   310,087   329,677   333,667   337,681   347,821   351,919
  360,187   368,551   372,769   374,887   377,011   383,419   387,721   398,581   407,377   423,001
  436,627   452,797   459,817   476,407   478,801   493,291   522,919   527,941   553,411   574,219
  584,767   590,077   592,741   595,411   603,457   608,851   611,557   619,711   627,919   650,071
  658,477   666,937   689,761   692,641   698,419   707,131   733,591   742,519   760,537   769,627
  772,669   784,897   791,047   812,761   825,301   837,937   847,477   863,497   879,667   886,177
  895,987   909,151   915,769   925,741   929,077   932,419   939,121   952,597   972,991   976,411
  986,707   990,151   997,057 1,021,417 1,024,921 1,035,469 1,074,607 1,085,407 1,110,817 1,114,471
1,125,469 1,155,061 1,177,507 1,181,269 1,215,397 1,253,887 1,281,187 1,285,111 1,324,681 1,328,671
1,372,957 1,409,731 1,422,097 1,426,231 1,442,827 1,451,161 1,480,519 1,484,737 1,527,247 1,570,357

1,792,617,147,127

REXX

<lang rexx>/*REXX program finds and displays a number of cuban primes or the Nth cuban prime. */ numeric digits 20 /*ensure enought decimal digits for #s.*/ parse arg N . /*obtain optional argument from the CL.*/ if N== | N=="," then N= 200 /*Not specified? Then use the default.*/ Nth= N<0 /*used for finding the Nth cuban prime.*/

  1. = 0 /*number of cuban primes found so far. */

!.=0;  !.0=1;  !.2=1;  !.3=1;  !.4=1;  !.5=1;  !.6=1;  !.8=1; s=41; $= sw= linesize() - 1; if sw<1 then sw=79 /*obtain width of the terminal screen. */ if (Nth & N==1) | N>0 then $= 7 /*handle case of the 1st cuban prime.*/ if (Nth & N==2) | N>1 then $= $ 19 /* " " " " 2nd " " */ if (Nth & N==3) | N>2 then $= $ 37 /* " " " " 3rd " " */

  1. = 3 /*above are needed for prime generation*/

N= abs(N) /*use absolute value of N from now on*/ if #<N then /* [↑] ending digs that aren't cubans.*/ 

    do j=4  until #=>N;          x= (j+1)**3  -  j**3        /*compute a possible cuban*/
    parse var  x      -1  _;   if !._       then iterate   /*check for the last digit*/
    if x// 3==0  then iterate;   if x// 7==0  then iterate   /*  "    " ÷ 3;  for ÷  7 */
    if x//11==0  then iterate;   if x//13==0  then iterate   /*  "    " ÷11;   "  ÷ 13 */
    if x//17==0  then iterate;   if x//19==0  then iterate   /*  "    " ÷17;   "  ÷ 19 */
    if x//23==0  then iterate;   if x//29==0  then iterate   /*  "    " ÷23;   "  ÷ 29 */
    if x//31==0  then iterate;   if x//37==0  then iterate   /*  "    " ÷31;   "  ÷ 37 */

             do k=s  by 6  until k*k>x          /*skip multiples of 3 (when using BY 6)*/
             if x// k    ==0  then iterate j    /*test for a "lower"   possible prime. */
             if x//(k+2) ==0  then iterate j    /*  "   "  " "higher"      "      "    */
             end   /*k*/
    #= # + 1                                    /*bump the number of cuban primes found*/
    if Nth then do;if #==N then do;say commas(x); leave j; end /*output is only 1 num. */
                               else iterate j                  /*keep searching.       */
                 end
    new=$ commas(x)                                            /*append a new cuban.   */
    if length(new)>sw  then do;  say $                         /*line too long, show #s*/
                                 $= commas(x)                  /*initialize next line. */
                            end
                       else $= new                             /*use this longer output*/
    end   /*j*/

if S\== then say $ /*check for any residual*/ exit /*stick a fork in it, we're all done. */ /*──────────────────────────────────────────────────────────────────────────────────────*/ commas: parse arg _; do jc=length(_)-3 to 1 by -3; _=insert(',', _, jc); end; return _</lang>

output   when using the default input of:     200 
7 19 37 61 127 271 331 397 547 631 919 1,657 1,801 1,951 2,269 2,437 2,791 3,169 3,571 4,219 4,447 5,167 5,419 6,211 7,057 7,351
8,269 9,241 10,267 11,719 12,097 13,267 13,669 16,651 19,441 19,927 22,447 23,497 24,571 25,117 26,227 27,361 33,391 35,317
42,841 45,757 47,251 49,537 50,311 55,897 59,221 60,919 65,269 70,687 73,477 74,419 75,367 81,181 82,171 87,211 88,237 89,269
92,401 96,661 102,121 103,231 104,347 110,017 112,327 114,661 115,837 126,691 129,169 131,671 135,469 140,617 144,541 145,861
151,201 155,269 163,567 169,219 170,647 176,419 180,811 189,757 200,467 202,021 213,067 231,019 234,361 241,117 246,247 251,431
260,191 263,737 267,307 276,337 279,991 283,669 285,517 292,969 296,731 298,621 310,087 329,677 333,667 337,681 347,821 351,919
360,187 368,551 372,769 374,887 377,011 383,419 387,721 398,581 407,377 423,001 436,627 452,797 459,817 476,407 478,801 493,291
522,919 527,941 553,411 574,219 584,767 590,077 592,741 595,411 603,457 608,851 611,557 619,711 627,919 650,071 658,477 666,937
689,761 692,641 698,419 707,131 733,591 742,519 760,537 769,627 772,669 784,897 791,047 812,761 825,301 837,937 847,477 863,497
879,667 886,177 895,987 909,151 915,769 925,741 929,077 932,419 939,121 952,597 972,991 976,411 986,707 990,151 997,057 1,021,417
1,024,921 1,035,469 1,074,607 1,085,407 1,110,817 1,114,471 1,125,469 1,155,061 1,177,507 1,181,269 1,215,397 1,253,887 1,281,187
1,285,111 1,324,681 1,328,671 1,372,957 1,409,731 1,422,097 1,426,231 1,442,827 1,451,161 1,480,519 1,484,737 1,527,247 1,570,357
output   when using the input of:     -100000 
1,792,617,147,127
Retrieved from "https://rosettacode.org/wiki/Cuban_primes?oldid=287674"