Jump to content

Super-Poulet numbers

From Rosetta Code
(Redirected from Super-Poulet number)
Task
Super-Poulet numbers
You are encouraged to solve this task according to the task description, using any language you may know.

A super-Poulet number is a Poulet number (or Fermat pseudoprime to base 2) whose every divisor d evenly divides 2d − 2.


Task
  • Find and display the first 20 super-Poulet numbers.


Stretch
  • Find and display the index and value of the first super-Poulet number greater than one million.


See also

ALGOL 68

Uses a sieve for the primes and caches the "2 ^ divisor mod divisor" values.
The powmod routine is based on a translation of the EasyLang code.
Note that if running this with ALGOL 68 Genie version 3, you will probably need to specify a large heap size

Note, the source of primes.incl.a68 (the RC ALGOL 68-primes library) is on a separate page on Rosetta Code - see the above link.

BEGIN # find some Super-Poulet numbers                                       #
    PR read "primes.incl.a68" PR                   # include prime utilities #
    []BOOL is prime = PRIMESIEVE 1 200 000; # sieve primes up to 1.2 million #

    # returns a in ^ b in MOD m                                              #
    PROC powmod = ( INT a in, b in, m )INT:
         BEGIN
            LONG INT r := 1, a := a in, b := b in;
            WHILE b >= 1 DO
                IF b MOD 2 = 1 THEN r *:= a MODAB m FI;
                b OVERAB 2;
                a *:= a MODAB m
            OD;
            SHORTEN r
         END # powmod # ;

    [ 1 : UPB is prime ]INT powmod2dd;         # cached powmod(2,d,d) values #
    FOR d TO UPB powmod2dd DO powmod2dd[ d ] := -1 OD;

    # returns powmod( 2, d, d ) using/storing the cache                      #
    PROC cached powmod 2dd = ( INT d )INT:
         IF   powmod2dd[ d ] >= 0
         THEN powmod2dd[ d ]
         ELSE powmod2dd[ d ] := powmod( 2, d, d )
         FI # cached powmod 2dd # ;

    # returns TRUE if n is a Super Poulet number, FALSE otherwise            #
    PROC is super poulet = ( INT n )BOOL:
         IF   is prime[ n ]
         THEN FALSE
         ELIF NOT ODD n
         THEN FALSE
         ELIF powmod( 2, n - 1, n ) /= 1
         THEN FALSE
         ELSE BOOL result := TRUE;
              FOR d FROM 3 BY 2 WHILE d * d <= n AND result DO
                  IF n MOD d = 0 THEN
                      IF result := cached powmod 2dd( d ) = 2 THEN
                          INT q = n OVER d;
                          IF q /= d THEN
                             result := cached powmod 2dd( q ) = 2
                          FI
                      FI
                  FI
              OD;
              result
         FI # is super poulet # ;

    INT  count     := 0;
    BOOL searching := TRUE;
    FOR pn FROM 3 BY 2 WHILE searching DO
        IF is super poulet( pn ) THEN
            IF   ( count +:= 1 ) <= 20 THEN
                print( ( " ", whole( pn, -5 ) ) );
                IF count MOD 10 = 0 THEN print( ( newline ) ) FI
            ELIF pn > 1 000 000 THEN
                print( ( newline, "Super-Poulet number ", whole( count, 0 ) ) );
                print( ( " (", whole( pn, 0 ), ") is the first over 1 000 000", newline ) )
                searching := FALSE
            FI
        FI
    OD
END
Output:
   341  1387  2047  2701  3277  4033  4369  4681  5461  7957
  8321 10261 13747 14491 15709 18721 19951 23377 31417 31609

Super-Poulet number 109 (1016801) is the first over 1 000 000

BASIC

FreeBASIC

Translation of: Nim
'#include "isprime.bas"

Function powMod (a As Ulongint, n As Ulongint, m As Ulongint) As Ulongint
    Dim As Ulongint vase = a Mod m
    Dim As Ulongint expo = n
    Dim As Ulongint result = 1
    
    While expo > 0
        If (expo And 1) Then result = (result * vase) Mod m
        expo Shr= 1
        vase = (vase * vase) Mod m
    Wend
    Return result
End Function

