# Wieferich primes

(Redirected from Weiferich primes)
Wieferich primes
You are encouraged to solve this task according to the task description, using any language you may know.
 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.

• 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.

Translation of: BASIC256

procedure Wieferich_Primes is

function Is_Prime (V : Positive) return Boolean is
D : Positive := 5;
begin
if V < 2       then return False; end if;
if V mod 2 = 0 then return V = 2; end if;
if V mod 3 = 0 then return V = 3; end if;
while D * D <= V loop
if V mod D = 0 then
return False;
end if;
D := D + 2;
end loop;
return True;
end Is_Prime;

function Is_Wieferich (N : Positive) return Boolean is
Q : Natural := 1;
begin
if not Is_Prime (N) then
return False;
end if;
for P in 2 .. N loop
Q := (2 * Q) mod N**2;
end loop;
return Q = 1;
end Is_Wieferich;

begin

for N in 1 .. 4999 loop
if Is_Wieferich (N) then
end if;
end loop;

end Wieferich_Primes;
Output:
Wieferich primes below 5000:
1093
3511

## ALGOL 68

Works with: ALGOL 68G version Any - tested with release 2.8.3.win32
BEGIN # find Wierferich Primes: primes p where p^2 evenly divides 2^(p-1)-1  #

INT max number = 5 000; # maximum number we will consider                #
# set precision of LONG LONG INT - p^5000 has over 1500 digits           #
PR precision 1600 PR
PR read "primes.incl.a68" PR                    # include prime utlities #
# get a list of primes up to max number                                  #
[]INT prime = EXTRACTPRIMESUPTO max number
FROMPRIMESIEVE PRIMESIEVE max number;

# find the primes                                                        #
INT           p pos            := LWB prime;
LONG LONG INT two to p minus 1 := 1;
INT           power            := 0;
INT           w count          := 0;
WHILE w count < 2 DO
INT p = prime[ p pos ];
WHILE power < ( p - 1 ) DO
two to p minus 1 *:= 2;
power            +:= 1
OD;
IF ( two to p minus 1 - 1 ) MOD ( p * p ) = 0 THEN
print( ( " ", whole( p, 0 ) ) );
w count +:= 1
FI;
p pos +:= 1
OD
END
Output:
1093 3511

## APL

Works in: Dyalog APL

⎕CY 'dfns' ⍝ import dfns namespace
⍝ pco ← prime finder
⍝ nats ← natural number arithmetic (uses strings)
⍝ Get all Wieferich primes below n:
wief{{/⍨{(,'0')(×)|nats 1 -nats 2 *nats -1}¨}1 pco}
wief 5000
1093 3511

## Arturo

wieferich?: function [n][
and? -> prime? n
-> zero? (dec 2 ^ n-1) % n ^ 2
]

print ["Wieferich primes less than 5000:" select 1..5000 => wieferich?]
Output:
Wieferich primes less than 5000: [1093 3511]

## AWK

# syntax: GAWK -f WIEFERICH_PRIMES.AWK
# converted from FreeBASIC
BEGIN {
start = 1
stop = 4999
for (i=start; i<=stop; i++) {
if (is_wieferich_prime(i)) {
printf("%d\n",i)
count++
}
}
printf("Wieferich primes %d-%d: %d\n",start,stop,count)
exit(0)
}
function is_prime(n,  d) {
d = 5
if (n < 2) { return(0) }
if (n % 2 == 0) { return(n == 2) }
if (n % 3 == 0) { return(n == 3) }
while (d*d <= n) {
if (n % d == 0) { return(0) }
d += 2
if (n % d == 0) { return(0) }
d += 4
}
return(1)
}
function is_wieferich_prime(p,  p2,q) {
if (!is_prime(p)) { return(0) }
q = 1
p2 = p^2
while (p > 1) {
q = (2*q) % p2
p--
}
return(q == 1)
}
Output:
1093
3511
Wieferich primes 1-4999: 2

