Super-Poulet numbers

Super-Poulet numbers 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.

A super-Poulet number is a Poulet number (or Fermat pseudoprime to base 2) whose every divisor d evenly divides 2^d − 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

• Task: Fermat pseudoprimes
• Wikipedia: super-Poulet number
• OEIS:A050217 - super-Poulet numbers

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

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

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

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].,
}

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

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.