Ultra useful primes

From Rosetta Code
Task
Ultra useful primes
You are encouraged to solve this task according to the task description, using any language you may know.

An ultra-useful prime is a member of the sequence where each a(n) is the smallest positive integer k such that 2(2n) - k is prime.

k must always be an odd number since 2 to any power is always even.


Task
  • Find and show here, on this page, the first 10 elements of the sequence.


Stretch
  • Find and show the next several elements. (The numbers get really big really fast. Only nineteen elements have been identified as of this writing.)


See also


ALGOL 68

Works with: ALGOL 68G version Any - tested with release 2.8.3.win32

Uses Algol 68G's LONG LONG INT which has programmer-specifiable precision. Uses Miller Rabin primality testing.

BEGIN # find members of the sequence a(n) = smallest k such that 2^(2^n) - k is prime #
    PR precision 650 PR # set number of digits for LONG LOMG INT       #
                        # 2^(2^10) has 308 digits but we need more for #
                        # Miller Rabin primality testing               #
    PR read "primes.incl.a68" PR # include the prime related utilities #
    FOR n TO 10 DO
        LONG LONG INT two up 2 up n = LONG LONG INT( 2 ) ^ ( 2 ^ n );
        FOR i BY 2
        WHILE IF is probably prime( two up 2 up n - i ) THEN
                  # found a sequence member #
                  print( ( " ", whole( i, 0 ) ) );
                  FALSE # stop looking #
              ELSE
                  TRUE # haven't found a sequence member yet #
              FI
        DO SKIP OD
    OD
END
Output:
 1 3 5 15 5 59 159 189 569 105

Arturo

ultraUseful: function [n][
    k: 1
    p: (2^2^n) - k
    while ø [
        if prime? p -> return k
        p: p-2
        k: k+2
    ]
]

print [pad "n" 3 "|" pad.right "k" 4]
print repeat "-" 10
loop 1..10 'x ->
    print [(pad to :string x 3) "|" (pad.right to :string ultraUseful x 4)]
Output:
  n | k    
----------
  1 | 1    
  2 | 3    
  3 | 5    
  4 | 15   
  5 | 5    
  6 | 59   
  7 | 159  
  8 | 189  
  9 | 569  
 10 | 105

C

Translation of: Wren
Library: GMP
#include <stdio.h>
#include <gmp.h>

int a(unsigned int n) {
    int k;
    mpz_t p;
    mpz_init_set_ui(p, 1);
    mpz_mul_2exp(p, p, 1 << n);
    mpz_sub_ui(p, p, 1);
    for (k = 1; ; k += 2) {
        if (mpz_probab_prime_p(p, 15) > 0) return k;
        mpz_sub_ui(p, p, 2);
    }
}

int main() {
    unsigned int n;
    printf(" n   k\n");
    printf("----------\n");
    for (n = 1; n < 15; ++n) printf("%2d   %d\n", n, a(n));
    return 0;
}
Output:
 n   k
----------
 1   1
 2   3
 3   5
 4   15
 5   5
 6   59
 7   159
 8   189
 9   569
10   105
11   1557
12   2549
13   2439
14   13797

Craft Basic

for n = 1 to 10

	let k = -1

	do

		let k = k + 2
		wait

	loop prime(2 ^ (2 ^ n) - k) = 0

	print "n = ", n, " k = ", k

next n

Factor

Works with: Factor version 0.99 2021-06-02
USING: io kernel lists lists.lazy math math.primes prettyprint ;

: useful ( -- list )
    1 lfrom
    [ 2^ 2^ 1 lfrom [ - prime? ] with lfilter car ] lmap-lazy ;

10 useful ltake [ pprint bl ] leach nl
Output:
1 3 5 15 5 59 159 189 569 105 

FreeBASIC

Translation of: Ring
#include "isprime.bas"

Dim As Longint n, k, limit = 10
Dim As ulongint num
For n = 1 To limit
    k = -1
    Do
        k += 2
        num = (2 ^ (2 ^ n)) - k
        If isPrime(num) Then
            Print "n = "; n; " k = "; k
            Exit Do
        End If
    Loop
Next

Sleep

Go

package main

import (
    "fmt"
    big "github.com/ncw/gmp"
)

var two = big.NewInt(2)

func a(n uint) int {
    one := big.NewInt(1)
    p := new(big.Int).Lsh(one, 1 << n)
    p.Sub(p, one)
    for k := 1; ; k += 2 {
        if p.ProbablyPrime(15) {
            return k
        }
        p.Sub(p, two)
    }
}