Sub getDivisors(n As Ulongint, divs() As Ulongint, Byref cnt As Integer)
    cnt = 0
    Redim divs(0)
    divs(cnt) = n
    cnt += 1
    
    Dim As Ulongint d, q
    d = 2
    While d * d <= n
        If n Mod d = 0 Then
            Redim Preserve divs(cnt)
            divs(cnt) = d
            cnt += 1
            q = n \ d
            If q <> d Then
                Redim Preserve divs(cnt)
                divs(cnt) = q
                cnt += 1
            End If
        End If
        d += 1
    Wend
End Sub

Function isSuperPouletNumber(n As Ulongint) As Boolean
    If isPrime(n) Or powMod(2, n - 1, n) <> 1 Then Return False
    
    Dim As Ulongint divs(100)
    Dim As Integer i, cnt
    getDivisors(n, divs(), cnt)
    
    For i = 0 To cnt - 1
        If powMod(2, divs(i), divs(i)) <> 2 Then Return False
    Next
    Return True
End Function

' Main program
Dim As Ulongint n = 2, cnt = 0
Dim As Boolean found1M = False, found10M = False

Print "First 20 super-Poulet numbers:"
While Not (found1M And found10M)
    If isSuperPouletNumber(n) Then
        cnt += 1
        If cnt <= 20 Then
            Print Using "######"; n;
            If cnt Mod 5 = 0 Then Print
        Elseif n > 1000000 And Not found1M Then
            Print Using !"\nFirst super-Poulet number greater than one million is & at index &"; n; cnt
            found1M = True
        Elseif n > 10000000 And Not found10M Then
            Print Using !"\nFirst super-Poulet number greater than ten million is & at index &"; n; cnt
            found10M = True
        End If
    End If
    n += 1
Wend

Sleep
Output:
First 20 super-Poulet numbers:
   341  1387  2047  2701  3277
  4033  4369  4681  5461  7957
  8321 10261 13747 14491 15709
 18721 19951 23377 31417 31609

First super-Poulet number greater than one million is 1016801 at index 109

First super-Poulet number greater than ten million is 10031653 at index 317

PureBasic

Translation of: FreeBASIC
XIncludeFile "isprime.pb"

Procedure.q powMod(a.q, n.q, m.q)
  bbase.q = a % m
  expo.q = n
  result.q = 1
  
  While expo > 0
    If expo & 1
      result = (result * bbase) % m
    EndIf
    expo >> 1
    bbase = (bbase * bbase) % m
  Wend
  ProcedureReturn result
EndProcedure

Structure Counter
  i.i
EndStructure

Procedure getDivisors(n.q, Array divs.q(1), *cnt.Counter)
  *cnt\i = 0
  ReDim divs(0)
  divs(0) = n
  *cnt\i + 1
  
  d.q = 2
  While d * d <= n
    If n % d = 0
      ReDim divs(*cnt\i)
      divs(*cnt\i) = d
      *cnt\i + 1
      q.q = n / d
      If q <> d
        ReDim divs(*cnt\i)
        divs(*cnt\i) = q
        *cnt\i + 1
      EndIf
    EndIf
    d + 1
  Wend
EndProcedure

Procedure.i isSuperPouletNumber(n.q)
  If isPrime(n) Or powMod(2, n - 1, n) <> 1
    ProcedureReturn #False
  EndIf
  
  Dim divs.q(100)
  cnt.i = 0
  getDivisors(n, divs(), @cnt)
  
  For i = 0 To cnt - 1
    If powMod(2, divs(i), divs(i)) <> 2
      ProcedureReturn #False
    EndIf
  Next
  ProcedureReturn #True
EndProcedure

OpenConsole()
n.q = 2
cnt.q = 0
found1M.i = #False
found10M.i = #False

