# Wagstaff primes

Wagstaff primes
You are encouraged to solve this task according to the task description, using any language you may know.
Definition

A Wagstaff prime is a prime number of the form (2^p + 1)/3 where the exponent p is an odd prime.

Example

(2^5 + 1)/3 = 11 is a Wagstaff prime because both 5 and 11 are primes.

Task

Find and show here the first 10 Wagstaff primes and their corresponding exponents p.

Stretch (requires arbitrary precision integers)

Find and show here the exponents p corresponding to the next 14 Wagstaff primes (not the primes themselves) and any more that you have the patience for.

When testing for primality, you may use a method which determines that a large number is probably prime with reasonable certainty.

Note

It can be shown (see talk page) that (2^p + 1)/3 is always integral if p is odd. So there's no need to check for that prior to checking for primality.

References

## ALGOL 68

Basic task - assumes LONG INT is at least 64 bit.

```BEGIN # find some Wagstaff primes: primes of the form ( 2^p + 1 ) / 3 #
#                            where p is an odd prime            #
INT max wagstaff = 10; # number of Wagstaff primes to find        #
INT w count     :=  0; # numbdr of Wagstaff primes found so far   #
# sieve the primes up to 200, hopefully enough...                 #
[ 0 : 200 ]BOOL primes;
primes[ 0 ] := primes[ 1 ] := FALSE;
primes[ 2 ] := TRUE;
FOR i FROM 3 BY 2 TO UPB primes DO primes[ i ] := TRUE  OD;
FOR i FROM 4 BY 2 TO UPB primes DO primes[ i ] := FALSE OD;
FOR i FROM 3 BY 2 TO ENTIER sqrt( UPB primes ) DO
IF primes[ i ] THEN
FOR s FROM i * i BY i + i TO UPB primes DO primes[ s ] := FALSE OD
FI
OD;
# attempt to find the Wagstaff primes                              #
LONG INT power of 2 := 2; # 2^1                                    #
FOR p FROM 3 BY 2 WHILE w count < max wagstaff DO
power of 2 *:= 4;
IF primes[ p ] THEN
LONG INT w := ( power of 2 + 1 ) OVER 3;
# check w is prime - trial division                        #
BOOL is prime := TRUE;
LONG INT n    := 3;
WHILE n * n <= w AND is prime DO
is prime := w MOD n /= 0;
n       +:= 2
OD;
IF is prime THEN
# have another Wagstaff prime                          #
w count +:= 1;
print( ( whole( w count, -2 )
, ": "
, whole( p, -4 )
, ": "
, whole( w, 0 )
, newline
)
)
FI
FI
OD
END```
Output:
``` 1:    3: 3
2:    5: 11
3:    7: 43
4:   11: 683
5:   13: 2731
6:   17: 43691
7:   19: 174763
8:   23: 2796203
9:   31: 715827883
10:   43: 2932031007403
```

## ALGOL W

As Algol W integers are limited to 32 bits, 64 bit reals are used which can accurately represent integers up to 2^53.

```begin % find some Wagstaff primes: primes of the form ( 2^p + 1 ) / 3        %
%                            where p is an odd prime                   %

integer   wCount;

% returns true if d exactly divides v, false otherwise                   %
logical procedure divides( long real value d, v ) ;
begin
long real q, p10;
q := v / d;
p10 := 1;
while p10 * 10 < q do begin
p10 := p10 * 10
end while_p10_lt_q ;
while p10 >= 1 do begin
while q >= p10 do q := q - p10;
p10 := p10 / 10
end while_p10_ge_1 ;
q = 0
end divides ;

% prints v as a 14 digit integer number                                  %
procedure printI14( long real value v ) ;
begin
integer   iv;
long real r;
r := abs( v );
if v < 0 then writeon( s_w := 0, "-" );
iv := truncate( r / 1000000 );
if iv < 1 then begin
writeon( i_w := 8, s_w := 0, "      ", truncate( r ) )
end
else begin
string(6) sValue;
writeon( i_w := 8, s_w := 0, iv );
iv := truncate( abs( r ) - ( iv * 1000000.0 ) );
for sPos := 5 step -1 until 0 do begin
sValue( sPos // 1 ) := code( ( iv rem 10 ) + decode( "0" ) );
iv := iv div 10
end for_sPos;
writeon( s_w := 0, sValue )
end if_iv_lt_1__
end printI ;

wCount     := 0;
% find the Wagstaff primes using long reals to hold the numbers, which   %
% accurately represent integers up to 2^53, so only consider primes < 53 %
for p := 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47 do begin
long real w, powerOfTwo;
logical   isPrime;
powerOfTwo := 1;
for i := 1 until p do powerOfTwo := powerOfTwo * 2;
w := ( powerOfTwo + 1 ) / 3;
isPrime := not divides( 2, w );
if isPrime then begin
integer   f, toNext;
long real f2;
f      := 3;
f2     := 9;
toNext := 16;
while isPrime and f2 <= w do begin
isPrime := not divides( f, w );
f       := f  + 2;
f2      := f2 + toNext;
toNext  := toNext + 8
end while_isPrime_and_f2_le_x
end if_isPrime ;
if isPrime then begin
wCount := wCount + 1;
write( i_w := 3, s_w := 0, wCount, ": ", p, "    " );
if p >= 32 then printI14( w )
else begin
integer iw;
iw := truncate( w );
writeon( i_w := 14, s_w := 0, iw )
end of_p_ge_32__ ;
if wCount >= 10 then goto done % stop at 10 Wagstaff primes       %
end if_isPrime
end for_p ;
done:
end.```
Output:
```  1:   3                 3
2:   5                11
3:   7                43
4:  11               683
5:  13              2731
6:  17             43691
7:  19            174763
8:  23           2796203
9:  31         715827883
10:  43     2932031007403
```

