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
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.