Super-Poulet numbers
Appearance
(Redirected from Super-Poulet number)
Super-Poulet numbers
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.
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
'#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
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
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
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
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
(*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
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
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].,
}
- 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
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
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.