Wagstaff primes: Difference between revisions
(Wagstaff primes in FreeBASIC) |
(Wagstaff primes in BASIC-256) |
||
Line 82:
10: 43: 2932031007403
</pre>
=={{header|BASIC256}}==
{{trans|FreeBASIC}}
<syntaxhighlight lang="BASIC256">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</syntaxhighlight>
=={{header|C++}}==
|
Revision as of 18:46, 28 September 2022
- 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
BASIC256
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++
#include <gmpxx.h>
#include <primesieve.hpp>
#include <iostream>
using big_int = mpz_class;
std::string to_string(const big_int& num, size_t n) {
std::string str = num.get_str();
size_t len = str.size();
if (len > n) {
str = str.substr(0, n / 2) + "..." + str.substr(len - n / 2);
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)
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
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
J
x:(,.f)p:I.1 p:(f=.3%~1+2^])p:i.14
3 3
5 11
7 43
11 683
13 2731
17 43691
19 174763
23 2796203
31 715827883
43 2932031007403
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
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...
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]
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 // 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)