## Arturo

```wagstaff?: function [e][
and? -> prime? e -> prime? (1+2^e)/3
]

summarize: function [n][
n: ~"|n|"
s: size n
if s > 20 -> n: ((take n 10)++"...")++drop.times:s-10 n
n ++ ~" (|s| digits)"
]

exponents: select.first:24 range.step:2 1 ∞ => wagstaff?
loop.with:'i exponents 'x -> print [
pad ~"|i+1|:" 3 pad ~"|x| -" 6 summarize (1+2^x)/3
]
```
Output:
```

p = 3, m = 3
p = 5, m = 11
p = 7, m = 43
p = 11, m = 683
p = 13, m = 2731
p = 17, m = 43691
p = 19, m = 174763
p = 23, m = 2796203
p = 31, m = 715827883
p = 43, m = 2932031007403
p = 61, m = 768614336404564651
p = 79, m = 2014876366...8195784363 (24 digits)
p = 101, m = 8451004001...1135470251 (30 digits)
p = 127, m = 5671372782...8628035243 (38 digits)
p = 167, m = 6235740319...3121833643 (50 digits)
p = 191, m = 1046183622...7339085483 (58 digits)
p = 199, m = 2678230073...8805883563 (60 digits)
p = 313, m = 5562466239...8130434731 (94 digits)
p = 347, m = 9556244233...1903606443 (104 digits)
p = 701, m = 3506757267...7823854251 (211 digits)
p = 1709, m = 9619252724...9070528171 (514 digits)
p = 2617, m = 2081504709...3435947691 (788 digits)
p = 3539, m = 7379609820...6486497963 (1065 digits)
p = 5807, m = 4018496237...6663568043 (1748 digits)

```

## BASIC256

Translation of: FreeBASIC
```function isPrime(v)
if v < 2 then return False
if v mod 2 = 0 then return v = 2
if v mod 3 = 0 then return v = 3
d = 5
while d * d <= v
if v mod d = 0 then return False else d += 2
end while
return True
end function

subroutine Wagstaff(num)
pri = 1
wcount = 0
wag = 0
while wcount < num
pri += 2
if isPrime(pri) then
wag = (2 ^ pri + 1) / 3
if isPrime(wag) then
wcount += 1
print rjust(wcount,2); ": "; rjust(pri,2); " => "; int(wag)
end if
end if
end while
end subroutine

call Wagstaff(9) #BASIC-256 does not allow larger numbers
end```

## C

Library: GMP

This takes about 24 seconds to find the first 24 Wagstaff primes but 589 seconds to find the first 29.

```#include <stdio.h>
#include <string.h>
#include <gmp.h>

int main() {
const int limit = 29;
int count = 0;
char tmp[40];
mpz_t p, w;
mpz_init_set_ui(p, 1);
mpz_init(w);
while (count < limit) {
mpz_nextprime(p, p);
mpz_set_ui(w, 1);
unsigned long ulp = mpz_get_ui(p);
mpz_mul_2exp(w, w, ulp);
mpz_add_ui(w, w, 1);
mpz_tdiv_q_ui(w, w, 3);
if (mpz_probab_prime_p(w, 15) > 0) {
++count;
char *ws = mpz_get_str(NULL, 10, w);
size_t le = strlen(ws);
if (le < 34) {
strcpy(tmp, ws);
} else {
strncpy(tmp, ws, 15);
strcpy(tmp + 15, "...");
strncpy(tmp + 18, ws + le - 15, 16);
}
printf("%5lu: %s", ulp, tmp);
if (le >=34) printf( " (%ld digits)", le);
printf("\n");
}
}
return 0;
}
```
Output:
```    3: 3
5: 11
7: 43
11: 683
13: 2731
17: 43691
19: 174763
23: 2796203
31: 715827883
43: 2932031007403
61: 768614336404564651
79: 201487636602438195784363
101: 845100400152152934331135470251
127: 567137278201564...101238628035243 (38 digits)
167: 623574031927851...838653121833643 (50 digits)
191: 104618362256444...574077339085483 (58 digits)
199: 267823007376498...963798805883563 (60 digits)
313: 556246623937737...099018130434731 (94 digits)
347: 955624423329196...712921903606443 (104 digits)
701: 350675726769891...862167823854251 (211 digits)
1709: 961925272487010...857299070528171 (514 digits)
2617: 208150470990263...847933435947691 (788 digits)
3539: 737960982013072...304686486497963 (1065 digits)
5807: 401849623734300...466686663568043 (1748 digits)
10501: 435374724597165...340431558126251 (3161 digits)
10691: 683222859828090...209259644365483 (3218 digits)
11279: 692149388152358...081084271176363 (3395 digits)
12391: 385083595083686...967789680077483 (3730 digits)
14479: 136831460986221...678801865288363 (4359 digits)
```

## C++

