I'm working on modernizing Rosetta Code's infrastructure. Starting with communications. Please accept this time-limited open invite to RC's Slack.. --Michael Mol (talk) 20:59, 30 May 2020 (UTC)

Wieferich primes

From Rosetta Code
Task
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.


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


Ada[edit]

Translation of: BASIC256
with Ada.Text_IO;
 
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
 
Ada.Text_IO.Put_Line ("Wieferich primes below 5000:");
for N in 1 .. 4999 loop
if Is_Wieferich (N) then
Ada.Text_IO.Put_Line (N'Image);
end if;
end loop;
 
end Wieferich_Primes;
Output:
Wieferich primes below 5000:
 1093
 3511

APL[edit]

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[edit]

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[edit]

 
# 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[edit]

BASIC256[edit]

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[edit]

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[edit]

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[edit]

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[edit]

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++[edit]

#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#[edit]

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) {
result.Add(p);
}
}
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

F#[edit]

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[edit]

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[edit]

 
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[edit]

Works with: Gforth
: prime? ( n -- ? ) here + [email protected] 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[edit]

 
#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[edit]

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

Haskell[edit]

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[edit]

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

About 12 times faster:

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

Java[edit]

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)
result.add(p);
}
return result;
}
}
Output:
Wieferich primes less than 5000:
1093
3511

jq[edit]

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

The task

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

Julia[edit]

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

Mathematica/Wolfram Language[edit]

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

Nim[edit]

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[edit]

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[edit]

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[edit]

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[edit]

(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[edit]

#!/usr/bin/python
 
def isPrime(n):
for i in range(2, int(n**0.5) + 1):
if n % i == 0:
return False
return True
 
def isWeiferich(p):
if not isPrime(p):
return False
q = 1
p2 = p ** 2
while p > 1:
q = (2 * q) % p2
p -= 1
if q == 1:
return True
else:
return False
 
 
if __name__ == '__main__':
print("Wieferich primes less than 5000: ")
for i in range(2, 5001):
if isWeiferich(i):
print(i)
Output:
Wieferich primes less than 5000: 
1093
3511

Quackery[edit]

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[edit]

#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[edit]

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[edit]

/*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 [email protected].#+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

Ruby[edit]

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

Rust[edit]

// [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[edit]

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[edit]

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[edit]

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