## BASIC

### BASIC256

Translation of: FreeBASIC
print "Wieferich primes less than 5000: "
for i = 1 to 5000
if isWeiferich(i) then print i
next i
end

function isWeiferich(p)
if not isPrime(p) then return False
q = 1
p2 = p ^ 2
while p > 1
q = (2 * q) mod p2
p -= 1
end while
if q = 1 then return True else return False
end function

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
Output:
Igual que la entrada de FreeBASIC.

### PureBasic

Translation of: FreeBASIC
Procedure.i isPrime(n)
Protected k

If n = 2 : ProcedureReturn #True
ElseIf n <= 1 Or n % 2 = 0 : ProcedureReturn #False
Else
For k = 3 To Int(Sqr(n)) Step 2
If n % k = 0
ProcedureReturn #False
EndIf
Next
EndIf

ProcedureReturn #True
EndProcedure

Procedure.i isWeiferich(p)
Protected q, p2

If Not isPrime(p) : ProcedureReturn #False : EndIf
q = 1
p2 = Pow(p, 2)
While p > 1
q = (2*q) % p2
p - 1
Wend
If q = 1
ProcedureReturn #True
Else
ProcedureReturn #False
EndIf
EndProcedure

OpenConsole()
PrintN("Wieferich primes less than 5000: ")
For i = 2 To 5000
If isWeiferich(i)
PrintN(Str(i))
EndIf
Next i
Input()
CloseConsole()
Output:
Igual que la entrada de FreeBASIC.

### Run BASIC

print "Wieferich primes less than 5000: "
for i = 1 to 5000
if isWeiferich(i) then print i
next i
end

function isPrime(n)
if n < 2       then isPrime = 0 : goto [exit]
if n = 2       then isPrime = 1 : goto [exit]
if n mod 2 = 0 then isPrime = 0 : goto [exit]
isPrime = 1
for i = 3 to int(n^.5) step 2
if n mod i = 0 then isPrime = 0 : goto [exit]
next i
[exit]
end function

function isWeiferich(p)
if isPrime(p) = 0 then isWeiferich = 0 : goto [exit]
q = 1
p2 = p^2
while p > 1
q = (2*q) mod p2
p = p - 1
wend
if q = 1 then
isWeiferich = 1 : goto [exit]
else
isWeiferich = 0 : goto [exit]
end if
[exit]
end function
Output:
Igual que la entrada de FreeBASIC.

### Yabasic

Translation of: FreeBASIC
print "Wieferich primes less than 5000: "
for i = 2 to 5000
if isWeiferich(i)  print i
next i
end

sub isWeiferich(p)
if not isPrime(p)  return False
q = 1
p2 = p ^ 2
while p > 1
q = mod((2*q), p2)
p = p - 1
wend
if q = 1 then return True else return False : fi
end sub

sub isPrime(v)
if v < 2  return False
if mod(v, 2) = 0  return v = 2
if mod(v, 3) = 0  return v = 3
d = 5
while d * d <= v
if mod(v, d) = 0 then return False else d = d + 2 : fi
wend
return True
end sub
Output:
Igual que la entrada de FreeBASIC.

## C

Translation of: C++
#include <stdbool.h>
#include <stdio.h>
#include <stdint.h>

#define LIMIT 5000
static bool PRIMES[LIMIT];

static void prime_sieve() {
uint64_t p;
int i;

PRIMES[0] = false;
PRIMES[1] = false;
for (i = 2; i < LIMIT; i++) {
PRIMES[i] = true;
}

for (i = 4; i < LIMIT; i += 2) {
PRIMES[i] = false;
}

for (p = 3;; p += 2) {
uint64_t q = p * p;
if (q >= LIMIT) {
break;
}
if (PRIMES[p]) {
uint64_t inc = 2 * p;
for (; q < LIMIT; q += inc) {
PRIMES[q] = false;
}
}
}
}

