Jump to content

Cubic special primes

From Rosetta Code
Cubic special 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.
Task

n   is smallest prime such that the difference of successive terms are the smallest cubics of positive integers, where     n   <   15000.

11l

Translation of: FreeBASIC
F is_prime(a)
   I a == 2
      R 1B
   I a < 2 | a % 2 == 0
      R 0B
   L(i) (3 .. Int(sqrt(a))).step(2)
      I a % i == 0
         R 0B
   R 1B

V p = 2
V n = 1
print(2, end' ‘ ’)

L p + n ^ 3 < 15000
   I is_prime(p + n ^ 3)
      p += n ^ 3
      n = 1
      print(p, end' ‘ ’)
   E
      n++
Output:
2 3 11 19 83 1811 2027 2243 2251 2467 2531 2539 3539 3547 4547 5059 10891 12619 13619 13627 13691 13907 14419 

Action!

INCLUDE "H6:SIEVE.ACT"

PROC Main()
  DEFINE MAX="14999"
  BYTE ARRAY primes(MAX+1)
  INT i=[2],next,count=[1],n=[1]

  Put(125) PutE() ;clear the screen
  Sieve(primes,MAX+1)
  PrintI(i)
  DO
    next=i+n*n*n
    IF next>=MAX THEN
      EXIT
    FI
    IF primes(next) THEN
      n=1 i=next count==+1
      Put(32) PrintI(i)
    ELSE
      n==+1
    FI
  OD
  PrintF("%E%EThere are %I cubic special primes",count)  
RETURN
Output:

Screenshot from Atari 8-bit computer

2 3 11 19 83 1811 2027 2243 2251 2467 2531 2539 3539 3547 4547 5059 10891 12619 13619 13627 13691 13907 14419

There are 23 cubic special primes

ALGOL 68

BEGIN # find a sequence of primes where the members differ by a cube #
    INT max prime = 15 000;
    # sieve the primes to max prime #
    PR read "primes.incl.a68" PR
    []BOOL prime = PRIMESIEVE max prime;
    # construct a table of cubes, we will need at most the cube root of max prime #
    [ 1 : ENTIER exp( ln( max prime ) / 3 ) ]INT cube;
    FOR i TO UPB cube DO cube[ i ] := i * i * i OD;
    # find the prime sequence #
    print( ( "2" ) );
    INT p := 2;
    WHILE p < max prime DO
        # find a prime that is p + a cube #
        INT q := 0;
        FOR i WHILE q := p + cube[ i ];
                    IF q > max prime THEN FALSE ELSE NOT prime[ q ] FI
        DO SKIP OD;
        IF q <= max prime THEN print( ( " ", whole( q, 0 ) ) ) FI;
        p := q
    OD;
    print( ( newline ) )
END
Output:
2 3 11 19 83 1811 2027 2243 2251 2467 2531 2539 3539 3547 4547 5059 10891 12619 13619 13627 13691 13907 14419

ALGOL W

begin % find a sequence of primes where the members differ by a cube %
    integer MAX_PRIME, MAX_CUBE;
    MAX_PRIME := 15000;
    MAX_CUBE  := entier( exp( ln( MAX_PRIME ) / 3 ) );
    begin
        logical array prime ( 1 :: MAX_PRIME );
        integer array cube  ( 1 :: MAX_CUBE  );
        integer p, q, c;
        % sieve the primes to MAX_PRIME %
        prime( 1 ) := false; prime( 2 ) := true;
        for i := 3 step 2 until MAX_PRIME do prime( i ) := true;
        for i := 4 step 2 until MAX_PRIME do prime( i ) := false;
        for i := 3 step 2 until entier( sqrt( MAX_PRIME ) ) do begin
            integer ii; ii := i + i;
            if prime( i ) then for s := i * i step ii until MAX_PRIME do prime( s ) := false
        end for_i;
        % construct a table of cubes, we will need at most the cube root of max prime %
        for i := 1 until MAX_CUBE do cube( i ) := i * i * i;
        % find the prime sequence %
        writeon( "2" );
        p := 2;
        while p < MAX_PRIME do begin
            % find a prime that is p + a cube %
            c := 1;
            q := p + cube( c );
            while q <= MAX_PRIME and not prime( q ) do begin
                c := c + 1;
                q := p + cube( c )
            end while_not_prime_q ;
            if q <= MAX_PRIME then writeon( i_w := 1, s_w := 0, " ", q );
            p := q
        end while_p_lt_MAX_PRIME
    end;
    write()
end.
Output:
2 3 11 19 83 1811 2027 2243 2251 2467 2531 2539 3539 3547 4547 5059 10891 12619 13619 13627 13691 13907 14419

Arturo

Translation of: FreeBASIC
[p n]: [2 1]
prints -> 2

until -> 15000 =< p + n^3 [
    if? prime? p + n^3 [
        'p + n^3
        n: 1
        prints -> p
    ]
    else -> 'n+1
]
Output:
2 3 11 19 83 1811 2027 2243 2251 2467 2531 2539 3539 3547 4547 5059 10891 12619 13619 13627 13691 13907 14419

AWK

# syntax: GAWK -f CUBIC_SPECIAL_PRIMES.AWK
# converted from FreeBASIC
BEGIN {
    start = p = 2
    stop = 15000
    n = 1
    printf("%5d ",p)
    count = 1
    do {
      if (is_prime(p + n^3)) {
        p += n^3
        n = 1
        printf("%5d%1s",p,++count%10?"":"\n")
      }
      else {
        n++
      }
    } while (p + n^3 <= stop)
    printf("\nCubic special primes %d-%d: %d\n",start,stop,count)
    exit(0)
}
function is_prime(x,  i) {
    if (x <= 1) {
      return(0)
    }
    for (i=2; i<=int(sqrt(x)); i++) {
      if (x % i == 0) {
        return(0)
      }
    }
    return(1)
}
Output:
    2     3    11    19    83  1811  2027  2243  2251  2467
 2531  2539  3539  3547  4547  5059 10891 12619 13619 13627
13691 13907 14419
Cubic special primes 2-15000: 23

C

Translation of: Wren
#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
#include <math.h>
#include <locale.h>

bool *sieve(int limit) {
    int i, p;
    limit++;
    // True denotes composite, false denotes prime.
    bool *c = calloc(limit, sizeof(bool)); // all false by default
    c[0] = true;
    c[1] = true;
    for (i = 4; i < limit; i += 2) c[i] = true;
    p = 3; // Start from 3.
    while (true) {
        int p2 = p * p;
        if (p2 >= limit) break;
        for (i = p2; i < limit; i += 2 * p) c[i] = true;
        while (true) {
            p += 2;
            if (!c[p]) break;
        }
    }
    return c;
}

bool isCube(int n) {
    int s = (int)cbrt((double)n);
    return s*s*s == n;
}

int main() {
    const int limit = 14999;
    int i, j, p, pc = 0, gap = 1, count = 1, lastCubicSpecial = 3;
    const char *fmt = "%'7d %'7d %'6d %'4d\n";
    bool *c = sieve(limit);
    for (i = 0; i < limit; ++i) {
        if (!c[i]) ++pc;
    }
    int *primes = (int *)malloc(pc * sizeof(int));
    for (i = 0, j = 0; i < limit; ++i) {
        if (!c[i]) primes[j++] = i;
    }
    free(c);
    printf("Cubic special primes under 15,000:\n");
    printf(" Prime1  Prime2    Gap  Cbrt\n");
    setlocale(LC_NUMERIC, "");
    printf(fmt, 2, 3, 1, 1);
    for (i = 2; i < pc; ++i) {
        p = primes[i];
        gap = p - lastCubicSpecial;
        if (isCube(gap)) {
            printf(fmt, lastCubicSpecial, p, gap, (int)cbrt((double)gap));
            lastCubicSpecial = p;
            ++count;
        }
    }
    printf("\n%d such primes found.\n", count+1);
    free(primes);
    return 0;
}
Output:
Same as Wren example.

C++

#include <cstdint>
#include <iostream>

bool is_prime(uint32_t number) {
	if ( number % 2 == 0 ) {
		return number == 2;
	}
	int k = 3;
	while ( k * k <= number ) {
		if ( number % k == 0 ) {
			return false;
		}
		k += 2;
	}
	return true;
}

int main() {
	uint32_t p = 2, n = 1;
	std::cout << 2 << " ";

	while ( p + n * n * n < 15'000 ) {
		if ( is_prime(p + n * n * n) ) {
			p += n * n * n;
			n = 1;
			std::cout << p << " ";
		} else {
			n += 1;
		}
	}
}
Output:
2 3 11 19 83 1811 2027 2243 2251 2467 2531 2539 3539 3547 4547 5059 10891 12619 13619 13627 13691 13907 14419 

Delphi

Works with: Delphi version 6.0

This program makes extensive use of the Delphi Prime-Generator Object

procedure CubicSpecialPrimes(Memo: TMemo);
const Limit = 15000;
var Sieve: TPrimeSieve;
var N,I,PP, Count: integer;
var S: string;
begin
Sieve:=TPrimeSieve.Create;
try
{Get all primes up to limit}
Sieve.Intialize(Limit);
N:=1;
Count:=0;
I:=1;
S:='';
while true do
	begin
	{Find first cube increment that is prime}
	PP:=N +I*I*I;
	if PP>Limit then break;
	if Sieve.Flags[PP] then
		begin
		Inc(Count);
		S:=S+Format('%6D',[PP]);
		if (Count mod 5) = 0 then S:=S+CRLF;
		{Step to next cube position}
		I:=1;
		N:=PP;
		end
	else Inc(I);
	end;
Memo.Lines.Add(Format('There are %d cubic special primes',[count]));
Memo.Lines.Add(S);
finally Sieve.Free; end;
end;
Output:
There are 23 cubic special primes
     2     3    11    19    83
  1811  2027  2243  2251  2467
  2531  2539  3539  3547  4547
  5059 10891 12619 13619 13627
 13691 13907 14419
Elapsed Time: 2.178 ms.

EasyLang

Translation of: 11l
fastfunc isprim num .
   i = 2
   while i <= sqrt num
      if num mod i = 0
         return 0
      .
      i += 1
   .
   return 1
.
p = 2
n = 1
write 2 & " "
repeat
   h = p + n * n * n
   until h >= 15000
   if isprim h = 1
      p = h
      n = 1
      write p & " "
   else
      n += 1
   .
.
Output:
2 3 11 19 83 1811 2027 2243 2251 2467 2531 2539 3539 3547 4547 5059 10891 12619 13619 13627 13691 13907 14419 

F#

// Cubic Special Primes: Nigel Galloway. March 30th., 2021
let fN=let N=[for n in [0..25]->n*n*n] in let mutable n=2 in (fun g->match List.contains(g-n)N with true->n<-g; true |_->false)
primes32()|>Seq.takeWhile((>)16000)|>Seq.filter fN|>Seq.iter(printf "%d "); printfn ""
Output:
2 3 11 19 83 1811 2027 2243 2251 2467 2531 2539 3539 3547 4547 5059 10891 12619 13619 13627 13691 13907 14419

Factor

Works with: Factor version 0.98
USING: fry io kernel lists lists.lazy math math.functions
math.primes prettyprint ;

2 [ 1 lfrom swap '[ 3 ^ _ + ] lmap-lazy [ prime? ] lfilter car ]
lfrom-by [ 15000 < ] lwhile [ pprint bl ] leach nl
Output:
2 3 11 19 83 1811 2027 2243 2251 2467 2531 2539 3539 3547 4547 5059 10891 12619 13619 13627 13691 13907 14419

FreeBASIC

#include once "isprime.bas"

dim as integer p = 2, n = 1

print 2;" ";

do
    if isprime(p + n^3) then
        p += n^3
        n = 1
        print p;" ";
    else
        n += 1
    end if
loop until p + n^3 >= 15000
print
end
Output:

2 3 11 19 83 1811 2027 2243 2251 2467 2531 2539 3539 3547 4547 5059 10891 12619 13619 13627 13691 13907 14419


FutureBasic

local fn IsPrime( n as NSUInteger ) as BOOL
—
  BOOL       isPrime = YES
  NSUInteger i
  
  if n < 2        then exit fn = NO
  if n = 2        then exit fn = YES
  if n mod 2 == 0 then exit fn = NO
  for i = 3 to int(n^.5) step 2
    if n mod i == 0 then exit fn = NO
  next
end fn = isPrime

int p, n, counter
p = 2
n = 1
counter = 2

printf @"%5d \b", 2

do
  if fn IsPrime ( p + n^3 )
    p += n^3
    n = 1
    printf @"%5d \b", p
    if counter mod 8 == 0 then print
    counter++
  else
    n++
  end if
until ( p + n^3 >= 15000 )

HandleEvents
Output:
    2     3    11    19    83  1811  2027  2243 
 2251  2467  2531  2539  3539  3547  4547  5059 
10891 12619 13619 13627 13691 13907 14419 

Go

Translation of: Wren
package main

import (
    "fmt"
    "math"
)

func sieve(limit int) []bool {
    limit++
    // True denotes composite, false denotes prime.
    c := make([]bool, limit) // all false by default
    c[0] = true
    c[1] = true
    // no need to bother with even numbers over 2 for this task
    p := 3 // Start from 3.
    for {
        p2 := p * p
        if p2 >= limit {
            break
        }
        for i := p2; i < limit; i += 2 * p {
            c[i] = true
        }
        for {
            p += 2
            if !c[p] {
                break
            }
        }
    }
    return c
}

func isCube(n int) bool {
    s := int(math.Cbrt(float64(n)))
    return s*s*s == n
}

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

func main() {
    c := sieve(14999)
    fmt.Println("Cubic special primes under 15,000:")
    fmt.Println(" Prime1  Prime2    Gap  Cbrt")
    lastCubicSpecial := 3
    gap := 1
    count := 1
    fmt.Printf("%7d %7d %6d %4d\n", 2, 3, 1, 1)
    for i := 5; i < 15000; i += 2 {
        if c[i] {
            continue
        }
        gap = i - lastCubicSpecial
        if isCube(gap) {
            cbrt := int(math.Cbrt(float64(gap)))
            fmt.Printf("%7s %7s %6s %4d\n", commas(lastCubicSpecial), commas(i), commas(gap), cbrt)
            lastCubicSpecial = i
            count++
        }
    }
    fmt.Println("\n", count+1, "such primes found.")
}
Output:
Same as Wren example.

J

cubes=. 3 ^~ 1 , 10 + i: 8j8
   
nextp=. ({.@#~ 1&p:)@(+&cubes)
   
nextp^:(14e3&>)^:a: 2
Output:
2 3 11 19 83 1811 2027 2243 2251 2467 2531 2539 3539 3547 4547 5059 10891 12619 13619 13627 13691 13907 14419

Java

public final class CubicSpecialPrimes {

	public static void main(String[] args) {
		int p = 2;
		int n = 1;
		System.out.print(2 + " ");

		while ( p + n * n * n < 15_000 ) {
		    if ( isPrime(p + n * n * n) ) {
		        p += n * n * n;
		        n = 1;
		        System.out.print(p + " ");
		    } else {
		        n += 1;
		    }
		} 
	}

	private static boolean isPrime(int number) {
	    if ( number % 2 == 0 ) {
	    	return number == 2;
	    }
	    
	    int k = 3;
	    while ( k * k <= number ) {
	    	if ( number % k == 0 ) {
	    		return false;
	    	}
	    	k += 2;
	    }
	    return true;
	}
	
}
Output:
2 3 11 19 83 1811 2027 2243 2251 2467 2531 2539 3539 3547 4547 5059 10891 12619 13619 13627 13691 13907 14419 

jq

Adapted from Julia

Works with: jq

Works with gojq, the Go implementation of jq

For the definition of `is_prime` used here, see https://rosettacode.org/wiki/Additive_primes

# Output: an array
def cubic_special_primes:
  . as $N
  | sqrt as $maxidx
  | {cprimes: [2], q: 0}
  | until( .cprimes[-1] >= $N or .q >= $N;
      label $out
      | foreach range(1; $maxidx + 1) as $i (.;
          .q = (.cprimes[-1] + ($i * $i * $i))
          | if .q >= $N
	    then .emit = true
            elif .q | is_prime
            then .cprimes = .cprimes + [.q]
            | .emit = true
            else .
	    end;
	    select(.emit)) | {cprimes, q}, break $out )
  | .cprimes ;	 
 
15000
| ("Cubic special primes < \(.):",
    cubic_special_primes)
Output:
[2,3,11,19,83,1811,2027,2243,2251,2467,2531,2539,3539,3547,4547,5059,10891,12619,13619,13627,13691,13907,14419]

Julia

using Primes

function cubicspecialprimes(N = 15000)
    pmask = primesmask(1, N)
    cprimes, maxidx = [2], isqrt(N)
    while (n = cprimes[end]) < N
        for i in 1:maxidx
            q = n + i * i * i
            if  q > N
                return cprimes
            elseif pmask[q]  # got next qprime
                push!(cprimes, q)
                break
            end
        end
    end
end

println("Cubic special primes < 15000:")
foreach(p -> print(rpad(p[2], 6), p[1] % 10 == 0 ? "\n" : ""), enumerate(cubicspecialprimes()))
Output:
Cubic special primes < 16000:
2     3     11    19    83    1811  2027  2243  2251  2467  
2531  2539  3539  3547  4547  5059  10891 12619 13619 13627 
13691 13907 14419

Lua

do  -- find a sequence of primes where the members differ by a cube

    local maxPrime = 15000
    -- sieve the primes to maxPrime
    local isPrime = {}
    for i = 1, maxPrime do isPrime[ i ] = i % 2 ~= 0 end
    isPrime[ 1 ]  = false
    isPrime[ 2 ]  = true
    for s = 3, math.floor( math.sqrt( maxPrime ) ), 2 do
        if isPrime[ s ] then
            for p = s * s, maxPrime, s do isPrime[ p ] = false end
        end
    end

    -- construct a table of cubes, we will need at most the cube root of maxPrime
    local cube = {}
    for i = 1, math.floor( ( maxPrime ^ ( 1 / 3 ) ) ) do cube[ i ] = i * i * i end

    -- find the prime sequence
    io.write( "2" )
    local p = 2
    while p < maxPrime do
        -- find a prime that is p + a cube
        local q, i = 0, 1
        repeat
            q, i = p + cube[ i ], i + 1
        until q > maxPrime or isPrime[ q ]
        if q <= maxPrime then io.write( " ", q ) end
        p = q
    end
    io.write( "\n" )
end
Output:
2 3 11 19 83 1811 2027 2243 2251 2467 2531 2539 3539 3547 4547 5059 10891 12619 13619 13627 13691 13907 14419

Mathematica /Wolfram Language

start = {2};
Do[
 If[PrimeQ[i],
  last = Last[start];
  If[IntegerQ[Surd[i - last, 3]],
   AppendTo[start, i]
   ]
  ]
 ,
 {i, 3, 15000}
 ]
start
Output:
{2,3,11,19,83,1811,2027,2243,2251,2467,2531,2539,3539,3547,4547,5059,10891,12619,13619,13627,13691,13907,14419}

Nim

Translation of: Go
import math, strformat

func sieve(limit: Positive): seq[bool] =
  # True denotes composite, false denotes prime.
  result = newSeq[bool](limit + 1)   # All false by default.
  result[0] = true
  result[1] = true
  # No need to bother with even numbers over 2 for this task.
  var p = 3
  while true:
    let p2 = p * p
    if p2 > limit: break
    for i in countup(p2, limit, 2 * p):
      result[i] = true
    while true:
      inc p, 2
      if not result[p]: break

func isCube(n: int): bool =
  let s = cbrt(n.toFloat).int
  result = s * s * s == n

let c = sieve(14999)
echo "Cubic special primes under 15_000:"
echo " Prime1  Prime2    Gap  Cbrt"
var lastCubicSpecial = 3
var count = 1
echo &"{2:7} {3:7} {1:6} {1:4}"
for n in countup(5, 14999, 2):
  if c[n]: continue
  let gap = n - lastCubicSpecial
  if gap.isCube:
    let gapCbrt = cbrt(gap.toFloat).int
    echo &"{lastCubicSpecial:7} {n:7} {gap:6} {gapCbrt:4}"
    lastCubicSpecial = n
    inc count
echo &"\n{count + 1} such primes found."
Output:
Cubic special primes under 15_000:
 Prime1  Prime2    Gap  Cbrt
      2       3      1    1
      3      11      8    2
     11      19      8    2
     19      83     64    4
     83    1811   1728   12
   1811    2027    216    6
   2027    2243    216    6
   2243    2251      8    2
   2251    2467    216    6
   2467    2531     64    4
   2531    2539      8    2
   2539    3539   1000   10
   3539    3547      8    2
   3547    4547   1000   10
   4547    5059    512    8
   5059   10891   5832   18
  10891   12619   1728   12
  12619   13619   1000   10
  13619   13627      8    2
  13627   13691     64    4
  13691   13907    216    6
  13907   14419    512    8

23 such primes found.

Perl

Library: ntheory
use strict;
use warnings;
use feature 'say';
use  ntheory 'is_prime';

my @sp = my $previous = 2;
do {
    my($next,$n);
    while () { last if is_prime( $next = $previous + ++$n**3 ) }
    push @sp, $next;
    $previous = $next;
} until $sp[-1] >= 15000;

pop @sp and say join ' ', @sp;
Output:
2 3 11 19 83 1811 2027 2243 2251 2467 2531 2539 3539 3547 4547 5059 10891 12619 13619 13627 13691 13907 14419

Phix

See Quadrat_Special_Primes#Phix


Python

#!/usr/bin/python

def isPrime(n):
    for i in range(2, int(n**0.5) + 1):
        if n % i == 0:
            return False        
    return True

if __name__ == '__main__':
    p = 2
    n = 1

    print("2",end = " ")
    while True:
        if isPrime(p + n**3):
            p += n**3
            n = 1
            print(p,end = " ")
        else:
            n += 1
        if p + n**3 >= 15000:
            break
Output:
2  3  11  19  83  1811  2027  2243  2251  2467  2531  2539  3539  3547  4547  5059  10891  12619  13619  13627  13691  13907  14419


Raku

A two character difference from the Quadrat Special Primes entry. (And it could have been one.)

my @sqp = 2, -> $previous {
    my $next;
    for (1..∞).map: *³ {
        $next = $previous + $_;
        last if $next.is-prime;
    }
    $next
} … *;

say "{+$_} matching numbers:\n", $_».fmt('%5d').batch(7).join: "\n" given
    @sqp[^(@sqp.first: * > 15000, :k)];
Output:
23 matching numbers:
    2     3    11    19    83  1811  2027
 2243  2251  2467  2531  2539  3539  3547
 4547  5059 10891 12619 13619 13627 13691
13907 14419

REXX

Version 1 Inline code

/*REXX program finds the smallest primes such that the difference of successive terms   */
/*───────────────────────────────────────────────────── are the smallest cubic numbers. */
parse arg hi cols .                              /*obtain optional argument from the CL.*/
if   hi=='' |   hi==","  then   hi= 15000        /* "      "         "   "   "     "    */
if cols=='' | cols==","  then cols=    10        /* "      "         "   "   "     "    */
call genP                                        /*build array of semaphores for primes.*/
w= 10                                            /*width of a number in any column.     */
                    title= 'the smallest primes  < '    commas(hi)     " such that the" ,
                           'difference of successive terma are the smallest cubic numbers'
if cols>0 then say ' index │'center(title,  1 + cols*(w+1)     )    /*display the title.*/
if cols>0 then say '───────┼'center(""   ,  1 + cols*(w+1), '─')    /*   "     "    sep.*/
found= 0;             idx= 1                     /*initialize number of  cbp  and index.*/
op= 1
$=                                               /*a list of  cubic  primes  (so far).  */
     do j=0  by 0
                      do k=1  until !.np
                      np= op + k**3              /*find the next cube plus the oldPrime.*/
                      if np>= hi  then leave j   /*Is newPrime too big?  Then quit.     */
                      end   /*k*/
     found= found + 1                            /*bump number of primes of this type.  */
     op= np                                      /*assign the newPrime  to the  oldPrime*/
     if cols<0           then iterate            /*Build the list  (to be shown later)? */
     c= commas(np)                               /*maybe add commas to the number.      */
     $= $ right(c, max(w, length(c) ) )          /*add a cubic prime──►list, allow big#.*/
     if found//cols\==0  then iterate            /*have we populated a line of output?  */
     say center(idx, 7)'│'  substr($, 2);   $=   /*display what we have so far  (cols). */
     idx= idx + cols                             /*bump the  index  count for the output*/
     end   /*j*/

if $\==''  then say center(idx, 7)"│"  substr($, 2)  /*possible display residual output.*/
if cols>0 then say '───────┴'center(""   ,  1 + cols*(w+1), '─')    /*display foot sep. */
say
say 'Found '       commas(found)         " of "       title
exit 0                                           /*stick a fork in it,  we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
commas: parse arg ?;  do jc=length(?)-3  to 1  by -3; ?=insert(',', ?, jc); end;  return ?
/*──────────────────────────────────────────────────────────────────────────────────────*/
genP: !.= 0                                      /*placeholders for primes (semaphores).*/
      @.1=2;  @.2=3;  @.3=5;  @.4=7;  @.5=11     /*define some low primes.              */
      !.2=1;  !.3=1;  !.5=1;  !.7=1;  !.11=1     /*   "     "   "    "     flags.       */
                        #=5;     s.#= @.# **2    /*number of primes so far;     prime². */
                                                 /* [↓]  generate more  primes  ≤  high.*/
        do j=@.#+2  by 2  to hi                  /*find odd primes from here on.        */
        parse var j '' -1 _; if     _==5  then iterate  /*J divisible by 5?  (right dig)*/
                             if j// 3==0  then iterate  /*"     "      " 3?             */
                             if j// 7==0  then iterate  /*"     "      " 7?             */
               do k=5  while s.k<=j              /* [↓]  divide by the known odd primes.*/
               if j // @.k == 0  then iterate j  /*Is  J ÷ X?  Then not prime.     ___  */
               end   /*k*/                       /* [↑]  only process numbers  ≤  √ J   */
        #= #+1;    @.#= j;    s.#= j*j;   !.j= 1 /*bump # of Ps; assign next P;  P²; P# */
        end          /*j*/;   return
output   when using the default inputs:
 index │  the smallest primes  <  15,000  such that the difference of successive terma are the smallest cubic numbers
───────┼───────────────────────────────────────────────────────────────────────────────────────────────────────────────
   1   │          2          3         11         19         83      1,811      2,027      2,243      2,251      2,467
  11   │      2,531      2,539      3,539      3,547      4,547      5,059     10,891     12,619     13,619     13,627
  21   │     13,691     13,907     14,419
───────┴───────────────────────────────────────────────────────────────────────────────────────────────────────────────

Found  23  of  the smallest primes  <  15,000  such that the difference of successive terma are the smallest cubic numbers
0.040 seconds

Version 2 Using REXX libraries and simplified

Libraries: How to use
Library: Sequences
Library: Functions
Library: Settings
Library: Abend

Below program shows how simple this task can be. Compare it with Version 1, which essentially uses the same technique! Procedure Xpon() (in Functions) returns the exponent of a number, used to control digits and layout. Procedure Primes() (in Sequences, using a sieve) returns an array with primes and a sparse array with numbers flagged as prime; the latter one is needed. So no functions like IsPrime? or IsCube? used! Some other entries employ the same fast solution. Version 2 works up to t = 100 million within acceptable time.

call Time('r')
say 'Cubic special primes - Using REXX libraries'
parse version version; say version; say
t = 150000; d = Xpon(t)+2; numeric digits d
call Primes t
call ShowCubics t,d
say
say Format(Time('e'),,3) 'seconds'
exit

ShowCubics:
procedure expose prim.
arg x,y
call Charout ,Right(2,y)
a = 2; b = 3; k = 1; n = 1
do while b <= x
   b = a + k*k*k
   if prim.flag.b = 1 then do
      n = n+1
      call Charout ,Right(b,y)
      if n//10 = 0 then
         say
      a = b; k = 1
   end
   k = k+1
end
say; say
say n 'cubic special primes found below' x
return

include Sequences
include Functions
Output:
REXX-Regina_3.9.6(MT) 5.00 29 Apr 2024
Cubic special primes

     2     3    11    19    83  1811  2027  2243  2251  2467
  2531  2539  3539  3547  4547  5059 10891 12619 13619 13627
 13691 13907 14419

23 cubic special primes found below 15000

0.006 seconds

Ring

Also see Quadrat_Special_Primes#Ring

load "stdlib.ring"

see "working..." + nl

Primes = []
limit1 = 50
oldPrime = 2
add(Primes,2)

for n = 1 to limit1
    nextPrime = oldPrime + pow(n,3)
    if isprime(nextPrime)
       n = 1
       add(Primes,nextPrime)
       oldPrime = nextPrime
    else
       nextPrime = nextPrime - oldPrime
    ok
next

see "prime1 prime2 Gap Cbrt" + nl
for n = 1 to Len(Primes)-1
    diff = Primes[n+1] - Primes[n]
    for m = 1 to diff
        if pow(m,3) = diff
           cbrt = m
           exit
        ok
    next
    see ""+ Primes[n] + "      " + Primes[n+1] + "    " + diff + "     " + cbrt + nl
next   
 
see "Found " + Len(Primes) + " of the smallest primes < 15,000  such that the difference of successive terma are the smallest cubic numbers" + nl  
 
see "done..." + nl
Output:
working...
prime1 prime2 Gap Cbrt
2      3    1     1
3      11    8     2
11      19    8     2
19      83    64     4
83      1811    1728     12
1811      2027    216     6
2027      2243    216     6
2243      2251    8     2
2251      2467    216     6
2467      2531    64     4
2531      2539    8     2
2539      3539    1000     10
3539      3547    8     2
3547      4547    1000     10
4547      5059    512     8
5059      10891    5832     18
10891      12619    1728     12
12619      13619    1000     10
13619      13627    8     2
13627      13691    64     4
13691      13907    216     6
13907      14419    512     8
Found 23 of the smallest primes < 15,000  such that the difference of successive terma are the smallest cubic numbers
done...

RPL

Works with: HP version 49
« { } 2 1 
  DO DUP2 3 ^ +
     IF DUP ISPRIME?    
     THEN 4 ROLL 4 ROLL + ROT DROP SWAP 1 
     ELSE DROP 1 + END
  UNTIL OVER 15000 > END
  DROP2
» 'TASK' STO
Output:
1: { 2 3 11 19 83 1811 2027 2243 2251 2467 2531 2539 3539 3547 4547 5059 10891 12619 13619 13627 13691 13907 14419 }

Ruby

require 'prime'

res = [2]

until res.last > 15000 do
  res << (1..).detect{|n| (res.last + n**3).prime? } ** 3 + res.last
end

puts res[..-2].join(" ")
Output:
2 3 11 19 83 1811 2027 2243 2251 2467 2531 2539 3539 3547 4547 5059 10891 12619 13619 13627 13691 13907 14419

Rust

Translation of: Python
use num_prime::{nt_funcs, Primality, PrimalityTestConfig};


const LIMIT: usize = 15_000;

fn main() {
    let mut p = 2_usize;
    let mut n = 1;

    print!("2 ");
    loop {
        let new_p = p + n * n * n;
        if nt_funcs::is_prime(&new_p, Some(PrimalityTestConfig::strict())) == Primality::Yes {
            p = new_p;
            n = 1;
            print!("{} ", p);
        } else {
            n += 1;
            if new_p >= LIMIT {
                break;
            }
        }
    }
}
Output:
2 3 11 19 83 1811 2027 2243 2251 2467 2531 2539 3539 3547 4547 5059 10891 12619 13619 13627 13691 13907 14419

Sidef

func cubic_primes(callback) {

    var prev = 2
    callback(prev)

    loop {
        var curr = (1..Inf -> lazy.map { prev + _**3 }.first { .is_prime })
        callback(curr)
        prev = curr
    }
}

say gather {
    cubic_primes({|k|
        break if (k >= 15000)
        take(k)
    })
}
Output:
[2, 3, 11, 19, 83, 1811, 2027, 2243, 2251, 2467, 2531, 2539, 3539, 3547, 4547, 5059, 10891, 12619, 13619, 13627, 13691, 13907, 14419]

Wren

Library: Wren-math
Library: Wren-fmt
import "./math" for Int
import "./fmt" for Fmt

var primes = Int.primeSieve(14999)
System.print("Cubic special primes under 15,000:")
System.print(" Prime1  Prime2    Gap  Cbrt")
var lastCubicSpecial = 3
var gap = 1
var count = 1
Fmt.print("$,7d $,7d $,6d $4d", 2, 3, 1, 1)
for (p in primes.skip(2)) {
    gap = p - lastCubicSpecial
    if (Int.isCube(gap)) {
        Fmt.print("$,7d $,7d $,6d $4d", lastCubicSpecial, p, gap, gap.cbrt.truncate)
        lastCubicSpecial = p
        count = count + 1
    }
}
System.print("\n%(count+1) such primes found.")
Output:
Cubic special primes under 15,000:
 Prime1  Prime2    Gap  Cbrt
      2       3      1    1
      3      11      8    2
     11      19      8    2
     19      83     64    4
     83   1,811  1,728   12
  1,811   2,027    216    6
  2,027   2,243    216    6
  2,243   2,251      8    2
  2,251   2,467    216    6
  2,467   2,531     64    4
  2,531   2,539      8    2
  2,539   3,539  1,000   10
  3,539   3,547      8    2
  3,547   4,547  1,000   10
  4,547   5,059    512    8
  5,059  10,891  5,832   18
 10,891  12,619  1,728   12
 12,619  13,619  1,000   10
 13,619  13,627      8    2
 13,627  13,691     64    4
 13,691  13,907    216    6
 13,907  14,419    512    8

23 such primes found.

XPL0

func IsPrime(N);        \Return 'true' if N is prime
int  N, I;
[if N <= 2 then return N = 2;
if (N&1) = 0 then \even >2\ return false;
for I:= 3 to sqrt(N) do
    [if rem(N/I) = 0 then return false;
    I:= I+1;
    ];
return true;
];

int N, C, T;
[N:= 2;
repeat  C:= 1;
        loop    [T:= N + C*C*C;
                if IsPrime(T) then
                        [IntOut(0, N);
                        ChOut(0, ^ );
                        N:= T;
                        quit;
                        ]
                else    C:= C+1;
                ];
until   N > 15000;
]
Output:
2 3 11 19 83 1811 2027 2243 2251 2467 2531 2539 3539 3547 4547 5059 10891 12619 13619 13627 13691 13907 14419 
Cookies help us deliver our services. By using our services, you agree to our use of cookies.