Wieferich primes

From Rosetta Code
Revision as of 23:27, 1 June 2021 by Thundergnat (talk | contribs) (Thundergnat moved page Weiferich primes to Wieferich primes: Misspelled name)
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++

<lang cpp>#include <cstdint>

  1. include <iostream>
  2. 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

Factor

Works with: Factor version 0.99 2021-02-05

<lang factor>USING: io kernel math math.functions math.primes prettyprint sequences ;

"Weiferich primes less than 5000:" print 5000 primes-upto [ [ 1 - 2^ 1 - ] [ sq divisor? ] bi ] filter .</lang>

Output:
Weiferich primes less than 5000:
V{ 1093 3511 }

Forth

Works with: Gforth

<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

Julia

<lang julia>using Primes

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

Perl

Library: ntheory

<lang perl>use feature 'say'; use bignum; use ntheory 'is_prime';

say 'Weiferich primes less than 5000: ' . join ', ', grep { is_prime($_) and not ( (2**($_-1) -1) % $_**2 ) } 1..5000;</lang>

Output:
Weiferich primes less than 5000: 1093, 3511

Phix

include mpfr.e
mpz z = mpz_init()
function weiferich(integer p)
    mpz_set_str(z,repeat('1',p-1),2)
    return mpz_fdiv_q_ui(z,z,p*p)=0
end function
printf(1,"Weiferich primes less than 5000: %V\n",{filter(get_primes_le(5000),weiferich)})
Output:
Weiferich primes less than 5000: {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 Weiferich 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= ' Weiferich 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 for the output. */ wp= 0 /*initialize number of Weiferich primes*/

     do j=1  to #;    p= @.j;     pm= p - 1     /*search for Weiferich primes in range.*/
     if (2**pm - 1) // p**2\==0  then iterate   /*P**2 not evenly divide  2**(P-1) - 1?*/
     wp= wp + 1                                 /*bump the counter of Weiferich primes.*/
     say center(wp, 7)'│'  center(commas(p), w) /*display the Weiferich prime to term. */
     end   /*j*/

say '───────┴'center("" , w, '─') /*display a foot sep for the output. */ say say 'Found ' commas(wp) title 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: !.= 0 /*placeholders for primes (semaphores).*/

     @.1=2;  @.2=3;  @.3=5;  @.4=7;  @.5=11     /*define some low primes.              */
     !.2=1;  !.3=1;  !.5=1;  !.7=1;  !.11=1     /*   "     "   "    "     flags.       */
                       #=5;     s.#= @.# **2    /*number of primes so far;     prime². */
                                                /* [↓]  generate more  primes  ≤  high.*/
       do j=@.#+2  by 2  to n-1                 /*find odd primes from here on.        */
       parse var j  -1 _; if     _==5  then iterate  /*J divisible by 5?  (right dig)*/
                            if j// 3==0  then iterate  /*"     "      " 3?             */
                            if j// 7==0  then iterate  /*"     "      " 7?             */
                                                /* [↑]  the above five lines saves time*/
              do k=5  while s.k<=j              /* [↓]  divide by the known odd primes.*/
              if j // @.k == 0  then iterate j  /*Is  J ÷ X?  Then not prime.     ___  */
              end   /*k*/                       /* [↑]  only process numbers  ≤  √ J   */
       #= #+1;    @.#= j;    s.#= j*j;   !.j= 1 /*bump # of Ps; assign next P;  P²; P# */
       end          /*j*/;   return</lang>
output   when using the default input:
 index │  Weiferich primes that are  <  5,000
───────┼──────────────────────────────────────
   1   │                 1,093
   2   │                 3,511
───────┴──────────────────────────────────────

Found  2  Weiferich primes that are  <  5,000

Wren

Library: Wren-math
Library: Wren-big

<lang ecmascript>import "/math" for Int import "/big" for BigInt

var primes = Int.primeSieve(5000) System.print("Weiferich 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:
Weiferich primes < 5000:
1093
3511