uint64_t modpow(uint64_t base, uint64_t exp, uint64_t mod) {
uint64_t result = 1;

if (mod == 1) {
return 0;
}

base %= mod;
for (; exp > 0; exp >>= 1) {
if ((exp & 1) == 1) {
result = (result * base) % mod;
}
base = (base * base) % mod;
}
return result;
}

void wieferich_primes() {
uint64_t p;

for (p = 2; p < LIMIT; ++p) {
if (PRIMES[p] && modpow(2, p - 1, p * p) == 1) {
printf("%lld\n", p);
}
}
}

int main() {
prime_sieve();

printf("Wieferich primes less than %d:\n", LIMIT);
wieferich_primes();

return 0;
}
Output:
Wieferich primes less than 5000:
1093
3511

## C++

#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';
}
Output:
Wieferich primes less than 5000:
1093
3511

## C#

Translation of: Java
using System;
using System.Collections.Generic;
using System.Linq;

namespace WieferichPrimes {
class Program {
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;
}

static bool[] PrimeSieve(int limit) {
bool[] sieve = Enumerable.Repeat(true, limit).ToArray();

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;
}

static List<int> WiefreichPrimes(int limit) {
bool[] sieve = PrimeSieve(limit);
List<int> result = new List<int>();
for (int p = 2; p < limit; p++) {
if (sieve[p] && ModPow(2, p - 1, p * p) == 1) {
}
}
return result;
}

static void Main() {
const int limit = 5000;
Console.WriteLine("Wieferich primes less that {0}:", limit);
foreach (int p in WiefreichPrimes(limit)) {
Console.WriteLine(p);
}
}
}
}
Output:
Wieferich primes less that 5000:
1093
3511

## EasyLang

Translation of: BASIC256
fastfunc isprim num .
i = 2
while i <= sqrt num
if num mod i = 0
return 0
.
i += 1
.
return 1
.
func wieferich p .
if isprim p = 0
return 0
.
q = 1
p2 = p * p
while p > 1
q = (2 * q) mod p2
p -= 1
.
if q = 1
return 1
.
.
print "Wieferich primes less than 5000: "
for i = 2 to 5000
if wieferich i = 1
print i
.
.
Output:
Wieferich primes less than 5000:
1093
3511

## F#

This task uses Extensible Prime Generator (F#)

// 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")
Output:
1093
3511
Real: 00:00:00.004

## Factor

Works with: Factor version 0.99 2021-02-05
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 .
Output:
Wieferich primes less than 5000:
V{ 1093 3511 }

## fermat

Func Iswief(p)=Isprime(p)*Divides(p^2, 2^(p-1)-1).
for i=2 to 5000 do if Iswief(i) then !!i fi od
Output:

1093 3511

## Forth

Works with: Gforth
: 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
Output:
Wieferich primes less than 5000:
1093
3511

## FreeBASIC

#include "isprime.bas"

function iswief( byval p as uinteger ) as boolean
if not isprime(p) then return 0
dim as integer q = 1, p2 = p^2
while p>1
q=(2*q) mod p2
p = p - 1
wend
if q=1 then return 1 else return 0
end function

for i as uinteger = 1 to 5000
if iswief(i) then print i
next i

## Go

Translation of: Wren
Library: Go-rcu
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))
}
}
}
Output:
Wieferich primes < 5,000:
1,093
3,511

isPrime :: Integer -> Bool
isPrime n
|n == 2 = True
|n == 1 = False
|otherwise = null \$ filter (\i -> mod n i == 0 ) [2 .. root]
where
root :: Integer
root = toInteger \$ floor \$ sqrt \$ fromIntegral n

isWieferichPrime :: Integer -> Bool
isWieferichPrime n =  isPrime n && mod ( 2 ^ ( n - 1 ) - 1 ) ( n ^ 2 ) == 0

solution :: [Integer]
solution = filter isWieferichPrime [2 .. 5000]

main :: IO ( )
main = do
putStrLn "Wieferich primes less than 5000:"
print solution
Output:
Wieferich primes less than 5000:
[1093,3511]

