Wieferich primes

From Rosetta Code
(Redirected from Weiferich primes)
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

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

ALGOL 68

Works with: ALGOL 68G version Any - tested with release 2.8.3.win32
Library: ALGOL 68-primes
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) {
                    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

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

Haskell 

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

About 12 times faster: 

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

The task 

# 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 5000
Numeric Digits max(9,length(2**n))               /*calculate nr of decimal digits needed*/
call genP                                        /*build array of semaphores for primes.*/
/****************
Do i=1 To nn
  Say i pr.i sq.i
  End
****************/
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 nn
  p=pr.j                                         /*search for Wieferich primes in range.*/
  If (2**(p-1)-1)//p**2=0 Then Do                /* P**2 evenly divides  2**(P-1)-1     */
    found=found+1                                /*bump the counter of Wieferich primes.*/
    Say center(found,7)'|' center(commas(p),w)   /*display the Wieferich prime.*/
    End
  End

Say '--------'center(""   , w, '-')              /*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:
/*****************************************************************************************
* Compute all primes less than n+1
* Output: nn - the number of primes found
*         pr.i The i-th prime
*         sq.i the square of pr.i
*****************************************************************************************/
  primes='2 3 5 7 11'
  nn=0
  Do while primes>''                             /* for efficiency note a few primes    */
    Parse Var primes w primes
    Call store w
    End
  Do j=pr.nn+2 by 2 To n                         /* look at odd numbers up to n         */
    If right(j,1)<>5 Then                        /* number does not end in 5            */
       If j//3<>0 Then                           /*    "   is not a multiple of 3       */
         If j//7<>0 Then Do                      /*    "   is not a multiple of 7       */
           Do k=5 while sq.k<=j                  /* check all primes up to sqrt(j)      */
             If j//pr.k==0 Then iterate j        /* if this is a divisor: j is not prime*/
             End
           Call store j                          /* j is a prime number. store it       */
           End
    End
  Return

store:
  Parse Arg w
  nn=nn+1
  pr.nn=w
  sq.nn=w**2
  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
» 'TASK' STO
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
Retrieved from "https://rosettacode.org/wiki/Wieferich_primes?oldid=371360"