De Polignac numbers: Difference between revisions

From Rosetta Code
Content added Content deleted
(Added Algol 68 and Algol W)
(Added PL/M)
Line 319: Line 319:
One thousandth: 31,941
One thousandth: 31,941
Ten thousandth: 273,421
Ten thousandth: 273,421
</pre>

=={{header|PL/M}}==
Based on the ALGOL 68 sample.
{{works with|8080 PL/M Compiler}} ... under CP/M (or an emulator)
<syntaxhighlight lang="plm">
100H: /* FIND SOME DE POLIGNAC NUMBERS - POSITIVE ODD NUMBERS THAT CAN'T BE */
/* WRITTEN AS P + 2**N FOR SOME PRIME P AND INTEGER N */

DECLARE FALSE LITERALLY '0', TRUE LITERALLY '0FFH';

/* CP/M SYSTEM CALL AND I/O ROUTINES */
BDOS: PROCEDURE( FN, ARG ); DECLARE FN BYTE, ARG ADDRESS; GOTO 5; END;
PR$CHAR: PROCEDURE( C ); DECLARE C BYTE; CALL BDOS( 2, C ); END;
PR$STRING: PROCEDURE( S ); DECLARE S ADDRESS; CALL BDOS( 9, S ); END;
PR$NL: PROCEDURE; CALL PR$CHAR( 0DH ); CALL PR$CHAR( 0AH ); END;
PR$NUMBER: PROCEDURE( N ); /* PRINTS A NUMBER IN THE MINIMUN FIELD WIDTH */
DECLARE N ADDRESS;
DECLARE V ADDRESS, N$STR ( 6 )BYTE, W BYTE;
V = N;
W = LAST( N$STR );
N$STR( W ) = '$';
N$STR( W := W - 1 ) = '0' + ( V MOD 10 );
DO WHILE( ( V := V / 10 ) > 0 );
N$STR( W := W - 1 ) = '0' + ( V MOD 10 );
END;
CALL PR$STRING( .N$STR( W ) );
END PR$NUMBER;
/* END SYSTEM CALL AND I/O ROUTINES */

DECLARE MAX$DP LITERALLY '4000', /* MAXIMUM NUMBER TO CONSIDER */
MAX$DP$PLUS$1 LITERALLY '4001'; /* MAX$DP + 1 FOR ARRAY BOUNDS */

/* SIEVE THE PRIMES TO MAX$DP */
DECLARE PRIME ( MAX$DP$PLUS$1 )BYTE;
DO;
DECLARE ( I, S ) ADDRESS;
PRIME( 0 ), PRIME( 1 ) = FALSE;
PRIME( 2 ) = TRUE;
DO I = 3 TO LAST( PRIME ) BY 2; PRIME( I ) = TRUE; END;
DO I = 4 TO LAST( PRIME ) BY 2; PRIME( I ) = FALSE; END;
DO I = 3 TO LAST( PRIME ) / 2 BY 2;
IF PRIME( I ) THEN DO;
DO S = I + I TO LAST( PRIME ) BY I; PRIME( S ) = FALSE; END;
END;
END;
END;

/* TABLE OF POWERS OF 2 UP TO MAX$DP */
DECLARE POWERS$OF$2 ( 12 )ADDRESS
INITIAL( 1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048 );

/* DISPLAYS A DE POLIGNAC NUMBER, IF NECESSARY */
SHOW$DE$POLIGNAC: PROCEDURE( DP$NUMBER );
DECLARE ( DP$NUMBER ) ADDRESS;
CALL PR$CHAR( ' ' );
IF DP$NUMBER < 10 THEN CALL PR$CHAR( ' ' );
IF DP$NUMBER < 100 THEN CALL PR$CHAR( ' ' );
IF DP$NUMBER < 1000 THEN CALL PR$CHAR( ' ' );
CALL PR$NUMBER( DP$NUMBER );
END SHOW$DE$POLIGNAC;

DECLARE ( I, P, COUNT ) ADDRESS;
DECLARE FOUND BYTE;