## J

I.(1&p: * 0=*: | _1+2x^<:) i.5000
1093 3511

p: I. (0=*:|_1+2x^<:) I.1 p: i.5000
1093 3511

## Java

Translation of: C++
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)
}
return result;
}
}
Output:
Wieferich primes less than 5000:
1093
3511

## jq

Works with gojq, the Go implementation of jq

gojq supports unbounded-precision integer arithmetic and so is up to this task.

def is_prime:
. as \$n
| if (\$n < 2)         then false
elif (\$n % 2 == 0)  then \$n == 2
elif (\$n % 3 == 0)  then \$n == 3
elif (\$n % 5 == 0)  then \$n == 5
elif (\$n % 7 == 0)  then \$n == 7
elif (\$n % 11 == 0) then \$n == 11
elif (\$n % 13 == 0) then \$n == 13
elif (\$n % 17 == 0) then \$n == 17
elif (\$n % 19 == 0) then \$n == 19
else {i:23}
| until( (.i * .i) > \$n or (\$n % .i == 0); .i += 2)
| .i * .i > \$n
end;

# Emit an array of primes less than `.`
def primes:
if . < 2 then []
else
[2] + [range(3; .; 2) | select(is_prime)]
end;

# for the sake of infinite-precision integer arithmetic
def power(\$b): . as \$a | reduce range(0; \$b) as \$i (1; .*\$a);

# Input: the limit
def wieferich:
primes[]
| . as \$p
| select( ( (2|power(\$p-1)) - 1) % (.*.) == 0);

5000 | wieferich
Output:
1093
3511

## Julia

using Primes

println(filter(p -> (big"2"^(p - 1) - 1) % p^2 == 0, primes(5000)))  # [1093, 3511]

## Mathematica/Wolfram Language

ClearAll[WieferichPrimeQ]
WieferichPrimeQ[n_Integer] := PrimeQ[n] && Divisible[2^(n - 1) - 1, n^2]
Select[Range[5000], WieferichPrimeQ]
Output:
{1093, 3511}

## Nim

Library: bignum
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
Output:
Wieferich primes less than 5000:
1093
3511

## PARI/GP

iswief(p)=if(isprime(p)&&(2^(p-1)-1)%p^2==0,1,0)
for(N=1,5000,if(iswief(N),print(N)))
Output:
1093
3511

## Perl

Library: ntheory
use feature 'say';
use ntheory qw(is_prime powmod);

say 'Wieferich primes less than 5000: ' . join ', ', grep { is_prime(\$_) and powmod(2, \$_-1, \$_*\$_) == 1 } 1..5000;
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)})

## PicoLisp