Library: GMP
Library: Primesieve
```#include <gmpxx.h>
#include <primesieve.hpp>

#include <iostream>

using big_int = mpz_class;

std::string to_string(const big_int& num, size_t max_digits) {
std::string str = num.get_str();
size_t len = str.size();
if (len > max_digits) {
str.replace(max_digits / 2, len - max_digits, "...");
str += " (";
str += std::to_string(len);
str += " digits)";
}
return str;
}

bool is_probably_prime(const big_int& n) {
return mpz_probab_prime_p(n.get_mpz_t(), 25) != 0;
}

int main() {
const big_int one(1);
primesieve::iterator pi;
pi.next_prime();
for (int i = 0; i < 24;) {
uint64_t p = pi.next_prime();
big_int n = ((one << p) + 1) / 3;
if (is_probably_prime(n))
std::cout << ++i << ": " << p << " - " << to_string(n, 30) << '\n';
}
}
```
Output:
```1: 3 - 3
2: 5 - 11
3: 7 - 43
4: 11 - 683
5: 13 - 2731
6: 17 - 43691
7: 19 - 174763
8: 23 - 2796203
9: 31 - 715827883
10: 43 - 2932031007403
11: 61 - 768614336404564651
12: 79 - 201487636602438195784363
13: 101 - 845100400152152934331135470251
14: 127 - 567137278201564...101238628035243 (38 digits)
15: 167 - 623574031927851...838653121833643 (50 digits)
16: 191 - 104618362256444...574077339085483 (58 digits)
17: 199 - 267823007376498...963798805883563 (60 digits)
18: 313 - 556246623937737...099018130434731 (94 digits)
19: 347 - 955624423329196...712921903606443 (104 digits)
20: 701 - 350675726769891...862167823854251 (211 digits)
21: 1709 - 961925272487010...857299070528171 (514 digits)
22: 2617 - 208150470990263...847933435947691 (788 digits)
23: 3539 - 737960982013072...304686486497963 (1065 digits)
24: 5807 - 401849623734300...466686663568043 (1748 digits)
```

## Craft Basic

```'Craft Basic may only calculate up to 9

let n = 9
let p = 1

do
let p = p + 2

if prime(p) then

let w = (2 ^ p + 1) / 3

if prime(w) then

let c = c + 1
print tab, c, tab, p, tab, w

endif

endif

title p

wait

loop c < n
```
Output:
```	1	3	3
2	5	11
3	7	43
4	11	683
5	13	2731
6	17	43691
7	19	174763
8	23	2796203
9	31	715827883
```

## Delphi

Works with: Delphi version 6.0

Finds Wagstaff primes up to 2^64, the biggest integer 32-bit Delphi supports.

```procedure WagstaffPrimes(Memo: TMemo);
{Finds Wagstaff primes up to 6^64}
var P,B,R: int64;
var Cnt: integer;
begin
{B is int64 to force 64-bit arithmetic}
P:=0; B:=1; Cnt:=0;
Memo.Lines.Add(' #:  P              (2^P + 1)/3');
while P<64 do
begin
R:=((B shl P) + 1) div 3;
if IsPrime(P) and IsPrime(R) then
begin
Inc(Cnt);
Memo.Lines.Add(Format('%2d: %2d %24.0n',[Cnt,P,R+0.0]));
end;
Inc(P);
end;
end;
```
Output:
``` #:  P              (2^P + 1)/3
1:  3                        3
2:  5                       11
3:  7                       43
4: 11                      683
5: 13                    2,731
6: 17                   43,691
7: 19                  174,763
8: 23                2,796,203
9: 31              715,827,883
10: 43        2,932,031,007,403
11: 61  768,614,336,404,564,651

Elapsed Time: 29.504 Sec.

```

## EasyLang

```func prime n .
if n mod 2 = 0 and n > 2
return 0
.
i = 3
while i <= sqrt n
if n mod i = 0
return 0
.
i += 2
.
return 1
.
pri = 1
while nwag <> 10
pri += 2
if prime pri = 1
wag = (pow 2 pri + 1) / 3
if prime wag = 1
nwag += 1
print pri & " => " & wag
.
.
.```

## F#