/* THE NUMBERS MUST BE ODD AND NOT OF THE FORM P + 2**N */
/* EITHER P IS ODD AND 2**N IS EVEN AND HENCE N > 0 AND P > 2 */
/* OR 2**N IS ODD AND P IS EVEN AND HENCE N = 0 AND P = 2 */

/* N = 0, P = 2 - THE ONLY POSSIBILITY IS 3 */
DO I = 1 TO 3 BY 2;
P = 2;
IF P + 1 <> I THEN DO;
COUNT = COUNT + 1;
CALL SHOW$DE$POLIGNAC( I );
END;
END;
/* N > 0, P > 2 */
I = 3;
DO WHILE I < MAX$DP AND COUNT < 50;
I = I + 2;
FOUND = FALSE;
P = 1;
DO WHILE NOT FOUND
AND P <= LAST( POWERS$OF$2 )
AND I > POWERS$OF$2( P );
FOUND = PRIME( I - POWERS$OF$2( P ) );
P = P + 1;
END;
IF NOT FOUND THEN DO;
CALL SHOW$DE$POLIGNAC( I );
IF ( COUNT := COUNT + 1 ) MOD 10 = 0 THEN CALL PR$NL;
END;
END;

EOF
</syntaxhighlight>
{{out}}
<pre>
1 127 149 251 331 337 373 509 599 701
757 809 877 905 907 959 977 997 1019 1087
1199 1207 1211 1243 1259 1271 1477 1529 1541 1549
1589 1597 1619 1649 1657 1719 1759 1777 1783 1807
1829 1859 1867 1927 1969 1973 1985 2171 2203 2213
</pre>
</pre>


Revision as of 21:02, 27 September 2022

De Polignac numbers 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.

Alphonse de Polignac, a French mathematician in the 1800s, conjectured that every positive odd integer could be formed from the sum of a power of 2 and a prime number.

He was subsequently proved incorrect.

The numbers that fail this condition are now known as de Polignac numbers.

Technically 1 is a de Polignac number, as there is no prime and power of 2 that sum to 1. De Polignac was aware but thought that 1 was a special case. However.

127 is also fails that condition, as there is no prime and power of 2 that sum to 127.

As it turns out, de Polignac numbers are not uncommon, in fact, there are an infinite number of them.


Task
  • Find and display the first fifty de Polignac numbers.


Stretch
  • Find and display the one thousandth de Polignac number.
  • Find and display the ten thousandth de Polignac number.


See also


Action!

Translation of: PL/M

which is based on the ALGOL 68 sample.

Library: Action! Sieve of Eratosthenes
;;; Find some De Polignac numbers: positive odd numbers that can't be
;;;      written as p + 2**n for some prime p and integer n

INCLUDE "H6:SIEVE.ACT"