(de **Mod (X Y N)
(let M 1
(loop
(when (bit? 1 Y)
(setq M (% (* M X) N)) )
(T (=0 (setq Y (>> 1 Y)))
M )
(setq X (% (* X X) N)) ) ) )
(let (D 2  L (1 2 2 . (4 2 4 2 4 6 2 6 .)))
(until (> D 5000)
(and
(=1 (**Mod 2 (dec D) (* D D)))
(println D) )
(inc 'D (++ L)) ) )
Output:
1093
3511

## Python

# Wieferich-Primzahlen
MAX: int = 5_000

# Berechnet a^n mod m
def pow_mod(a: int, n: int, m: int) -> int:
assert n >= 0 and m != 0, "pow_mod(a, n, m), n >= 0, m <> 0"
res: int = 1
a %= m
while n > 0:
if n%2:
res = (res*a)%m
n -= 1
else:
a = (a*a)%m
n //= 2
return res%m

def is_prime(n: int) -> bool:
for i in range(2, int(n**0.5) + 1):
if n%i == 0:
return False
return True

def is_wieferich(p: int) -> True:
if is_prime(p) == False:
return False
if pow_mod(2, p - 1, p*p) == 1:
return True
else:
return False

if __name__ == '__main__':
print(f"Wieferich primes less than {MAX}:")
for i in range(2, MAX + 1):
if is_wieferich(i):
print(i)
Output:
Wieferich primes less than 5000:
1093
3511

## Quackery

eratosthenes and isprime are defined at Sieve of Eratosthenes#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 ] ]
Output:
1093
3511

## Racket

#lang typed/racket
(require math/number-theory)

(: wieferich-prime? (-> Positive-Integer Boolean))

(define (wieferich-prime? p)
(and (prime? p)
(divides? (* p p) (sub1 (expt 2 (sub1 p))))))

(module+ main
(define wieferich-primes<5000
(for/list : (Listof Integer) ((p (sequence-filter wieferich-prime?
(in-range 1 5000))))
p))
wieferich-primes<5000)
Output:
'(1093 3511)

## Raku

put "Wieferich primes less than 5000: ", join ', ', ^5000 .grep: { .is-prime and not ( exp(\$_-1, 2) - 1 ) % .² };
Output:
Wieferich primes less than 5000: 1093, 3511

## 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   "   "     "       */
found= 0                                         /*initialize number of Wieferich primes*/
do j=1  to #;    p= @.j                 /*search for Wieferich primes in range.*/
if (2**(p-1)-1)//p**2\==0  then iterate /*P**2 not evenly divide  2**(P-1) - 1?*/       /* ◄■■■■■■■ the filter.*/
found= found + 1                        /*bump the counter of Wieferich primes.*/
say center(found, 7)'│'  center(commas(p), w)    /*display the Wieferich prime.*/
end   /*j*/

say '───────┴'center(""   , w, '─');   say       /*display a  foot sep  for the output. */
say 'Found '       commas(found)       title     /*   "    "  summary    "   "     "    */
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:        @.1=2; @.2=3; @.3=5; @.4=7; @.5=11  /*define some low primes  (index-1).   */
!.=0;  !.2=1; !.3=1; !.5=1; !.7=1; !.11=1  /*   "     "   "    "     (semaphores).*/
#= 5;  sq.#= @.# ** 2 /*number of primes so far;     prime². */
do j=@.#+2  by 2  to n-1;  parse var j '' -1 _   /*get right decimal digit of J.*/
if    _==5  then iterate                               /*J ÷ by 5?    Yes, skip.*/
if j//3==0  then iterate;  if j//7==0  then iterate    /*" "  " 3?    J ÷ by 7? */
do k=5  while sq.k<=j             /* [↓]  divide by the known odd primes.*/
if j//@.k==0  then iterate j      /*Is J ÷ a P?   Then not prime.   ___  */
end   /*k*/                       /* [↑]  only process numbers  ≤  √ J   */
#= #+1;    @.#= j;   sq.#= j*j;   !.j= 1 /*bump # Ps; assign next P; P sqare; P.*/
end          /*j*/;               return
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

## RPL

Works with: RPL version HP49-C
« { } 2
WHILE DUP 5000 < REPEAT
IF 2 OVER 1 - ^ 1 - OVER SQ MOD NOT THEN SWAP OVER + SWAP END
NEXTPRIME
END DROP
Output:
1: { 1093 3511 }

## Ruby

require "prime"

puts Prime.each(5000).select{|p| 2.pow(p-1 ,p*p) == 1 }
Output:
1093
3511

## 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);
}
}
Output:
Wieferich primes less than 5000:
1093
3511

## Sidef

func is_wieferich_prime(p, base=2) {
powmod(base, p-1, p**2) == 1
}

say ("Wieferich primes less than 5000: ", 5000.primes.grep(is_wieferich_prime))
Output:
Wieferich primes less than 5000: [1093, 3511]

## Swift

Translation of: C++
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)
}
Output:
Wieferich primes less than 5000:
1093
3511

## Wren

Library: Wren-math
Library: Wren-big
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)
}
Output:
Wieferich primes < 5000:
1093
3511