func main() {
    fmt.Println(" n   k")
    fmt.Println("----------")
    for n := uint(1); n < 14; n++ {
        fmt.Printf("%2d   %d\n", n, a(n))
    }
}
Output:
 n   k
----------
 1   1
 2   3
 3   5
 4   15
 5   5
 6   59
 7   159
 8   189
 9   569
10   105
11   1557
12   2549
13   2439

J

Implementation:

uup=: {{
  ref=. 2^2x^y+1
  k=. 1
  while. -. 1 p: ref-k do. k=. k+2 end.
}}"0

I don't have the patience to get this little laptop to compute the first 10 such elements, so here I only show the first five:

   uup i.10
1 3 5 15 5 59 159 189 569 105


Java

import java.math.BigInteger;

public final class UltraUsefulPrimes {

	public static void main(String[] args) {
		for ( int n = 1; n <= 10; n++ ) {
			showUltraUsefulPrime(n);
		}		
	}
	
	private static void showUltraUsefulPrime(int n) {
		BigInteger prime = BigInteger.ONE.shiftLeft(1 << n);
		BigInteger k = BigInteger.ONE; 
		while ( ! prime.subtract(k).isProbablePrime(20) ) {
			k = k.add(BigInteger.TWO);
		}
		
		System.out.print(k + " ");		   
	}

}
Output:
1 3 5 15 5 59 159 189 569 105 

Julia

using Primes

nearpow2pow2prime(n) = findfirst(k -> isprime(2^(big"2"^n) - k), 1:10000)

@time println([nearpow2pow2prime(n) for n in 1:12])
Output:
[1, 3, 5, 15, 5, 59, 159, 189, 569, 105, 1557, 2549]
  3.896011 seconds (266.08 k allocations: 19.988 MiB, 1.87% compilation time)

Mathematica/Wolfram Language

ClearAll[FindUltraUsefulPrimeK]
FindUltraUsefulPrimeK[n_] := Module[{num, tmp},
  num = 2^(2^n);
  Do[
   If[PrimeQ[num - k],
    tmp = k;
    Break[];
    ]
   ,
   {k, 1, \[Infinity], 2}
   ];
  tmp
  ]
res = FindUltraUsefulPrimeK /@ Range[13];
TableForm[res, TableHeadings -> Automatic]
Output:
1	1
2	3
3	5
4	15
5	5
6	59
7	159
8	189
9	569
10	105
11	1557

Nim

Library: Nim-Integers
import std/strformat
import integers

let One = newInteger(1)

echo " n  k"
var count = 1
var n = 1
while count <= 13:
  var k = 1
  var p = One shl (1 shl n) - k
  while not p.isPrime:
    p -= 2
    k += 2
  echo &"{n:2}  {k}"
  inc n
  inc count
Output:
 n  k
 1  1
 2  3
 3  5
 4  15
 5  5
 6  59
 7  159
 8  189
 9  569
10  105
11  1557
12  2549
13  2439

Perl

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

sub useful {
    my @n = @_;
    my @u;
    for my $n (@n) {
        my $p = 2**(2**$n);
        LOOP: for (my $k = 1; $k < $p; $k += 2) {
            is_prime($p-$k) and push @u, $k and last LOOP;
       }
    }
    @u
}

say join ' ', useful 1..13;
Output:
1 3 5 15 5 59 159 189 569 105 1557 2549 2439

Phix

with javascript_semantics
atom t0 = time()
include mpfr.e
mpz p = mpz_init()
 
function a(integer n)
    mpz_ui_pow_ui(p,2,power(2,n))
    mpz_sub_si(p,p,1)
    integer k = 1
    while not mpz_prime(p) do
        k += 2
        mpz_sub_si(p,p,2)
    end while
    return k
end function
 
for i=1 to 10 do
    printf(1,"%d ",a(i))
end for
if machine_bits()=64 then
    ?elapsed(time()-t0)
    for i=11 to 13 do
        printf(1,"%d ",a(i))
    end for
end if
?elapsed(time()-t0)
Output:
1 3 5 15 5 59 159 189 569 105 "0.0s"
1557 2549 2439 "1 minute and 1s"

Raku

The first 10 take less than a quarter second. 11 through 13, a little under 30 seconds. Drops off a cliff after that.

sub useful ($n) {
    (|$n).map: {
        my $p = 1 +< ( 1 +< $_ );
        ^$p .first: ($p - *).is-prime
    }
}

put useful 1..10;

put useful 11..13;
Output:
1 3 5 15 5 59 159 189 569 105
1557 2549 2439

Python

