Wieferich primes
Appearance
Wieferich primes
You are encouraged to solve this task according to the task description, using any language you may know.
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
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
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
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
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
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
#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#
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
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
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
: 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
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
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
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
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
Ring
load "stdlib.ring"
see "working..." + nl
for i = 1 to 5000
if isWeiferich(i)
see "" + i + nl
ok
next
see "done..." + nl
function isWeiferich(p)
if not isPrime(p)
return False
ok
q = 1
p2 = pow(p,2)
while p > 1
q = (2 * q) % p2
p -= 1
end
if q = 1
return True
else
return False
ok
- Output:
working... 1093 3511 done...
RPL
« { } 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
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
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
Categories:
- Programming Tasks
- Prime Numbers
- WikipediaSourced
- Ada
- ALGOL 68
- ALGOL 68-primes
- APL
- Arturo
- AWK
- BASIC
- BASIC256
- PureBasic
- Run BASIC
- Yabasic
- C
- C++
- C sharp
- EasyLang
- F Sharp
- Factor
- Fermat
- Forth
- FreeBASIC
- Go
- Go-rcu
- Haskell
- J
- Java
- Jq
- Julia
- Mathematica
- Wolfram Language
- Nim
- Bignum
- PARI/GP
- Perl
- Ntheory
- Phix
- PicoLisp
- Python
- Quackery
- Racket
- Raku
- REXX
- Ring
- RPL
- Ruby
- Rust
- Sidef
- Swift
- Wren
- Wren-math
- Wren-big