Wieferich primes
| This page uses content from Wikipedia. The original article was at Wieferich prime. The list of authors can be seen in the page history. As with Rosetta Code, the text of Wikipedia is available under the GNU FDL. (See links for details on variance) |
In number theory, a Wieferich prime is a prime number p such that p2 evenly divides 2(p − 1) − 1 .
It is conjectured that there are infinitely many Wieferich primes, but as of March 2021,only two have been identified.
- Task
- Write a routine (function procedure, whatever) to find Wieferich primes.
- Use that routine to identify and display all of the Wieferich primes less than 5000.
- See also
C++
<lang cpp>#include <cstdint>
- include <iostream>
- include <vector>
std::vector<bool> prime_sieve(uint64_t limit) {
std::vector<bool> sieve(limit, true);
if (limit > 0)
sieve[0] = false;
if (limit > 1)
sieve[1] = false;
for (uint64_t i = 4; i < limit; i += 2)
sieve[i] = false;
for (uint64_t p = 3; ; p += 2) {
uint64_t q = p * p;
if (q >= limit)
break;
if (sieve[p]) {
uint64_t inc = 2 * p;
for (; q < limit; q += inc)
sieve[q] = false;
}
}
return sieve;
}
uint64_t modpow(uint64_t base, uint64_t exp, uint64_t mod) {
if (mod == 1)
return 0;
uint64_t result = 1;
base %= mod;
for (; exp > 0; exp >>= 1) {
if ((exp & 1) == 1)
result = (result * base) % mod;
base = (base * base) % mod;
}
return result;
}
std::vector<uint64_t> wieferich_primes(uint64_t limit) {
std::vector<uint64_t> result;
std::vector<bool> sieve(prime_sieve(limit));
for (uint64_t p = 2; p < limit; ++p)
if (sieve[p] && modpow(2, p - 1, p * p) == 1)
result.push_back(p);
return result;
}
int main() {
const uint64_t limit = 5000;
std::cout << "Wieferich primes less than " << limit << ":\n";
for (uint64_t p : wieferich_primes(limit))
std::cout << p << '\n';
}</lang>
- Output:
Wieferich primes less than 5000: 1093 3511
F#
This task uses Extensible Prime Generator (F#) <lang fsharp> // Weiferich primes: Nigel Galloway. June 2nd., 2021 primes32()|>Seq.takeWhile((>)5000)|>Seq.filter(fun n->(2I**(n-1)-1I)%(bigint(n*n))=0I)|>Seq.iter(printfn "%d") </lang>
- Output:
1093 3511 Real: 00:00:00.004
Factor
<lang factor>USING: io kernel math math.functions math.primes prettyprint sequences ;
"Wieferich primes less than 5000:" print 5000 primes-upto [ [ 1 - 2^ 1 - ] [ sq divisor? ] bi ] filter .</lang>
- Output:
Wieferich primes less than 5000:
V{ 1093 3511 }
Forth
<lang forth>: prime? ( n -- ? ) here + c@ 0= ;
- notprime! ( n -- ) here + 1 swap c! ;
- prime_sieve { n -- }
here n erase
0 notprime!
1 notprime!
n 4 > if
n 4 do i notprime! 2 +loop
then
3
begin
dup dup * n <
while
dup prime? if
n over dup * do
i notprime!
dup 2* +loop
then
2 +
repeat
drop ;
- modpow { c b a -- a^b mod c }
c 1 = if 0 exit then
1
a c mod to a
begin
b 0>
while
b 1 and 1 = if
a * c mod
then
a a * c mod to a
b 2/ to b
repeat ;
- wieferich_prime? { p -- ? }
p prime? if p p * p 1- 2 modpow 1 = else false then ;
- wieferich_primes { n -- }
." Wieferich primes less than " n 1 .r ." :" cr
n prime_sieve
n 0 do
i wieferich_prime? if
i 1 .r cr
then
loop ;
5000 wieferich_primes bye</lang>
- Output:
Wieferich primes less than 5000: 1093 3511
Go
<lang go>package main
import (
"fmt" "math/big" "rcu"
)
func main() {
primes := rcu.Primes(5000)
zero := new(big.Int)
one := big.NewInt(1)
num := new(big.Int)
fmt.Println("Wieferich primes < 5,000:")
for _, p := range primes {
num.Set(one)
num.Lsh(num, uint(p-1))
num.Sub(num, one)
den := big.NewInt(int64(p * p))
if num.Rem(num, den).Cmp(zero) == 0 {
fmt.Println(rcu.Commatize(p))
}
}
}</lang>
- Output:
Wieferich primes < 5,000: 1,093 3,511
Java
<lang java>import java.util.*;
public class WieferichPrimes {
public static void main(String[] args) {
final int limit = 5000;
System.out.printf("Wieferich primes less than %d:\n", limit);
for (Integer p : wieferichPrimes(limit))
System.out.println(p);
}
private static boolean[] primeSieve(int limit) {
boolean[] sieve = new boolean[limit];
Arrays.fill(sieve, true);
if (limit > 0)
sieve[0] = false;
if (limit > 1)
sieve[1] = false;
for (int i = 4; i < limit; i += 2)
sieve[i] = false;
for (int p = 3; ; p += 2) {
int q = p * p;
if (q >= limit)
break;
if (sieve[p]) {
int inc = 2 * p;
for (; q < limit; q += inc)
sieve[q] = false;
}
}
return sieve;
}
private static long modpow(long base, long exp, long mod) {
if (mod == 1)
return 0;
long result = 1;
base %= mod;
for (; exp > 0; exp >>= 1) {
if ((exp & 1) == 1)
result = (result * base) % mod;
base = (base * base) % mod;
}
return result;
}
private static List<Integer> wieferichPrimes(int limit) {
boolean[] sieve = primeSieve(limit);
List<Integer> result = new ArrayList<>();
for (int p = 2; p < limit; ++p) {
if (sieve[p] && modpow(2, p - 1, p * p) == 1)
result.add(p);
}
return result;
}
}</lang>
- Output:
Wieferich primes less than 5000: 1093 3511
Julia
<lang julia>using Primes
println(filter(p -> (big"2"^(p - 1) - 1) % p^2 == 0, primes(5000))) # [1093, 3511] </lang>
Nim
<lang Nim>import math import bignum
func isPrime(n: Positive): bool =
if n mod 2 == 0: return n == 2 if n mod 3 == 0: return n == 3 var d = 5 while d <= sqrt(n.toFloat).int: if n mod d == 0: return false inc d, 2 if n mod d == 0: return false inc d, 4 result = true
echo "Wieferich primes less than 5000:" let two = newInt(2) for p in 2u..<5000:
if p.isPrime:
if exp(two, p - 1, p * p) == 1: # Modular exponentiation.
echo p</lang>
- Output:
Wieferich primes less than 5000: 1093 3511
Perl
<lang perl>use feature 'say'; use bignum; use ntheory 'is_prime';
say 'Wieferich primes less than 5000: ' . join ', ', grep { is_prime($_) and not ( (2**($_-1) -1) % $_**2 ) } 1..5000;</lang>
- Output:
Wieferich primes less than 5000: 1093, 3511
Phix
with javascript_semantics include mpfr.e function weiferich(integer p) mpz p2pm1m1 = mpz_init() mpz_ui_pow_ui(p2pm1m1,2,p-1) mpz_sub_ui(p2pm1m1,p2pm1m1,1) return mpz_fdiv_q_ui(p2pm1m1,p2pm1m1,p*p)=0 end function printf(1,"Weiferich primes less than 5000: %V\n",{filter(get_primes_le(5000),weiferich)})
- Output:
Wieferich primes less than 5000: {1093,3511}
alternative (same results), should be significantly faster, in the (largely pointless!) hunt for larger numbers.
with javascript_semantics include mpfr.e mpz base = mpz_init(2), {modulus, z} = mpz_inits(2) function weiferich(integer p) mpz_set_si(modulus,p*p) mpz_powm_ui(z, base, p-1, modulus) return mpz_cmp_si(z,1)=0 end function printf(1,"Weiferich primes less than 5000: %V\n",{filter(get_primes_le(5000),weiferich)})
Quackery
eratosthenes and isprime are defined at Sieve of Eratosthenes#Quackery.
<lang Quackery> 5000 eratosthenes
[ dup isprime iff
[ dup 1 - bit 1 -
swap dup * mod
0 = ]
else [ drop false ] ] is wieferich ( n --> b )
5000 times [ i^ wieferich if [ i^ echo cr ] ]</lang>
- Output:
1093 3511
Raku
<lang perl6>put "Wieferich primes less than 5000: ", join ', ', ^5000 .grep: { .is-prime and not ( exp($_-1, 2) - 1 ) % .² };</lang>
- Output:
Wieferich primes less than 5000: 1093, 3511
REXX
<lang rexx>/*REXX program finds and displays Wieferich primes which are under a specified limit N*/ parse arg n . /*obtain optional argument from the CL.*/ if n== | n=="," then n= 5000 /*Not specified? Then use the default.*/ numeric digits 3000 /*bump # of dec. digs for calculation. */ numeric digits max(9, length(2**n) ) /*calculate # of decimal digits needed.*/ call genP /*build array of semaphores for primes.*/ title= ' Wieferich primes that are < ' commas(n) /*title for the output. */ w= length(title) + 2 /*width of field for the primes listed.*/ say ' index │'center(title, w) /*display the title for the output. */ say '───────┼'center("" , w, '─') /* " a sep for the output. */ wp= 0 /*initialize number of Wieferich primes*/
do j=1 to #; p= @.j; pm= p - 1 /*search for Wieferich primes in range.*/
if (2**pm - 1) // p**2\==0 then iterate /*P**2 not evenly divide 2**(P-1) - 1?*/
wp= wp + 1 /*bump the counter of Wieferich primes.*/
say center(wp, 7)'│' center(commas(p), w) /*display the Wieferich prime to term. */
end /*j*/
say '───────┴'center("" , w, '─') /*display a foot sep for the output. */ say say 'Found ' commas(wp) title exit 0 /*stick a fork in it, we're all done. */ /*──────────────────────────────────────────────────────────────────────────────────────*/ commas: parse arg ?; do jc=length(?)-3 to 1 by -3; ?=insert(',', ?, jc); end; return ? /*──────────────────────────────────────────────────────────────────────────────────────*/ genP: !.= 0 /*placeholders for primes (semaphores).*/
@.1=2; @.2=3; @.3=5; @.4=7; @.5=11 /*define some low primes. */
!.2=1; !.3=1; !.5=1; !.7=1; !.11=1 /* " " " " flags. */
#=5; s.#= @.# **2 /*number of primes so far; prime². */
/* [↓] generate more primes ≤ high.*/
do j=@.#+2 by 2 to n-1 /*find odd primes from here on. */
parse var j -1 _; if _==5 then iterate /*J divisible by 5? (right dig)*/
if j// 3==0 then iterate /*" " " 3? */
if j// 7==0 then iterate /*" " " 7? */
/* [↑] the above five lines saves time*/
do k=5 while s.k<=j /* [↓] divide by the known odd primes.*/
if j // @.k == 0 then iterate j /*Is J ÷ X? Then not prime. ___ */
end /*k*/ /* [↑] only process numbers ≤ √ J */
#= #+1; @.#= j; s.#= j*j; !.j= 1 /*bump # of Ps; assign next P; P²; P# */
end /*j*/; return</lang>
- output when using the default input:
index │ Wieferich primes that are < 5,000 ───────┼────────────────────────────────────── 1 │ 1,093 2 │ 3,511 ───────┴────────────────────────────────────── Found 2 Wieferich primes that are < 5,000
Rust
<lang rust>// [dependencies] // primal = "0.3" // mod_exp = "1.0"
fn wieferich_primes(limit: usize) -> impl std::iter::Iterator<Item = usize> {
primal::Primes::all()
.take_while(move |x| *x < limit)
.filter(|x| mod_exp::mod_exp(2, *x - 1, *x * *x) == 1)
}
fn main() {
let limit = 5000;
println!("Wieferich primes less than {}:", limit);
for p in wieferich_primes(limit) {
println!("{}", p);
}
}</lang>
- Output:
Wieferich primes less than 5000: 1093 3511
Swift
<lang swift>func primeSieve(limit: Int) -> [Bool] {
guard limit > 0 else {
return []
}
var sieve = Array(repeating: true, count: limit)
sieve[0] = false
if limit > 1 {
sieve[1] = false
}
if limit > 4 {
for i in stride(from: 4, to: limit, by: 2) {
sieve[i] = false
}
}
var p = 3
while true {
var q = p * p
if q >= limit {
break
}
if sieve[p] {
let inc = 2 * p
while q < limit {
sieve[q] = false
q += inc
}
}
p += 2
}
return sieve
}
func modpow(base: Int, exponent: Int, mod: Int) -> Int {
if mod == 1 {
return 0
}
var result = 1
var exp = exponent
var b = base
b %= mod
while exp > 0 {
if (exp & 1) == 1 {
result = (result * b) % mod
}
b = (b * b) % mod
exp >>= 1
}
return result
}
func wieferichPrimes(limit: Int) -> [Int] {
let sieve = primeSieve(limit: limit)
var result: [Int] = []
for p in 2..<limit {
if sieve[p] && modpow(base: 2, exponent: p - 1, mod: p * p) == 1 {
result.append(p)
}
}
return result
}
let limit = 5000 print("Wieferich primes less than \(limit):") for p in wieferichPrimes(limit: limit) {
print(p)
}</lang>
- Output:
Wieferich primes less than 5000: 1093 3511
Wren
<lang ecmascript>import "/math" for Int import "/big" for BigInt
var primes = Int.primeSieve(5000) System.print("Wieferich primes < 5000:") for (p in primes) {
var num = (BigInt.one << (p - 1)) - 1 var den = p * p if (num % den == 0) System.print(p)
}</lang>
- Output:
Wieferich primes < 5000: 1093 3511