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
Wieferich primes is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.
This page uses content from Wikipedia. The original article was at Wieferich prime. The list of authors can be seen in the page history. As with Rosetta Code, the text of Wikipedia is available under the GNU FDL. (See links for details on variance)


In number theory, a Wieferich prime is a prime number p such that p2 evenly divides 2(p − 1) − 1 .


It is conjectured that there are infinitely many Wieferich primes, but as of March 2021,only two have been identified.


Task
  • Write a routine (function procedure, whatever) to find Wieferich primes.
  • Use that routine to identify and display all of the Wieferich primes less than 5000.


See also


C[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 }

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

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

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

Julia[edit]

using Primes
 
println(filter(p -> (big"2"^(p - 1) - 1) % p^2 == 0, primes(5000))) # [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

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

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


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?*/
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). */
 !.=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 _ /*find odd primes from here on.*/
if _==5 then iterate /*get right digit; J ÷ by 5? */
if j//3==0 then iterate; if j//7==0 then iterate /*J ÷ by 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

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