Wieferich primes: Difference between revisions
m →{{header|REXX}}: simplified program, added/changed whitespace and comments. |
m →{{header|REXX}}: changed some comments and whitespace. |
||
Line 558: | Line 558: | ||
!.=0; !.2=1; !.3=1; !.5=1; !.7=1; !.11=1 /* " " " " (semaphores).*/ |
!.=0; !.2=1; !.3=1; !.5=1; !.7=1; !.11=1 /* " " " " (semaphores).*/ |
||
#= 5; sq.#= @.# **2 /*number of primes so far; prime². */ |
#= 5; sq.#= @.# **2 /*number of primes so far; prime². */ |
||
do j=@.#+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? */ |
|||
if j// 3==0 then iterate /*" " " 3? */ |
|||
if j// 7==0 then iterate /*" " " 7? */ |
|||
do k=5 while sq.k<=j /* [↓] divide by the known odd primes.*/ |
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. ___ */ |
if j//@.k==0 then iterate j /*Is J ÷ a P? Then not prime. ___ */ |
Revision as of 04:01, 8 July 2021
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
<lang 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;
}</lang>
- Output:
Wieferich primes less than 5000: 1093 3511
C++
<lang cpp>#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';
}</lang>
- Output:
Wieferich primes less than 5000: 1093 3511
C#
<lang csharp>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); } } }
}</lang>
- Output:
Wieferich primes less that 5000: 1093 3511
F#
This task uses Extensible Prime Generator (F#) <lang fsharp> // 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") </lang>
- Output:
1093 3511 Real: 00:00:00.004
Factor
<lang 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 .</lang>
- Output:
Wieferich primes less than 5000: V{ 1093 3511 }
Forth
<lang 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</lang>
- Output:
Wieferich primes less than 5000: 1093 3511
Go
<lang 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)) } }
}</lang>
- Output:
Wieferich primes < 5,000: 1,093 3,511
Java
<lang 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; }
}</lang>
- Output:
Wieferich primes less than 5000: 1093 3511
Julia
<lang julia>using Primes
println(filter(p -> (big"2"^(p - 1) - 1) % p^2 == 0, primes(5000))) # [1093, 3511] </lang>
Nim
<lang 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</lang>
- Output:
Wieferich primes less than 5000: 1093 3511
Perl
<lang 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;</lang>
- 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)})
Quackery
eratosthenes
and isprime
are defined at Sieve of Eratosthenes#Quackery.
<lang 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 ] ]</lang>
- Output:
1093 3511
Raku
<lang perl6>put "Wieferich primes less than 5000: ", join ', ', ^5000 .grep: { .is-prime and not ( exp($_-1, 2) - 1 ) % .² };</lang>
- Output:
Wieferich primes less than 5000: 1093, 3511
REXX
<lang 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?*/ 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, '─') /*display a foot sep for the output. */ say say 'Found ' commas(found) title /*display a summary for the output. */ 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 j=@.#+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</lang>
- 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
<lang 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); }
}</lang>
- Output:
Wieferich primes less than 5000: 1093 3511
Sidef
<lang ruby>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))</lang>
- Output:
Wieferich primes less than 5000: [1093, 3511]
Swift
<lang 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)
}</lang>
- Output:
Wieferich primes less than 5000: 1093 3511
Wren
<lang ecmascript>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)
}</lang>
- Output:
Wieferich primes < 5000: 1093 3511