Wagstaff primes
A Wagstaff prime is a prime number of the form (2^p + 1)/3 where the exponent p is an odd prime.
Wagstaff 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.
- Definition
- 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.
- 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
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
Wren
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
if (!w.isInteger || !Int.isPrime(w)) return [false, null]
return [true, [p, w]]
}
var isBigWagstaff = Fn.new { |p|
var w = Mpz.one.lsh(p).add(1)
if (!w.isDivisibleUi(3)) return [false, null]
w.div(3)
if (w.probPrime(15) == 0) return [false, null]
return [true, p]
}
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
}
var res = isBigWagstaff.call(p)
if (res[0]) {
Fmt.print("$5d ($3d secs)", res[1], (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.056335 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 (198 secs) 10691 (210 secs) 11279 (251 secs) 12391 (348 secs) 14479 (599 secs)