This task uses Extensible Prime Generator (F#)

```// Wagstaff primes. Nigel Galloway: September 15th., 2022
let isWagstaff n=let mutable g=(1I+2I**n)/3I in if Open.Numeric.Primes.MillerRabin.IsProbablePrime &g then Some (n,g) else None
primes32()|>Seq.choose isWagstaff|>Seq.take 10|>Seq.iter(fun (n,g)->printfn \$"%d{n}->%A{g}")
primes32()|>Seq.choose isWagstaff|>Seq.skip 10|>Seq,take 14|>Seq.iter(fun(n,_)->printf \$"%d{n} "); printfn ""
```
Output:
```3->3
5->11
7->43
11->683
13->2731
17->43691
19->174763
23->2796203
31->715827883
43->2932031007403

61 79 101 127 167 191 199 313 347 701 1709 2617 3539 5807
```

## FreeBASIC

Translation of: Python
```Function isPrime(Byval num As Ulongint) As Boolean
If num < 2 Then Return False
If num Mod 2 = 0 Then Return num = 2
If num Mod 3 = 0 Then Return num = 3
Dim d As Uinteger = 5
While d * d <= num
If num Mod d = 0 Then Return False Else d += 2
Wend
Return True
End Function

Sub Wagstaff(num As Ulongint)
Dim As Ulongint pri = 1, wcount = 0, wag
While wcount < num
pri += 2
If isPrime(pri) Then
wag = (2 ^ pri + 1) / 3
If isPrime(wag) Then
wcount += 1
Print Using "###: ### => ##,###,###,###,###"; wcount; pri; wag
End If
End If
Wend
End Sub

Wagstaff(10)
Sleep```
Output:
```  1:   3 =>                  3
2:   5 =>                 11
3:   7 =>                 43
4:  11 =>                683
5:  13 =>              2,731
6:  17 =>             43,691
7:  19 =>            174,763
8:  23 =>          2,796,203
9:  31 =>        715,827,883
10:  43 =>  2,932,031,007,403```

## Gambas

Translation of: FreeBASIC
```Use "isprime.bas"

Public Sub Main()

Wagstaff(10)

End

Sub Wagstaff(num As Long)

Dim pri As Long = 1
Dim wcount As Long = 0
Dim wag As Long

While wcount < num
pri = pri + 2
If isPrime(pri) Then
wag = (2 ^ pri + 1) / 3
If isPrime(wag) Then
wcount += 1
Print Format\$(Str(wcount), "###"); ": "; Format\$(Str(pri), "###"); " => "; Int(wag)
End If
End If
Wend

End Sub
```
Output:
`Similar to FreeBASIC entry.`

## Go

Translation of: C
Library: Go-rcu

Timings: about 25 seconds for first 24 and 596 seconds for first 29.

```package main

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

func main() {
const limit = 29
count := 0
p := 1
one := big.NewInt(1)
three := big.NewInt(3)
w := new(big.Int)
for count < limit {
for {
p += 2
if rcu.IsPrime(p) {
break
}
}
w.SetUint64(1)
w.Lsh(w, uint(p))
w.Add(w, one)
w.Quo(w, three)
if w.ProbablyPrime(15) {
count++
ws := w.String()
le := len(ws)
if le >= 34 {
ws = ws[0:15] + "..." + ws[le-15:]
}
fmt.Printf("%5d: %s", p, ws)
if le >= 34 {
fmt.Printf(" (%d digits)", le)
}
println()
}
}
}
```
Output:
```Same as C entry
```

## J

```   (,.f)p:I.1 p:(f=.3%~1+2x^])p:i.60
3                                                            3
5                                                           11
7                                                           43
11                                                          683
13                                                         2731
17                                                        43691
19                                                       174763
23                                                      2796203
31                                                    715827883
43                                                2932031007403
61                                           768614336404564651
79                                     201487636602438195784363
101                               845100400152152934331135470251
127                       56713727820156410577229101238628035243
167           62357403192785191176690552862561408838653121833643
191   1046183622564446793972631570534611069350392574077339085483
199 267823007376498379256993682056860433753700498963798805883563
```

Stretch goal (times are elapsed times, machine was a 1.2GHz intel i3):

```{{
T0=. 6!:1''
f=. 3%~1+2x^]
c=. 0
echo 'count power seconds'
for_p. p:i.1e4 do.
if. 1 p: f p do.
c=. c+1
echo 5 6 8j3":c, p, (6!:1'')-T0
if. 24 <: c do. return. end.
end.
end.
}}_
count power seconds
1     3   0.004
2     5   0.006
3     7   0.008
4    11   0.013
5    13   0.018
6    17   0.021
7    19   0.025
8    23   0.028
9    31   0.030
10    43   0.033
11    61   0.035
12    79   0.039
13   101   0.044
14   127   0.052
15   167   0.062
16   191   0.075
17   199   0.089
18   313   0.120
19   347   0.167
20   701   0.436
21  1709   4.140
22  2617  16.035
23  3539  45.089
24  5807 181.280
```

## Julia

• Requires Primes.jl
```using Primes

function wagstaffpair(p::Integer)
isodd(p) || return (false, nothing)
isprime(p) || return (false, nothing)

m = (2^big(p) + 1) ÷ 3

isprime(m) || return (false, nothing)

return (true, m)
end

function findn_wagstaff_pairs(n_to_find::T) where T <: Integer
pairs = Tuple{T, BigInt}[]
count = 0
i = 2
while count < n_to_find
iswag, m = wagstaffpair(i)
iswag && push!(pairs, (i, m))
count += iswag
i += 1
end
return pairs
end

function println_wagstaff(pair; max_digit_display::Integer=20)
p, m = pair
mstr = string(m)

if length(mstr) > max_digit_display
umiddle = cld(max_digit_display, 2)
lmiddle = fld(max_digit_display, 2)
mstr = join((mstr[1:umiddle], "...", mstr[end-lmiddle+1:end],
" (\$(length(mstr)) digits)"))
end

println("p = \$p, m = \$mstr")
end

foreach(println_wagstaff, findn_wagstaff_pairs(24))
```
```p = 3, m = 3
p = 5, m = 11
p = 7, m = 43
p = 11, m = 683
p = 13, m = 2731
p = 17, m = 43691
p = 19, m = 174763
p = 23, m = 2796203
p = 31, m = 715827883
p = 43, m = 2932031007403
p = 61, m = 768614336404564651
p = 79, m = 2014876366...8195784363 (24 digits)
p = 101, m = 8451004001...1135470251 (30 digits)
p = 127, m = 5671372782...8628035243 (38 digits)
p = 167, m = 6235740319...3121833643 (50 digits)
p = 191, m = 1046183622...7339085483 (58 digits)
p = 199, m = 2678230073...8805883563 (60 digits)
p = 313, m = 5562466239...8130434731 (94 digits)
p = 347, m = 9556244233...1903606443 (104 digits)
p = 701, m = 3506757267...7823854251 (211 digits)
p = 1709, m = 9619252724...9070528171 (514 digits)
p = 2617, m = 2081504709...3435947691 (788 digits)
p = 3539, m = 7379609820...6486497963 (1065 digits)
p = 5807, m = 4018496237...6663568043 (1748 digits)

```

## Java

```import java.math.BigInteger;

public class Main {
public static void main(String[] args) {
BigInteger d = new BigInteger("3"), a;
int lmt = 25, sl, c = 0;
for (int i = 3; i < 5808; ) {
a = BigInteger.ONE.shiftLeft(i).add(BigInteger.ONE).divide(d);
if (a.isProbablePrime(1)) {
System.out.printf("%2d %4d ", ++c, i);
String s = a.toString(); sl = s.length();
if (sl < lmt) System.out.println(a);
else System.out.println(s.substring(0, 11) + ".." + s.substring(sl - 11, sl) + " " + sl + " digits");
}
i = BigInteger.valueOf(i).nextProbablePrime().intValue();
}
}
}
```
Output:

Takes under 30 seconds @ Tio.run

``` 1    3 3
2    5 11
3    7 43
4   11 683
5   13 2731
6   17 43691
7   19 174763
8   23 2796203
9   31 715827883
10   43 2932031007403
11   61 768614336404564651
12   79 201487636602438195784363
13  101 84510040015..31135470251 30 digits
14  127 56713727820..38628035243 38 digits
15  167 62357403192..53121833643 50 digits
16  191 10461836225..77339085483 58 digits
17  199 26782300737..98805883563 60 digits
18  313 55624662393..18130434731 94 digits
19  347 95562442332..21903606443 104 digits
20  701 35067572676..67823854251 211 digits
21 1709 96192527248..99070528171 514 digits
22 2617 20815047099..33435947691 788 digits
23 3539 73796098201..86486497963 1065 digits
24 5807 40184962373..86663568043 1748 digits
```

## Nim

Library: Nim-Integers
```import std/strformat
import integers

func compress(str: string; size: int): string =
if str.len <= 2 * size: str
else: &"{str[0..<size]}...{str[^size..^1]} ({str.len} digits)"

echo "First 24 Wagstaff primes:"
let One = newInteger(1)
var count = 0
var p = 3
while count < 24:
if p.isPrime:
let n = (One shl p + 1) div 3
if n.isPrime:
inc count
echo &"{p:4}: {compress(\$n, 15)}"
inc p, 2
```
Output:
```First 24 Wagstaff primes:
3: 3
5: 11
7: 43
11: 683
13: 2731
17: 43691
19: 174763
23: 2796203
31: 715827883
43: 2932031007403
61: 768614336404564651
79: 201487636602438195784363
101: 845100400152152934331135470251
127: 567137278201564...101238628035243 (38 digits)
167: 623574031927851...838653121833643 (50 digits)
191: 104618362256444...574077339085483 (58 digits)
199: 267823007376498...963798805883563 (60 digits)
313: 556246623937737...099018130434731 (94 digits)
347: 955624423329196...712921903606443 (104 digits)
701: 350675726769891...862167823854251 (211 digits)
1709: 961925272487010...857299070528171 (514 digits)
2617: 208150470990263...847933435947691 (788 digits)
3539: 737960982013072...304686486497963 (1065 digits)
5807: 401849623734300...466686663568043 (1748 digits)
```

## OCaml

```let is_prime n =
let rec test x =
let q = n / x in x > q || x * q <> n && n mod (x + 2) <> 0 && test (x + 6)
in if n < 5 then n lor 1 = 3 else n land 1 <> 0 && n mod 3 <> 0 && test 5

let is_wagstaff n =
let w = succ (1 lsl n) / 3 in
if is_prime n && is_prime w then Some (n, w) else None

let () =
let show (p, w) = Printf.printf "%u -> %u%!\n" p w in
Seq.(ints 3 |> filter_map is_wagstaff |> take 11 |> iter show)
```
Output:
```3 -> 3
5 -> 11
7 -> 43
11 -> 683
13 -> 2731
17 -> 43691
19 -> 174763
23 -> 2796203
31 -> 715827883
43 -> 2932031007403
61 -> 768614336404564651
```

## PARI/GP

wprp.gp

```\\ compiler: gp2c option: gp2c-run -g wprp.gp

/* wprp(p): input odd prime p > 5 . */
/* returns 1 if (2^p+1)/3 is a Wagstaff probable prime. */
wprp(p) = {

/* trial division up to a reasonable depth (time ratio tdiv/llr ≈ 0.2) */
my(l=log(p), ld=log(l));
forprimestep(q = 1, sqr(ld)^(l/log(2))\4, p+p,
if(Mod(2,q)^p == -1, return)
);

/*            From  R. & H. LIFCHITZ July 2000           */
/*   see: http://primenumbers.net/Documents/TestNP.pdf   */
/* if (2^p+1)/3 is prime ==> 25^2^(p-1) ≡ 25 (mod 2^p+1) */
/* Lucas-Lehmer-Riesel test with fast modular reduction. */
my(s=25, m=2^p-1);
for(i = 2, p,
s = sqr(s);
s = bitand(s,m) - s>>p
);
s==25
};      /* end wprp */

/* get exponents of Wagstaff prps in range [a,b] */
wprprun(a, b) = {
my(t=0, c=0, thr=default(nbthreads));
a = max(a,3);
gettime();
if(a <= 5,
if(a == 3, c++; p = 3; printf("#%d\tW%d\t%2dmin, %2d,%03d ms\n", c, p, t\60000%60, t\1000%60, t%1000));
c++; p = 5; printf("#%d\tW%d\t%2dmin, %2d,%03d ms\n", c, p, t\60000%60, t\1000%60, t%1000);
a = 7
);
parforprime(p= a, b, wprp(p), d,       /* wprp(p) -> d  copy from parallel world into real world. */
if(d,
t += gettime()\thr;
c++;
printf("#%d\tW%d\t%2dmin, %2d,%03d ms\n", c, p, t\60000%60, t\1000%60, t%1000)
)
)
};   /* end wprprun */

/* if running wprp as script */
\\ export(wprp);

wprprun(2, 42737)
```
Output:

Done on Intel(R) Core(TM) i5-8250U CPU @ 1.60GHz, 4 hyperthreaded cores.

```#1  	W3       0min,  0,000 ms
#2  	W5       0min,  0,000 ms
#3  	W7       0min,  0,000 ms
#4  	W11      0min,  0,000 ms
#5  	W13      0min,  0,000 ms
#6  	W17      0min,  0,000 ms
#7  	W19      0min,  0,000 ms
#8  	W23      0min,  0,000 ms
#9   	W31      0min,  0,000 ms
#10  	W43      0min,  0,000 ms
#11 	W61      0min,  0,000 ms
#12 	W79      0min,  0,000 ms
#13 	W101     0min,  0,000 ms
#14 	W127	 0min,  0,000 ms
#15 	W167	 0min,  0,000 ms
#16 	W191	 0min,  0,000 ms
#17 	W199	 0min,  0,000 ms
#18 	W313	 0min,  0,000 ms
#19 	W347	 0min,  0,000 ms
#20 	W701	 0min,  0,002 ms
#21 	W1709	 0min,  0,018 ms
#22 	W2617	 0min,  0,055 ms
#23 	W3539	 0min,  0,113 ms
#24 	W5807    0min,  0,422 ms
#25 	W10501   0min,  2,545 ms
#26 	W10691   0min,  2,612 ms
#27 	W11279   0min,  3,113 ms
#28 	W12391   0min,  4,252 ms
#29 	W14479   0min,  7,719 ms
#30 	W42737   7min, 33,111 ms
? ##
***   last result: cpu time 1h, 8,965 ms, real time 7min, 33,449 ms.```

## Perl

First 20 terms are very fast, and 30 terms in just the time it takes to make a cup of coffee! (caveat: IFF you allow for the travel time of a flight to Hawaii to pick up some nice Kona beans first). Seriously, don't try this at home; bigint is mandatory here, which makes it glacially slow...

Library: ntheory
```use v5.36;
use bigint;
use ntheory 'is_prime';

sub abbr (\$d) { my \$l = length \$d; \$l < 61 ? \$d : substr(\$d,0,30) . '..' . substr(\$d,-30) . " (\$l digits)" }

my(\$p,@W) = 2;
until (@W == 30) {
next unless 0 != ++\$p % 2;
push @W, \$p if is_prime(\$p) and is_prime((2**\$p + 1)/3)
}

printf "%2d: %5d - %s\n", \$_+1, \$W[\$_], abbr( (2**\$W[\$_] + 1) / 3) for 0..\$#W;
```
Output:
``` 1:     3 - 3
2:     5 - 11
3:     7 - 43
4:    11 - 683
5:    13 - 2731
6:    17 - 43691
7:    19 - 174763
8:    23 - 2796203
9:    31 - 715827883
10:    43 - 2932031007403
11:    61 - 768614336404564651
12:    79 - 201487636602438195784363
13:   101 - 845100400152152934331135470251
14:   127 - 56713727820156410577229101238628035243
15:   167 - 62357403192785191176690552862561408838653121833643
16:   191 - 1046183622564446793972631570534611069350392574077339085483
17:   199 - 267823007376498379256993682056860433753700498963798805883563
18:   313 - 556246623937737000623703569314..370363869220418099018130434731 (94 digits)
19:   347 - 955624423329196463171175373042..686867002917439712921903606443 (104 digits)
20:   701 - 350675726769891567149399325525..838589548478180862167823854251 (211 digits)
21:  1709 - 961925272487010690059185752191..922686975212036857299070528171 (514 digits)
22:  2617 - 208150470990263232578066751439..039652079499645847933435947691 (788 digits)
23:  3539 - 737960982013072251717827114247..199733732972086304686486497963 (1065 digits)
24:  5807 - 401849623734300407681642626119..205019459748642466686663568043 (1748 digits)
25: 10501 - 435374724597165148564719869736..420302813453630340431558126251 (3161 digits)
26: 10691 - 683222859828090890299791032836..667114601633276209259644365483 (3218 digits)
27: 11279 - 692149388152358271880180910368..100257387805213081084271176363 (3395 digits)
28: 12391 - 385083595083686155852404805983..016288359843183967789680077483 (3730 digits)
29: 14479 - 136831460986221458568512534755..409643689259682678801865288363 (4359 digits)
30: 42737 - 438332262203687089339145725663..235631349264841566953606392491 (12865 digits)```

## Phix

```with javascript_semantics
atom t0 = time()
include mpfr.e

function isWagstaff(integer p, bool bLenOnly=false)
mpz w = mpz_init()
mpz_ui_pow_ui(w,2,p)
mpz_add_si(w,w,1)
assert(mpz_fdiv_q_ui(w,w,3)=0)
return iff(mpz_prime(w)?iff(bLenOnly?mpz_sizeinbase(w,10):mpz_get_str(w,comma_fill:=true)):"")
end function

integer p, pdx = 2
sequence wagstaff = {}
while length(wagstaff)<10 do
p = get_prime(pdx)
string res = isWagstaff(p)
if length(res) then wagstaff &= {{p,res}} end if
pdx += 1
end while
printf(1,"First 10 Wagstaff primes for the values of 'p' shown:\n")
papply(true,printf,{1,{"%2d: %s\n"},wagstaff})
printf(1,"\nTook %s\n\n",elapsed(time()-t0))

printf(1,"Wagstaff primes that we can find in 5 seconds:\n")
while time()<t0+5 do
p = get_prime(pdx)
object res = isWagstaff(p,true)
if integer(res) then
printf(1,"%5d (%,d digits, %s)\n", {p,res,elapsed(time()-t0)})
end if
pdx += 1
end while
```
Output:
```First 10 Wagstaff primes for the values of 'p' shown:
3: 3
5: 11
7: 43
11: 683
13: 2,731
17: 43,691
19: 174,763
23: 2,796,203
31: 715,827,883
43: 2,932,031,007,403

Took 0.1s

Wagstaff primes that we can find in 5 seconds:
61 (18 digits, 0.1s)
79 (24 digits, 0.1s)
101 (30 digits, 0.1s)
127 (38 digits, 0.1s)
167 (50 digits, 0.1s)
191 (58 digits, 0.1s)
199 (60 digits, 0.1s)
313 (94 digits, 0.1s)
347 (104 digits, 0.1s)
701 (211 digits, 0.2s)
1709 (514 digits, 1.1s)
2617 (788 digits, 4.4s)
```

## Python

```""" Rosetta code Wagstaff_primes task """

from sympy import isprime

def wagstaff(N):
""" find first N Wagstaff primes """
pri, wcount = 1, 0
while wcount < N:
pri += 2
if isprime(pri):
wag = (2**pri + 1) // 3
if isprime(wag):
wcount += 1
print(f'{wcount: 3}: {pri: 5} => ',
f'{wag:,}' if wcount < 11 else f'[{len(str(wag))} digit number]')

wagstaff(24)
```
Output:
```  1:     3 =>  3
2:     5 =>  11
3:     7 =>  43
4:    11 =>  683
5:    13 =>  2,731
6:    17 =>  43,691
7:    19 =>  174,763
8:    23 =>  2,796,203
9:    31 =>  715,827,883
10:    43 =>  2,932,031,007,403
11:    61 =>  [18 digit number]
12:    79 =>  [24 digit number]
13:   101 =>  [30 digit number]
14:   127 =>  [38 digit number]
15:   167 =>  [50 digit number]
16:   191 =>  [58 digit number]
17:   199 =>  [60 digit number]
18:   313 =>  [94 digit number]
19:   347 =>  [104 digit number]
20:   701 =>  [211 digit number]
21:  1709 =>  [514 digit number]
22:  2617 =>  [788 digit number]
23:  3539 =>  [1065 digit number]
24:  5807 =>  [1748 digit number]
```

## Quackery

`from`, `index`, `incr`, and `end` are defined at Loops/Increment loop index within loop body#Quackery.

`prime` is defined at Miller–Rabin primality test#Quackery.

```  [ bit 1+ 3 / ] is wagstaff ( n --> n )

[] 3 from
[ 2 incr
index prime while
index wagstaff prime while
index join
dup size 24 = if end ]
10 split swap
witheach
[ dup say "p = " echo
wagstaff
say ", w(p) = " echo cr ]
witheach
[ say "p = " echo cr ]```
Output:
```p = 3, w(p) = 3
p = 5, w(p) = 11
p = 7, w(p) = 43
p = 11, w(p) = 683
p = 13, w(p) = 2731
p = 17, w(p) = 43691
p = 19, w(p) = 174763
p = 23, w(p) = 2796203
p = 31, w(p) = 715827883
p = 43, w(p) = 2932031007403
p = 61
p = 79
p = 101
p = 127
p = 167
p = 191
p = 199
p = 313
p = 347
p = 701
p = 1709
p = 2617
p = 3539
p = 5807```

## Raku

```# First 20

my @wagstaff = (^∞).grep: { .is-prime && ((1 + 1 +< \$_)/3).is-prime };

say ++\$ ~ ": \$_ - {(1 + 1 +< \$_)/3}" for @wagstaff[^20];

say .fmt("\nTotal elapsed seconds: (%.2f)\n") given (my \$elapsed = now) - INIT now;

# However many I have patience for

my atomicint \$count = 20;

hyper for @wagstaff[20] .. * {
next unless .is-prime;
say ++⚛\$count ~ ": \$_ ({sprintf "%.2f", now - \$elapsed})" and \$elapsed = now if is-prime (1 + 1 +< \$_)/3;
}
```
Output:

Turns out, "however many I have patience for" is 9 (29).

```1: 3 - 3
2: 5 - 11
3: 7 - 43
4: 11 - 683
5: 13 - 2731
6: 17 - 43691
7: 19 - 174763
8: 23 - 2796203
9: 31 - 715827883
10: 43 - 2932031007403
11: 61 - 768614336404564651
12: 79 - 201487636602438195784363
13: 101 - 845100400152152934331135470251
14: 127 - 56713727820156410577229101238628035243
15: 167 - 62357403192785191176690552862561408838653121833643
16: 191 - 1046183622564446793972631570534611069350392574077339085483
17: 199 - 267823007376498379256993682056860433753700498963798805883563
18: 313 - 5562466239377370006237035693149875298444543026970449921737087520370363869220418099018130434731
19: 347 - 95562442332919646317117537304253622533190207882011713489066201641121786503686867002917439712921903606443
20: 701 - 3506757267698915671493993255253419110366893201882115906332187268712488102864720548075777690851725634151321979237452145142012954833455869242085209847101253899970149114085629732127775838589548478180862167823854251

Total elapsed seconds: (0.09)

21: 1709 (0.87)
22: 2617 (0.92)
23: 3539 (2.03)
24: 5807 (13.18)
25: 10501 (114.14)
26: 10691 (114.14)
27: 11279 (48.13)
28: 12391 (92.99)
29: 14479 (179.87)
^C```

## RPL

Works with: HP version 49
```≪ { } 2
WHILE OVER SIZE 10 < REPEAT
NEXTPRIME
2 OVER ^ 1 + 3 /
IF DUP ISPRIME? THEN
OVER →STR ":" + SWAP + ROT SWAP + SWAP
ELSE DROP END
END DROP
≫ 'WAGSTAFF' STO
```
Output:
```1: { "3:3" "5:11" "7:43" "11:683" "13:2731" "17:43691" "19:174763" "23:2796203" "31:715827883" "43:2932031007403" }
```

## Ruby

For the first 12 the Ruby prime library is sufficiently fast. Here, the external GMP gem is used in order to do more.

```require 'prime'
require 'gmp'

wagstaffs = Enumerator.new do |y|
odd_primes = Prime.each
odd_primes.next #skip 2
loop do
p = odd_primes.next
candidate = (2 ** p + 1)/3
y << [p, candidate] unless GMP::Z.new(candidate).probab_prime?.zero?
end
end

10.times{puts "%5d - %s" % wagstaffs.next}
14.times{puts "%5d" % wagstaffs.next.first}
```
Output:
```    3 - 3
5 - 11
7 - 43
11 - 683
13 - 2731
17 - 43691
19 - 174763
23 - 2796203
31 - 715827883
43 - 2932031007403
61
79
101
127
167
191
199
313
347
701
1709
2617
3539
5807
```

## Wren

Library: Wren-math
Library: Wren-gmp
Library: Wren-fmt

An embedded solution so we can use GMP to test for probable primeness in a reasonable time.

Even so, as far as the stretch goal is concerned, takes around 25 seconds to find the next 14 Wagstaff primes but almost 10 minutes to find the next 19 on my machine (Core i7). I gave up after that.

```import "./math" for Int
import "./gmp" for Mpz
import "./fmt" for Fmt

var isWagstaff = Fn.new { |p|
var w = (2.pow(p) + 1) / 3  // always integral
if (!Int.isPrime(w)) return [false, null]
return [true, [p, w]]
}

var isBigWagstaff = Fn.new { |p|
var w = Mpz.one.lsh(p).add(1).div(3)
return w.probPrime(15) > 0
}

var start = System.clock
var p = 1
var wagstaff = []
while (wagstaff.count < 10) {
while (true) {
p = p + 2
if (Int.isPrime(p)) break
}
var res = isWagstaff.call(p)
if (res[0]) wagstaff.add(res[1])
}
System.print("First 10 Wagstaff primes for the values of 'p' shown:")
for (i in 0..9) Fmt.print("\$2d: \$d", wagstaff[i][0], wagstaff[i][1])
System.print("\nTook %(System.clock - start) seconds")

var limit = 19
var count = 0
System.print("\nValues of 'p' for the next %(limit) Wagstaff primes and")
System.print("overall time taken to reach them to higher second:")
while (count < limit) {
while (true) {
p = p + 2
if (Int.isPrime(p)) break
}
if (isBigWagstaff.call(p)) {
Fmt.print("\$5d (\$3d secs)", p, (System.clock - start).ceil)
count = count + 1
}
}
```
Output:
```First 10 Wagstaff primes for the values of 'p' shown:
3: 3
5: 11
7: 43
11: 683
13: 2731
17: 43691
19: 174763
23: 2796203
31: 715827883
43: 2932031007403

Took 0.055678 seconds

Values of 'p' for the next 19 Wagstaff primes and
overall time taken to reach them to higher second:
61 (  1 secs)
79 (  1 secs)
101 (  1 secs)
127 (  1 secs)
167 (  1 secs)
191 (  1 secs)
199 (  1 secs)
313 (  1 secs)
347 (  1 secs)
701 (  1 secs)
1709 (  1 secs)
2617 (  2 secs)
3539 (  5 secs)
5807 ( 25 secs)
10501 (193 secs)
10691 (205 secs)
11279 (244 secs)
12391 (341 secs)
14479 (593 secs)
```

## XPL0

```func IsPrime(N);        \Return 'true' if N is prime
real N;  int  I;
[if N <= 2. then return N = 2.;
if Mod(N, 2.) = 0. then \even\ return false;
for I:= 3 to fix(sqrt(N)) do
[if Mod(N, float(I)) = 0. then return false;
I:= I+1;
];
return true;
];

real P, Q;  int C;
[P:= 2.;  C:= 0;
Format(1, 0);
repeat  if IsPrime(P) then
[Q:= Pow(2., P) + 1.;
if Mod(Q, 3.) = 0. and IsPrime(Q/3.) then
[Text(0, "(2^^");
RlOut(0, P);
Text(0, " + 1)/3 = ");
RlOut(0, Q/3.);
CrLf(0);
C:= C+1;
];
];
P:= P+1.;
until   C >= 10;
]```
Output:
```(2^3 + 1)/3 = 3
(2^5 + 1)/3 = 11
(2^7 + 1)/3 = 43
(2^11 + 1)/3 = 683
(2^13 + 1)/3 = 2731
(2^17 + 1)/3 = 43691
(2^19 + 1)/3 = 174763
(2^23 + 1)/3 = 2796203
(2^31 + 1)/3 = 715827883
(2^43 + 1)/3 = 2932031007403
```