PROC showDePolignac( CARD dpNumber )
  Put(' )
  IF dpNumber <   10 THEN Put(' ) FI
  IF dpNumber <  100 THEN Put(' ) FI
  IF dpNumber < 1000 THEN Put(' ) FI
  PrintC( dpNumber )
RETURN

PROC Main()
  DEFINE MAX_DP = "4000", MAX_POWER_OF_2 = "11"

  BYTE ARRAY primes(MAX_DP+1)
  CARD ARRAY powersOf2(MAX_POWER_OF_2+1) =
                 [1 2 4 8 16 32 64 128 256 512 1024 2048]
  CARD i, p, count
  BYTE found
 
  Sieve(primes,MAX_DP+1)

  ; the numbers must be odd and not of the form p + 2**n
  ; either p is odd and 2**n is even and hence n > 0 and p > 2
  ;     or 2**n is odd and p is even and hence n = 0 and p = 2

  ; n = 0, p = 2 - the only possibility is 3
  FOR i = 1 TO 3 STEP 2 DO
    p = 2
    IF p + 1 <> i THEN
      count ==+ 1
      showDePolignac( i )
    FI
  OD
  ; n > 0, p > 2
  i = 3
  WHILE i < MAX_DP AND count < 50 DO
    i ==+ 2
    found = 0
    p = 1
    WHILE found = 0 AND p <= MAX_POWER_OF_2 AND i > powersOf2( p ) DO
       found = primes( i - powersOf2( p ) )
       p ==+ 1
    OD
    IF found = 0 THEN
      count ==+ 1
      showDePolignac( i )
      IF count MOD 10 = 0 THEN PutE() FI
    FI
  OD

RETURN
Output:
    1  127  149  251  331  337  373  509  599  701
  757  809  877  905  907  959  977  997 1019 1087
 1199 1207 1211 1243 1259 1271 1477 1529 1541 1549
 1589 1597 1619 1649 1657 1719 1759 1777 1783 1807
 1829 1859 1867 1927 1969 1973 1985 2171 2203 2213

ALGOL 68

BEGIN # find some De Polignac Numbers - positive odd numbers that can't be    #
      # written as p + 2^n for some prime p and integer n                     #
    INT max number = 500 000; # maximum number we will consider               #
    # sieve the primes to max number                                          #
    [ 0 : max number ]BOOL prime;
    prime[ 0 ] := prime[ 1 ] := FALSE;
    prime[ 2 ] := TRUE;
    FOR i FROM 3 BY 2 TO UPB prime DO prime[ i ] := TRUE  OD;
    FOR i FROM 4 BY 2 TO UPB prime DO prime[ i ] := FALSE OD;
    FOR i FROM 3 BY 2 TO ENTIER sqrt( UPB prime ) DO
        IF prime[ i ] THEN
            FOR s FROM i * i BY i + i TO UPB prime DO prime[ s ] := FALSE OD
        FI
    OD;
    # table of powers of 2 greater than 2^0 ( up to around 2 000 000 )        #
    #       increase the table size if max number > 2 000 000                 #
    [ 1 : 20 ]INT powers of 2;
    BEGIN
        INT p2 := 1;
        FOR i TO UPB powers of 2 DO
            powers of 2[ i ] := ( p2 *:= 2 )
        OD
    END;
    # the numbers must be odd and not of the form p + 2^n                     #
    # either p is odd and 2^n is even and hence n > 0 and p > 2               #
    #     or 2^n is odd and p is even and hence n = 0 and p = 2               #
    INT dp count := 0;
    # n = 0, p = 2 - the only possibility is 3                                #
    FOR i BY 2 TO 3 DO
        INT p = 2;
        IF p + 1 /= i THEN
            print( ( whole( i, -5 ) ) );
            dp count +:= 1
        FI
    OD;
    # n > 0, p > 2                                                            #
    FOR i FROM 5 BY 2 TO max number DO
        BOOL found := FALSE;
        FOR p TO UPB powers of 2 WHILE NOT found AND i > powers of 2[ p ] DO
            found := prime[ i - powers of 2[ p ] ]
        OD;
        IF NOT found THEN
            IF ( dp count +:= 1 ) <= 50 THEN
                print( ( whole( i, -5 ) ) );
                IF dp count MOD 10 = 0 THEN print( ( newline ) ) FI
            ELIF dp count = 1 000 OR dp count = 10 000 THEN
                print( ( "The ", whole( dp count, -5 )
                       , "th De Polignac number is ", whole( i, -7 )
                       , newline
                       )
                     )
            FI
        FI
    OD;
    print( ( "Found ", whole( dp count, 0 )
           , " De Polignac numbers up to ", whole( max number, 0 )
           , newline
           )
         )
END
Output:
    1  127  149  251  331  337  373  509  599  701
  757  809  877  905  907  959  977  997 1019 1087
 1199 1207 1211 1243 1259 1271 1477 1529 1541 1549
 1589 1597 1619 1649 1657 1719 1759 1777 1783 1807
 1829 1859 1867 1927 1969 1973 1985 2171 2203 2213
The  1000th De Polignac number is   31941
The 10000th De Polignac number is  273421
Found 19075 De Polignac numbers up to 500000

ALGOL W

Translation of: ALGOL 68
begin % find some De Polignac numbers - positive odd numbers that can't be    %
      % written as p + 2^n for some prime p and integer n                     %
    integer MAX_NUMBER;                     % maximum number we will consider %
    integer MAX_POWER;                      % maximum powerOf2 < MAX_NUMBER   %
    MAX_NUMBER := 500000;
    MAX_POWER  :=     20;
    begin
        logical array prime      ( 0 :: MAX_NUMBER );
        integer array powersOf2  ( 1 :: MAX_POWER  );
        integer dpCount;
        % sieve the primes to MAX_NUMBER                                      %
        begin
            integer rootMaxNumber;
            rootMaxNumber:= entier( sqrt( MAX_NUMBER ) );
            prime( 0 ) := prime( 1 ) := false;
            prime( 2 ) := true;
            for i := 3 step 2 until MAX_NUMBER do prime( i ) := true;
            for i := 4 step 2 until MAX_NUMBER do prime( i ) := false;
            for i := 3 step 2 until rootMaxNumber do begin
                if prime( i ) then begin
                    integer i2;
                    i2 := i + i;
                    for s := i * i step i2 until MAX_NUMBER do prime( s ) := false
                end if_prime_i_and_i_lt_rootMaxNumber
            end for_i
        end;
        % table of powers of 2 greater than 2^0 ( up to around 2 000 000 )    %
        %       increase the table size if MAX_NUMBER > 2 000 000             %
        begin
            integer p2;
            p2 := 1;
            for i := 1 until MAX_POWER do begin
                p2             := p2 * 2;
                powersOf2( i ) := p2
            end for_i
        end;
        % the numbers must be odd and not of the form p + 2^n                 %
        % either p is odd and 2^n is even and hence n > 0 and p > 2           %
        %     or 2^n is odd and p is even and hence n = 0 and p = 2           %
        % n = 0, p = 2 - the only possibility is 3                            %
        dpCount := 0;
        for i := 1 step 2 until 3 do begin
            integer p;
            p := 2;
            if p + 1 not = i then begin
                dpCount := dpCount + 1;
                writeon( i_w := 5, i )
            end if_p_plus_1_ne_i
        end for_i ;
        % n > 0, p > 2                                                        %
        for i := 5 step 2 until MAX_NUMBER do begin
            logical found;
            integer p;
            found := false;
            p := 1;
            while p <= MAX_POWER and not found and i > powersOf2( p ) do begin
                found := prime( i - powersOf2( p ) );
                p     := p + 1
            end while_not_found_and_have_a_suitable_power ;
            if not found then begin
                dpCount := dpCount + 1;
                if dpCount <= 50 then begin
                    writeon( i_w := 5, i );
                    if dpCount rem 10 = 0 then write()
                    end
                else if dpCount = 1000 or dpCount = 10000 then begin
                    write( i_w := 5, s_w := 0, "The ", dpCount
                         , "th De Polignac number is ", i
                         )
                end if_various_DePolignac_numbers
            end if_not_found
        end for_i;
        write( i_w := 1, s_w := 0
             , "Found ", dpCount, " De Polignac numbers up to ", MAX_NUMBER
             )
    end
end.
Output:
    1  127  149  251  331  337  373  509  599  701
  757  809  877  905  907  959  977  997 1019 1087
 1199 1207 1211 1243 1259 1271 1477 1529 1541 1549
 1589 1597 1619 1649 1657 1719 1759 1777 1783 1807
 1829 1859 1867 1927 1969 1973 1985 2171 2203 2213
The  1000th De Polignac number is   31941
The 10000th De Polignac number is  273421
Found 19075 De Polignac numbers up to 500000

J

Implementation: 

depolignac=: (1+i.@>.&.-:) (-.,) i.@>.&.(2&^.) +/ i.&.(p:inv)

Task example: 

   5 10$depolignac 3000
   1  127  149  251  331  337  373  509  599  701
 757  809  877  905  907  959  977  997 1019 1087
1199 1207 1211 1243 1259 1271 1477 1529 1541 1549
1589 1597 1619 1649 1657 1719 1759 1777 1783 1807
1829 1859 1867 1927 1969 1973 1985 2171 2203 2213

Stretch: 

   T=: depolignac 3e5
   999 { T
31941
   9999 { T
273421

(We use 999 here for the 1000th number because 0 is the first J index.)

Phix

Translation of: Python
with javascript_semantics
constant MAX_VALUE = 1_000_000,
        all_primes = get_primes_le(MAX_VALUE),
       powers_of_2 = sq_power(2,tagset(floor(log2(MAX_VALUE)),0))
sequence allvalues = repeat(true,MAX_VALUE),
        dePolignac = {}

for i in all_primes do
    for j in powers_of_2 do
        if i + j < MAX_VALUE then
            allvalues[i + j] = false
        end if
    end for
end for
for n=1 to MAX_VALUE by 2 do
    if allvalues[n] then dePolignac &= n end if
end for

printf(1,"First fifty de Polignac numbers:\n%s\n",{join_by(dePolignac[1..50],1,10,"",fmt:="%5d")})
printf(1,"One thousandth: %,d\n",dePolignac[1000])
printf(1,"Ten thousandth: %,d\n",dePolignac[10000])
Output:
First fifty de Polignac numbers:
    1  127  149  251  331  337  373  509  599  701
  757  809  877  905  907  959  977  997 1019 1087
 1199 1207 1211 1243 1259 1271 1477 1529 1541 1549
 1589 1597 1619 1649 1657 1719 1759 1777 1783 1807
 1829 1859 1867 1927 1969 1973 1985 2171 2203 2213

One thousandth: 31,941
Ten thousandth: 273,421

PL/M

Based on the ALGOL 68 sample.

Works with: 8080 PL/M Compiler

... under CP/M (or an emulator) 

100H: /* FIND SOME DE POLIGNAC NUMBERS - POSITIVE ODD NUMBERS THAT CAN'T BE  */
      /* WRITTEN AS P + 2**N FOR SOME PRIME P AND INTEGER N                  */

   DECLARE FALSE LITERALLY '0', TRUE LITERALLY '0FFH';

   /* CP/M SYSTEM CALL AND I/O ROUTINES                                      */
   BDOS:      PROCEDURE( FN, ARG ); DECLARE FN BYTE, ARG ADDRESS; GOTO 5; END;
   PR$CHAR:   PROCEDURE( C ); DECLARE C BYTE;    CALL BDOS( 2, C );  END;
   PR$STRING: PROCEDURE( S ); DECLARE S ADDRESS; CALL BDOS( 9, S );  END;
   PR$NL:     PROCEDURE;   CALL PR$CHAR( 0DH ); CALL PR$CHAR( 0AH ); END;
   PR$NUMBER: PROCEDURE( N ); /* PRINTS A NUMBER IN THE MINIMUN FIELD WIDTH  */
      DECLARE N ADDRESS;
      DECLARE V ADDRESS, N$STR ( 6 )BYTE, W BYTE;
      V = N;
      W = LAST( N$STR );
      N$STR( W ) = '$';
      N$STR( W := W - 1 ) = '0' + ( V MOD 10 );
      DO WHILE( ( V := V / 10 ) > 0 );
         N$STR( W := W - 1 ) = '0' + ( V MOD 10 );
      END;
      CALL PR$STRING( .N$STR( W ) );
   END PR$NUMBER;
   /* END SYSTEM CALL AND I/O ROUTINES                                       */

   DECLARE MAX$DP        LITERALLY '4000', /* MAXIMUM NUMBER TO CONSIDER     */
           MAX$DP$PLUS$1 LITERALLY '4001'; /* MAX$DP + 1 FOR ARRAY BOUNDS    */

   /* SIEVE THE PRIMES TO MAX$DP                                             */
   DECLARE PRIME ( MAX$DP$PLUS$1 )BYTE;
   DO;
      DECLARE ( I, S ) ADDRESS;
      PRIME( 0 ),  PRIME( 1 ) = FALSE;
      PRIME( 2 ) = TRUE;
      DO I = 3 TO LAST( PRIME ) BY 2; PRIME( I ) = TRUE;  END;
      DO I = 4 TO LAST( PRIME ) BY 2; PRIME( I ) = FALSE; END;
      DO I = 3 TO LAST( PRIME ) / 2 BY 2;
         IF PRIME( I ) THEN DO;
            DO S = I + I TO LAST( PRIME ) BY I; PRIME( S ) = FALSE; END;
         END;
      END;
   END;

   /* TABLE OF POWERS OF 2 UP TO MAX$DP                                      */
   DECLARE POWERS$OF$2 ( 12 )ADDRESS
           INITIAL( 1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048 );

   /* DISPLAYS A DE POLIGNAC NUMBER, IF NECESSARY                            */
   SHOW$DE$POLIGNAC: PROCEDURE( DP$NUMBER );
      DECLARE ( DP$NUMBER ) ADDRESS;
      CALL PR$CHAR( ' ' );
      IF DP$NUMBER <   10 THEN CALL PR$CHAR( ' ' );
      IF DP$NUMBER <  100 THEN CALL PR$CHAR( ' ' );
      IF DP$NUMBER < 1000 THEN CALL PR$CHAR( ' ' );
      CALL PR$NUMBER( DP$NUMBER );
   END SHOW$DE$POLIGNAC;

   DECLARE ( I, P, COUNT ) ADDRESS;
   DECLARE FOUND BYTE;

   /* THE NUMBERS MUST BE ODD AND NOT OF THE FORM P + 2**N                */
   /* EITHER P IS ODD AND 2**N IS EVEN AND HENCE N > 0 AND P > 2          */
   /*     OR 2**N IS ODD AND P IS EVEN AND HENCE N = 0 AND P = 2          */

   /* N = 0, P = 2 - THE ONLY POSSIBILITY IS 3                            */
   DO I = 1 TO 3 BY 2;
      P = 2;
      IF P + 1 <> I THEN DO;
         COUNT = COUNT + 1;
         CALL SHOW$DE$POLIGNAC( I );
      END;
   END;
   /* N > 0, P > 2                                                        */
   I = 3;
   DO WHILE I < MAX$DP AND COUNT < 50;
      I = I + 2;
      FOUND = FALSE;
      P = 1;
      DO WHILE NOT FOUND
           AND P <= LAST( POWERS$OF$2 )
           AND I > POWERS$OF$2( P );
         FOUND = PRIME( I - POWERS$OF$2( P ) );
         P     = P + 1;
      END;
      IF NOT FOUND THEN DO;
         CALL SHOW$DE$POLIGNAC( I );
         IF ( COUNT := COUNT + 1 ) MOD 10 = 0 THEN CALL PR$NL;
      END;
   END;

EOF
Output:
    1  127  149  251  331  337  373  509  599  701
  757  809  877  905  907  959  977  997 1019 1087
 1199 1207 1211 1243 1259 1271 1477 1529 1541 1549
 1589 1597 1619 1649 1657 1719 1759 1777 1783 1807
 1829 1859 1867 1927 1969 1973 1985 2171 2203 2213

Python

''' Rosetta code rosettacode.org/wiki/De_Polignac_numbers '''

from sympy import isprime
from math import log
from numpy import ndarray

max_value = 1_000_000

all_primes = [i for i in range(max_value + 1) if isprime(i)]
powers_of_2 = [2**i for i in range(int(log(max_value, 2)))]

allvalues = ndarray(max_value, dtype=bool)
allvalues[:] = True

for i in all_primes:
    for j in powers_of_2:
        if i + j < max_value:
            allvalues[i + j] = False
        
dePolignac = [n for n in range(1, max_value) if n & 1 == 1 and allvalues[n]]

print('First fifty de Polignac numbers:')
for i, n in enumerate(dePolignac[:50]):
    print(f'{n:5}', end='\n' if (i + 1) % 10 == 0 else '')
    
print(f'\nOne thousandth: {dePolignac[999]:,}')
print(f'\nTen thousandth: {dePolignac[9999]:,}')
Output:
First fifty de Polignac numbers:
    1  127  149  251  331  337  373  509  599  701
  757  809  877  905  907  959  977  997 1019 1087
 1199 1207 1211 1243 1259 1271 1477 1529 1541 1549
 1589 1597 1619 1649 1657 1719 1759 1777 1783 1807
 1829 1859 1867 1927 1969 1973 1985 2171 2203 2213

One thousandth: 31,941

Ten thousandth: 273,421

Raku

use List::Divvy;
use Lingua::EN::Numbers;
constant @po2 = (1..∞).map: 2 ** *;
my @dePolignac = lazy 1, |(2..∞).hyper.map(* × 2 + 1).grep: -> $n { all @po2.&upto($n).map: { !is-prime $n - $_ } };

say "First fifty de Polignac numbers:\n" ~ @dePolignac[^50]».&comma».fmt("%5s").batch(10).join: "\n";
say "\nOne thousandth: " ~ comma @dePolignac[999];
say "\nTen thousandth: " ~ comma @dePolignac[9999];
Output:
First fifty de Polignac numbers:
    1   127   149   251   331   337   373   509   599   701
  757   809   877   905   907   959   977   997 1,019 1,087
1,199 1,207 1,211 1,243 1,259 1,271 1,477 1,529 1,541 1,549
1,589 1,597 1,619 1,649 1,657 1,719 1,759 1,777 1,783 1,807
1,829 1,859 1,867 1,927 1,969 1,973 1,985 2,171 2,203 2,213

One thousandth: 31,941

Ten thousandth: 273,421

Wren

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

var pows2 = (0..19).map { |i| 1 << i }.toList
var dp = [1]
var dp1000
var dp10000
var count = 1
var n = 3
while (true) {
    var found = false
    for (pow in pows2) {
        if (pow > n) break
        if (Int.isPrime(n-pow)) {
            found = true
            break
        }
    }
    if (!found) {
        count = count + 1
        if (count <= 50) {
            dp.add(n)
        } else if (count == 1000) {
            dp1000 = n
        } else if (count == 10000) {
            dp10000 = n
            break
        }
    }
    n = n + 2
}
System.print("First 50 De Polignac numbers:")
Fmt.tprint("$,5d", dp, 10)
Fmt.print("\nOne thousandth: $,d", dp1000)
Fmt.print("\nTen thousandth: $,d", dp10000)
Output:
First 50 De Polignac numbers:
    1   127   149   251   331   337   373   509   599   701 
  757   809   877   905   907   959   977   997 1,019 1,087 
1,199 1,207 1,211 1,243 1,259 1,271 1,477 1,529 1,541 1,549 
1,589 1,597 1,619 1,649 1,657 1,719 1,759 1,777 1,783 1,807 
1,829 1,859 1,867 1,927 1,969 1,973 1,985 2,171 2,203 2,213 

One thousandth: 31,941

Ten thousandth: 273,421

XPL0

Slow. Takes 225 seconds on Pi4. 

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;
];

func IsDePolignac(N);   \Return 'true' if N is a de Polignac number
int  N, P, T;
[for P:= N-1 downto 2 do
    if IsPrime(P) then
        [T:= 1;
        repeat  if T+P = N then return false;
                T:= T<<1;
        until   T+P > N;
        ];
return true;
];

int N, C;
[Format(5, 0);
N:= 1;  C:= 0;
loop    [if IsDePolignac(N) then
            [C:= C+1;
            if C<=50 or C=1000 or C=10000 then
                [RlOut(0, float(N));
                if rem(C/10) = 0 then CrLf(0);
                if C = 10000 then quit;
                ];
            ];
        N:= N+2;
        ];
]
Output:
    1  127  149  251  331  337  373  509  599  701
  757  809  877  905  907  959  977  997 1019 1087
 1199 1207 1211 1243 1259 1271 1477 1529 1541 1549
 1589 1597 1619 1649 1657 1719 1759 1777 1783 1807
 1829 1859 1867 1927 1969 1973 1985 2171 2203 2213
31941
273421
Retrieved from "https://rosettacode.org/wiki/De_Polignac_numbers?oldid=330740"