PrintN("First 20 super-Poulet numbers:")
While Not (found1M And found10M)
  If isSuperPouletNumber(n)
    cnt + 1
    If cnt <= 20
      Print(RSet(Str(n), 6))
      If cnt % 5 = 0 : PrintN("") : EndIf
    ElseIf n > 1000000 And Not found1M
      PrintN(#CRLF$ + "First super-Poulet number greater than one million is " + Str(n) + " at index " + Str(cnt))
      found1M = #True
    ElseIf n > 10000000 And Not found10M
      PrintN(#CRLF$ + "First super-Poulet number greater than ten million is " + Str(n) + " at index " + Str(cnt))
      found10M = #True
    EndIf
  EndIf
  n + 1
Wend

PrintN(#CRLF$ + "Press ENTER to exit"): Input()
CloseConsole()
Output:
Same as FreeBASIC entry.

QB64

Translation of: FreeBASIC
Dim n As _Unsigned _Integer64: n = 2
Dim cnt As _Unsigned _Integer64: cnt = 0
Dim found1M As Integer: found1M = 0
Dim found10M As Integer: found10M = 0

Print "First 20 super-Poulet numbers:"

While Not (found1M And found10M)
    If isSuperPouletNumber(n) Then
        cnt = cnt + 1
        If cnt <= 20 Then
            Print Using "######"; n;
            If cnt Mod 5 = 0 Then Print
        ElseIf n > 1000000 And Not found1M Then
            Print
            Print "First super-Poulet number greater than one million is"; n; "at index"; cnt
            found1M = -1
        ElseIf n > 10000000 And Not found10M Then
            Print "First super-Poulet number greater than ten million is"; n; "at index"; cnt
            found10M = -1
        End If
    End If
    n = n + 1
Wend
End

Function isPrime% (num As _Unsigned _Integer64)
    If num < 2 Then Exit Function
    If num Mod 2 = 0 Then isPrime = 0
    If num Mod 3 = 0 Then isPrime = -1

    Dim d As _Unsigned _Integer64: d = 5
    While d * d <= num
        If num Mod d = 0 Then Exit Function
        d = d + 2
    Wend
    isPrime = -1
End Function

Function powMod (a As _Unsigned _Integer64, n As _Unsigned _Integer64, m As _Unsigned _Integer64)
    Dim bbase As _Unsigned _Integer64: bbase = a Mod m
    Dim expo As _Unsigned _Integer64: expo = n
    Dim result As _Unsigned _Integer64: result = 1

    While expo > 0
        If (expo And 1) Then result = (result * bbase) Mod m
        expo = _SHR(expo, 1)
        bbase = (bbase * bbase) Mod m
    Wend
    powMod = result
End Function

Sub getDivisors (n As _Unsigned _Integer64, divs() As _Unsigned _Integer64, cnt As Integer)
    cnt = 0
    ReDim divs(0)
    divs(cnt) = n
    cnt = cnt + 1

    Dim d As _Unsigned _Integer64: d = 2
    Dim q As _Unsigned _Integer64

    While d * d <= n
        If n Mod d = 0 Then
            ReDim _Preserve divs(cnt)
            divs(cnt) = d
            cnt = cnt + 1
            q = n \ d
            If q <> d Then
                ReDim _Preserve divs(cnt)
                divs(cnt) = q
                cnt = cnt + 1
            End If
        End If
        d = d + 1
    Wend
End Sub

Function isSuperPouletNumber% (n As _Unsigned _Integer64)
    If isPrime(n) Or powMod(2, n - 1, n) <> 1 Then Exit Function

    Dim divs(100) As _Unsigned _Integer64
    Dim i As Integer, cnt As Integer

    Call getDivisors(n, divs(), cnt)

    For i = 0 To cnt - 1
        If powMod(2, divs(i), divs(i)) <> 2 Then Exit Function
    Next
    isSuperPouletNumber = -1
End Function
Output:
Same as FreeBASIC entry.

C++

#include <algorithm>
#include <iostream>
#include <vector>

std::vector<unsigned int> divisors(unsigned int n) {
    std::vector<unsigned int> result{1};
    unsigned int power = 2;
    for (; (n & 1) == 0; power <<= 1, n >>= 1)
        result.push_back(power);
    for (unsigned int p = 3; p * p <= n; p += 2) {
        size_t size = result.size();
        for (power = p; n % p == 0; power *= p, n /= p)
            for (size_t i = 0; i != size; ++i)
                result.push_back(power * result[i]);
    }
    if (n > 1) {
        size_t size = result.size();
        for (size_t i = 0; i != size; ++i)
            result.push_back(n * result[i]);
    }
    return result;
}

unsigned long long modpow(unsigned long long base, unsigned int exp,
                          unsigned int mod) {
    if (mod == 1)
        return 0;
    unsigned long long result = 1;
    base %= mod;
    for (; exp != 0; exp >>= 1) {
        if ((exp & 1) == 1)
            result = (result * base) % mod;
        base = (base * base) % mod;
    }
    return result;
}

bool is_prime(unsigned int n) {
    if (n < 2)
        return false;
    if (n % 2 == 0)
        return n == 2;
    if (n % 3 == 0)
        return n == 3;
    if (n % 5 == 0)
        return n == 5;
    static constexpr unsigned int wheel[] = {4, 2, 4, 2, 4, 6, 2, 6};
    for (unsigned int p = 7;;) {
        for (unsigned int w : wheel) {
            if (p * p > n)
                return true;
            if (n % p == 0)
                return false;
            p += w;
        }
    }
}

bool is_poulet_number(unsigned int n) {
    return modpow(2, n - 1, n) == 1 && !is_prime(n);
}

bool is_sp_num(unsigned int n) {
    if (!is_poulet_number(n))
        return false;
    auto div = divisors(n);
    return all_of(div.begin() + 1, div.end(),
                  [](unsigned int d) { return modpow(2, d, d) == 2; });
}

int main() {
    std::cout.imbue(std::locale(""));
    unsigned int n = 1, count = 0;
    std::cout << "First 20 super-Poulet numbers:\n";
    for (; count < 20; n += 2) {
        if (is_sp_num(n)) {
            ++count;
            std::cout << n << ' ';
        }
    }
    std::cout << '\n';
    for (unsigned int limit = 1000000; limit <= 10000000; limit *= 10) {
        for (;;) {
            n += 2;
            if (is_sp_num(n)) {
                ++count;
                if (n > limit)
                    break;
            }
        }
        std::cout << "\nIndex and value of first super-Poulet greater than "
                  << limit << ":\n#" << count << " is " << n << '\n';
    }
}
Output:
First 20 super-Poulet numbers:
341 1,387 2,047 2,701 3,277 4,033 4,369 4,681 5,461 7,957 8,321 10,261 13,747 14,491 15,709 18,721 19,951 23,377 31,417 31,609 

Index and value of first super-Poulet greater than 1,000,000:
#109 is 1,016,801

Index and value of first super-Poulet greater than 10,000,000:
#317 is 10,031,653

EasyLang

Translation of: Nim
fastfunc powmod a b m .
   r = 1
   while b >= 1
      if b mod 2 = 1
         r = r * a mod m
      .
      b = b div 2
      a = a * a mod m
   .
   return r
.
fastfunc isprim num .
   i = 2
   while i <= sqrt num
      if num mod i = 0
         return 0
      .
      i += 1
   .
   return 1
.
func[] divisors n .
   d = 2
   while d * d <= n
      if n mod d = 0
         r[] &= d
         q = n div d
         if q <> d
            r[] &= q
         .
      .
      d += 1
   .
   return r[]
.
func is_super_poulet n .
   if isprim n = 1 or powmod 2 (n - 1) n <> 1
      return 0
   .
   for d in divisors n
      if powmod 2 d d <> 2
         return 0
      .
   .
   return 1
.
n = 3
while 1 = 1
   if is_super_poulet n = 1
      cnt += 1
      if cnt <= 20
         write n & " "
      .
      if n > 1e6
         break 1
      .
   .
   n += 2
.
print ""
print ""
print cnt & " " & n
Output:
341 1387 2047 2701 3277 4033 4369 4681 5461 7957 8321 10261 13747 14491 15709 18721 19951 23377 31417 31609 

109 1016801

Factor

Works with: Factor version 0.99 2022-04-03
USING: combinators.short-circuit io kernel lists lists.lazy math
math.functions math.primes math.primes.factors prettyprint
sequences ;

: super-poulet? ( n -- ? )
    {
        [ prime? not ]
        [ [ 1 - 2^ ] keep mod 1 = ]
        [ divisors [ [ 2^ 2 - ] keep divisor? ] all? ]
    } 1&& ;

: super-poulets ( -- list )
    1 lfrom [ super-poulet? ] lfilter ;

20 super-poulets ltake [ pprint bl ] leach nl
Output:
341 1387 2047 2701 3277 4033 4369 4681 5461 7957 8321 10261 13747 14491 15709 18721 19951 23377 31417 31609 

J

Implementation (only good for the first 60 super-poulet numbers):

spou1=: {{ 2 = 2x(y&|)@^ y }}

is_super_poulet=: {{
   if. 2~:#q=. q: y do. 0 return. end.
   if. spou1 {. q do.
     if. spou1 {: q do.
       if. spou1 y do. 1 return. end.
     end.
   end.
   0
}}"0

Task example:

   20{. (#~ is_super_poulet) 1+i.1e5
341 1387 2047 2701 3277 4033 4369 4681 5461 7957 8321 10261 13747 14491 15709 18721 19951 23377 31417 31609

Java

import java.util.Set;
import java.util.TreeSet;

public final class SuperPouletNumbers {

	public static void main(String[] args) {
		int n = 2;
		int count = 0;
		boolean searching = true;
		while ( searching ) {
		    if ( isSuperPouletNumber(n) ) {
		    	count += 1;
		    	if ( count <= 20 ) {
		    		System.out.print(String.format("%6d%s", n, ( count % 5 == 0 ? "\n" : " " )));
		    	} else if ( n > 1_000_000 ) {
		    		System.out.println(System.lineSeparator() +
		    			"The first super-Poulet number greater than one million is " + n + " at index " + count);
		    		searching = false;
		    	}
		    }
		    n += 1;
		}		
	}	
	
	private static boolean isSuperPouletNumber(int number) {
		if ( isPrime(number) || modulusPower(2, number - 1, number) != 1 ) {
			return false;
		}
		
		for ( int divisor : divisors(number) ) {
			if ( modulusPower(2, divisor, divisor) != 2 ) {
				return false;
			}
		}
		return true;
	}
	
	// Return the divisors of the given number, excluding 1
	private static Set<Integer> divisors(int number) {
		Set<Integer> result = new TreeSet<Integer>();
		for ( int d = 2; d * d <= number; d++ ) { 
	        if ( number % d == 0 ) { 
	        	result.add(d);  
	            result.add(number / d);
	        } 
	    } 
		result.add(number);
		return result;
	} 
	
	private static long modulusPower(long base, long exponent, long modulus) {
	    if ( modulus == 1 ) {
	        return 0;
	    }	    
	    
	    base %= modulus;
	    long result = 1;
	    while ( exponent > 0 ) {
	        if ( ( exponent & 1 ) == 1 ) {
	            result = ( result * base ) % modulus;
	        }
	        base = ( base * base ) % modulus;
	        exponent >>= 1;
	    }
	    return result;
	}
	
	private static boolean isPrime(int number) {
		if ( number < 2 ) {
			return false;
		}
		if ( number % 2 == 0 ) {
			return number == 2;
		}
		if ( number % 3 == 0 ) {
			return number == 3;
		}
		
		int delta = 2;
		int k = 5;
		while ( k * k <= number ) {
		    if ( number % k == 0 ) {
		    	return false;
		    }
		    k += delta;
		    delta = 6 - delta;
		}
		return true;
	}

}
Output:
   341   1387   2047   2701   3277
  4033   4369   4681   5461   7957
  8321  10261  13747  14491  15709
 18721  19951  23377  31417  31609

The first super-Poulet number greater than one million is 1016801 at index 109

Julia

using Primes

divisors(n) = @view sort!(vec(map(prod, Iterators.product((p.^(0:m) for (p, m) in eachfactor(n))...))))[begin:end-1]

issuperPoulet(n) = !isprime(n) && big"2"^(n - 1) % n == 1 && all(d -> (big"2"^d - 2) % d == 0, divisors(n))

spoulets = filter(issuperPoulet, 1:12_000_000)

println("The first 20 super-Poulet numbers are: ", spoulets[1:20])

idx1m = findfirst(>(1_000_000), spoulets)
idx10m = findfirst(>(10_000_000), spoulets)

println("The first super-Poulet number over 1 million is the ", idx1m, "th one, which is ", spoulets[idx1m])
println("The first super-Poulet number over 10 million is the ", idx10m, "th one, which is ", spoulets[idx10m])
Output:
The first 20 super-Poulet numbers are: [341, 1387, 2047, 2701, 3277, 4033, 4369, 4681, 5461, 7957, 8321, 10261, 13747, 14491, 15709, 18721, 19951, 23377, 31417, 31609]
The first super-Poulet number over 1 million is the 109th one, which is 1016801
The first super-Poulet number over 10 million is the 317th one, which is 10031653

Mathematica /Wolfram Language

Translation of: Python
(*Define the super-Poulet number test*)
IsSuperPoulet[n_Integer] := 
  Module[{divisors, condition1, condition2}, divisors = Divisors[n];
   condition1 = Not[PrimeQ[n]] && Mod[2^(n - 1), n] == 1;
   condition2 = AllTrue[Most@divisors, Mod[2^# - 2, #] == 0 &];
   condition1 && condition2];

(*Generate super-Poulet numbers up to 1,100,000*)
superPoulets = Select[Range[2, 1100000], IsSuperPoulet];

(*Print the first 20 super-Poulet numbers*)
Print["The first 20 super-Poulet numbers are: ", 
  superPoulets[[1 ;; 20]]];

(*Find and print the first super-Poulet number over 1 million*)
firstOverOneMillion = First@Select[superPoulets, # > 1000000 &];
idx1m = Position[superPoulets, firstOverOneMillion][[1, 1]];
Print["The first super-Poulet number over 1 million is the ", idx1m, 
  "th one, which is ", firstOverOneMillion];
Output:
The first 20 super-Poulet numbers are: {341, 1387, 2047, 2701, 3277, 4033, 4369, 4681, 5461, 7957, 8321, 10261, 13747, 14491, 15709, 18721, 19951, 23377, 31417, 31609}
The first super-Poulet number over 1 million is the 109th one, which is 1016801


Nim

import std/[strformat, strutils]

proc powMod*(a, n, m: int): int =
  ## Return "a^n mod m".
  var a = a mod m
  var n = n
  if a > 0:
    result = 1
    while n > 0:
      if (n and 1) != 0:
        result = (result * a) mod m
      n = n shr 1
      a = (a * a) mod m

func isPrime(n: Natural): bool =
  ## Return true if "n" is prime.
  if n < 2: return false
  if (n and 1) == 0: return n == 2
  if n mod 3 == 0: return n == 3
  var k = 5
  var delta = 2
  while k * k <= n:
    if n mod k == 0: return false
    inc k, delta
    delta = 6 - delta
  result = true

iterator divisors(n: Positive): int =
  ## Yield the divisors of "n", except 1.
  yield n
  var d = 2
  while d * d <= n:
    if n mod d == 0:
      let q = n div d
      yield d
      if q != d:
        yield q
    inc d

func isSuperPouletNumber(n: int): bool =
  ## Return true is "x" is a Fermat pseudoprime to base "a".
  if n.isPrime or powMod(2, n - 1, n) != 1: return false
  for d in n.divisors:
    if powMod(2, d, d) != 2:
      return false
  result = true

var n = 2
var count = 0
while true:
  if n.isSuperPouletNumber:
    inc count
    if count <= 20:
      stdout.write &"{n:5}"
      stdout.write if count mod 5 == 0: '\n' else: ' '
      if count == 20: echo()
    elif n > 1_000_000:
      echo &"First super-Poulet number greater than one million is {insertSep($n)} at index {count}."
      break
  inc n
Output:
  341  1387  2047  2701  3277
 4033  4369  4681  5461  7957
 8321 10261 13747 14491 15709
18721 19951 23377 31417 31609

First super-Poulet number greater than one million is 1_016_801 at index 109.

Perl

Library: ntheory
use v5.36;
use experimental <builtin for_list>;
use List::AllUtils <firstidx all>;
use ntheory <is_prime divisors powmod>;

sub comma { reverse ((reverse shift) =~ s/(.{3})/$1,/gr) =~ s/^,//r }

my @poulet       = grep { !is_prime($_) && (1 == powmod 2, $_ - 1, $_) } 2 .. 1.1e7;
my @super_poulet = grep { all { 2 == powmod 2, $_, $_ } grep { $_ > 1 } divisors $_ } @poulet;

say "First 20 super-Poulet numbers:\n" . join ' ', @super_poulet[0..19];
for my($i,$j) (1, 1e6, 10, 1e7) {
    my $index = firstidx { $_ > $j } @super_poulet;
    say "\nIndex and value of first super-Poulet greater than $i million:";
    say "#@{[1+$index]} is " . comma $super_poulet[$index];
}
Output:
First 20 super-Poulet numbers:
341 1387 2047 2701 3277 4033 4369 4681 5461 7957 8321 10261 13747 14491 15709 18721 19951 23377 31417 31609

Index and value of first super-Poulet greater than 1 million:
#109 is 1,016,801

Index and value of first super-Poulet greater than 10 million:
#317 is 10,031,653

Phix

with javascript_semantics
include mpfr.e 
constant two = mpz_init(2)

function super_poulet(integer x)
    if is_prime(x) or x<=1 then return false end if
    mpz z = mpz_init(x)
    mpz_powm_ui(z, two, x-1, z)
    if mpz_cmp_si(z,1)!=0 then return false end if
    sequence f = factors(x)
    for i=1 to length(f) do
        mpz_set_si(z,f[i])
        mpz_powm(z, two, z, z)
        if mpz_cmp_si(z,2)!=0 then return false end if
    end for
    return true
end function
 
integer count = 0, x = 1
sequence first20 = {}
while count<20 do
    if super_poulet(x) then
        if count<20 then first20 &= x end if
        count += 1
    end if
    x += 2
end while
printf(1,"First 20 super-Poulet numbers:\n%v\n\n",{first20})
for lim in {1e6,1e7} do
    while true do
        x += 2
        if super_poulet(x) then
            count += 1
            if x>lim then exit end if
        end if
    end while
    printf(1,"The %d%s super-Poulet number is %,d and the first greater than %,d\n",
             {count,ord(count),x,lim})
end for
Output:
First 20 super-Poulet numbers:
{341,1387,2047,2701,3277,4033,4369,4681,5461,7957,8321,10261,13747,14491,15709,18721,19951,23377,31417,31609}

The 109th super-Poulet number is 1,016,801 and the first greater than 1,000,000
The 317th super-Poulet number is 10,031,653 and the first greater than 10,000,000

Python

Translation of: Julia
from sympy import isprime, divisors
 
def is_super_Poulet(n):
    return not isprime(n) and 2**(n - 1) % n == 1 and all((2**d - 2) % d == 0 for d in divisors(n))

spoulets = [n for n in range(1, 1_100_000) if is_super_Poulet(n)]

print('The first 20 super-Poulet numbers are:', spoulets[:20])

idx1m, val1m = next((i, v) for i, v in enumerate(spoulets) if v > 1_000_000)
print(f'The first super-Poulet number over 1 million is the {idx1m}th one, which is {val1m}')
Output:
The first 20 super-Poulet numbers are: [341, 1387, 2047, 2701, 3277, 4033, 4369, 4681, 5461, 7957, 8321, 10261, 13747, 14491, 15709, 18721, 19951, 23377, 31417, 31609]
The first super-Poulet number over 1 million is the 108th one, which is 1016801


Raku

use Prime::Factor;
use Lingua::EN::Numbers;

my @poulet = lazy (2..*).hyper(:2000batch).grep: { !.is-prime && (1 == expmod 2, $_ - 1, $_) }
my @super-poulet = @poulet.grep: { all .&divisors.skip(1).map: { 2 == expmod 2, $_, $_ } }

say "First 20 super-Poulet numbers:\n" ~ @super-poulet[^20].gist;

for 1e6.Int, 1e7.Int -> $threshold {
    say "\nIndex and value of first super-Poulet greater than {$threshold.&cardinal}:";
    my $index = @super-poulet.first: * > $threshold, :k;
    say "{(1+$index).&ordinal-digit} super-Poulet number == " ~ @super-poulet[$index].&comma;
}
Output:
First 20 super-Poulet numbers:
(341 1387 2047 2701 3277 4033 4369 4681 5461 7957 8321 10261 13747 14491 15709 18721 19951 23377 31417 31609)

Index and value of first super-Poulet greater than one million:
109th super-Poulet number == 1,016,801

Index and value of first super-Poulet greater than ten million:
317th super-Poulet number == 10,031,653

Ruby

require 'prime'

class Integer
  def proper_divisors
    return [] if self == 1
    primes = prime_division.flat_map{|prime, freq| [prime] * freq}
    (1...primes.size).each_with_object([1]) do |n, res|
      primes.combination(n).map{|combi| res << combi.inject(:*)}
    end.flatten.uniq.sort
  end
end

super_poulets = (1..).lazy.select do |n|
  n.prime? == false  && 
  2.pow(n-1, n) == 1 && # modular exponentiation
  n.proper_divisors[1..].all?{|d| 2.pow(d, d) == 2} # again
end 

m = 20
puts "First #{m} super-Poulet numbers:\n#{super_poulets.first(m).join(-", ") }"

[1_000_000, 10_000_000].each do |m|
  puts "\nValue and index of first super-Poulet number greater than #{m}:"
  puts "%d is #%d" % super_poulets.with_index(1).detect{|n, i| n > m }
end
Output:
First 20 super-Poulet numbers:
341, 1387, 2047, 2701, 3277, 4033, 4369, 4681, 5461, 7957, 8321, 10261, 13747, 14491, 15709, 18721, 19951, 23377, 31417, 31609

Value and index of first super-Poulet number greater than 1000000:
1016801 is #109

Value and index of first super-Poulet number greater than 10000000:
10031653 is #317

Swift

Translation of: C++
func divisors(number: UInt32) -> [UInt32] {
    var result: [UInt32] = [1]
    var power: UInt32 = 2
    var n = number
    while (n & 1) == 0 {
        result.append(power)
        power <<= 1
        n >>= 1
    }
    var p: UInt32 = 3
    while p * p <= n {
        let size = result.count
        power = p
        while n % p == 0 {
            for i in 0..<size {
                result.append(power * result[i])
            }
            n /= p
            power *= p
        }
        p += 2
    }
    if n > 1 {
        let size = result.count
        for i in 0..<size {
            result.append(n * result[i])
        }
    }
    return result
}

func modPow(base: UInt64, exponent: UInt32, mod: UInt32) -> UInt64 {
    if mod == 1 {
        return 0
    }
    var result: UInt64 = 1
    var exp = exponent
    var b = base
    let m = UInt64(mod)
    b %= m
    while exp > 0 {
        if (exp & 1) == 1 {
            result = (result * b) % m
        }
        b = (b * b) % m
        exp >>= 1
    }
    return result
}

func isPrime(number: UInt32) -> Bool {
    if number < 2 {
        return false
    }
    if number % 2 == 0 {
        return number == 2
    }
    if number % 3 == 0 {
        return number == 3
    }
    if number % 5 == 0 {
        return number == 5
    }
    var p: UInt32 = 7
    let wheel: [UInt32] = [4,2,4,2,4,6,2,6]
    while true {
        for w in wheel {
            if p * p > number {
                return true
            }
            if number % p == 0 {
                return false
            }
            p += w
        }
    }
}

func isPouletNumber(_ n: UInt32) -> Bool {
    return modPow(base: 2, exponent: n - 1, mod: n) == 1 && !isPrime(number: n)
}

func isSuperPouletNumber(_ n: UInt32) -> Bool {
    if (!isPouletNumber(n)) {
        return false
    }
    let div = divisors(number: n)
    return div[1...].allSatisfy({modPow(base: 2, exponent: $0, mod: $0) == 2})
}

var n: UInt32 = 1
var count = 0

print("First 20 super-Poulet numbers:")
while count < 20 {
    if (isSuperPouletNumber(n)) {
        count += 1
        print("\(n)", terminator: " ")
    }
    n += 2
}
print()

var limit = 1000000
while limit <= 10000000 {
    while true {
        n += 2
        if (isSuperPouletNumber(n)) {
            count += 1
            if (n > limit) {
                break
            }
        }
    }
    print("\nIndex and value of first super-Poulet greater than \(limit):\n#\(count) is \(n)")
    limit *= 10
}
Output:
First 20 super-Poulet numbers:
341 1387 2047 2701 3277 4033 4369 4681 5461 7957 8321 10261 13747 14491 15709 18721 19951 23377 31417 31609 

Index and value of first super-Poulet greater than 1000000:
#109 is 1016801

Index and value of first super-Poulet greater than 10000000:
#317 is 10031653

Wren

Library: Wren-math
Library: Wren-gmp
Library: Wren-fmt
import "./math" for Int
import "./gmp" for Mpz
import "./fmt" for Fmt

var one = Mpz.one

var isSuperPoulet = Fn.new { |x|
    if (Int.isPrime(x)) return false
    var bx = Mpz.from(x)
    if (Mpz.two.modPow(x-1, bx) != one) return false
    var t = Mpz.new()
    return Int.divisors(x).skip(1).all { |d| t.uiPow(2, d).sub(2).isDivisibleUi(d) }
}

var count = 0
var first20 = List.filled(20, 0)
var x = 3
while (count < 20) {
    if (isSuperPoulet.call(x)) {
        first20[count] = x
        count = count + 1
    }
    x = x + 2  // Poulet numbers are always odd
}
System.print("The first 20 super-Poulet numbers are:")
System.print(first20)
System.print()
var limit = 1e6
while (true) {
   if (isSuperPoulet.call(x)) {
        count = count + 1
        if (x > limit) {
            Fmt.print("The $r super-Poulet number and the first over $,d is $,d.", count, limit, x)
            if (limit == 1e6) limit = 1e7 else return
        }
   }
   x = x + 2
}
Output:
The first 20 super-Poulet numbers are:
[341, 1387, 2047, 2701, 3277, 4033, 4369, 4681, 5461, 7957, 8321, 10261, 13747, 14491, 15709, 18721, 19951, 23377, 31417, 31609]

The 109th super-Poulet number and the first over 1,000,000 is 1,016,801.
The 317th super-Poulet number and the first over 10,000,000 is 10,031,653.
Cookies help us deliver our services. By using our services, you agree to our use of cookies.