# useful.py by xing216
from gmpy2 import is_prime
def useful(n):
    k = 1
    is_useful = False
    while is_useful == False:
        if is_prime(2**(2**n) - k):
            is_useful = True
            break
        k += 2
    return k
if __name__ == "__main__":
    print("n | k")
    for i in range(1,14):
        print(f"{i:<4}{useful(i)}")
Output:
n | K
1   1
2   3
3   5
4   15
5   5
6   59
7   159
8   189
9   569
10  105
11  1557
12  2549
13  2439

Ring

see "works..." + nl
limit = 10

for n = 1 to limit
    k = -1
    while true
          k = k + 2
          num = pow(2,pow(2,n)) - k
          if isPrime(num)
             ? "n = " + n + " k = " + k
             exit
          ok
    end
next
see "done.." + nl

func isPrime num
     if (num <= 1) return 0 ok
     if (num % 2 = 0 and num != 2) return 0 ok
     for i = 3 to floor(num / 2) -1 step 2
         if (num % i = 0) return 0 ok
     next
     return 1
Output:
works...
n = 1 k = 1
n = 2 k = 3
n = 3 k = 5
n = 4 k = 15
n = 5 k = 5
n = 6 k = 59
n = 7 k = 159
n = 8 k = 189
n = 9 k = 569
n = 10 k = 105
done...

Ruby

The 'prime?' method in Ruby's OpenSSL (standard) lib is much faster than the 'prime' lib. 0.05 sec. for this, about a minute for the next three.

require 'openssl'

(1..10).each do |n|
   pow = 2 ** (2 ** n)
   print "#{n}:\t"
   puts (1..).step(2).detect{|k| OpenSSL::BN.new(pow-k).prime?}
end
Output:
1:      1
2:      3
3:      5
4:      15
5:      5
6:      59
7:      159
8:      189
9:      569
10:     105

Sidef

say(" n   k")
say("----------")

for n in (1..13) {
    var t = 2**(2**n)
    printf("%2d   %d\n", n, {|k| t - k -> is_prob_prime }.first)
}
Output:
 n   k
----------
 1   1
 2   3
 3   5
 4   15
 5   5
 6   59
 7   159
 8   189
 9   569
10   105
11   1557
12   2549
13   2439

(takes ~20 seconds)

V (Vlang)

import math

fn main() {
	limit := 10 // depending on computer, higher numbers = longer times
	mut num, mut k := i64(0), i64(0)
	println("n   k\n-------")
	for n in 1..limit {
		k = -1
		for n < limit {
			  k = k + 2
			  num = math.powi(2, math.powi(2 , n)) - k
			  if is_prime(num) == true {
				 println("${n}   ${k}")
				 break
			  }
		}
	}
}

fn is_prime(num i64) bool {
     if num <= 1 {return false}
     if num % 2 == 0 && num != 2 {return false}
	 for idx := 3; idx <= math.floor(num / 2) - 1; idx += 2 {
         if num % idx == 0 {return false}
     }
     return true
}
Output:
n   k
-------
1   1
2   3
3   5
4   15
5   5
6   59
7   159
8   189
9   569
10  105

Wren

CLI

Library: Wren-big
Library: Wren-fmt

Manages to find the first ten but takes 84 seconds to do so.

import "./big" for BigInt
import "./fmt" for Fmt

var one = BigInt.one
var two = BigInt.two

var a = Fn.new { |n|
    var p = (BigInt.one << (1 << n)) - one
    var k = 1
    while (true) {
        if (p.isProbablePrime(5)) return k
        p = p - two
        k = k + 2
    }
}

System.print(" n   k")
System.print("----------")
for (n in 1..10) Fmt.print("$2d   $d", n, a.call(n))
Output:
 n   k
----------
 1   1
 2   3
 3   5
 4   15
 5   5
 6   59
 7   159
 8   189
 9   569
10   105

Embedded

Library: Wren-gmp

The following takes about 17 seconds to get to n = 13 but 7 minutes 10 seconds to reach n = 14. I didn't bother after that.

import "./gmp" for Mpz
import "./fmt" for Fmt

var one = Mpz.one
var two = Mpz.two

var a = Fn.new { |n|
    var p = Mpz.one.lsh(1 << n).sub(one)
    var k = 1
    while (true) {
        if (p.probPrime(15) > 0) return k
        p.sub(two)
        k = k + 2
    }
}

System.print(" n   k")
System.print("----------")
for (n in 1..14) Fmt.print("$2d   $d", n, a.call(n))
Output:
 n   k
----------
 1   1
 2   3
 3   5
 4   15
 5   5
 6   59
 7   159
 8   189
 9   569
10   105
11   1557
12   2549
13   2439
14   13797