1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980 |
- #!/usr/bin/perl
- # Generate Super-Poulet pseudoprimes to base 2, using prime factors bellow a certain limit.
- use 5.020;
- use warnings;
- use experimental qw(signatures);
- use IO::Handle;
- use ntheory qw(forcomb forprimes divisors powmod znorder);
- use Math::Prime::Util::GMP;
- sub super_poulet_pseudoprimes ($limit, $callback) {
- my %common_divisors;
- warn ":: Sieving...\n";
- forprimes {
- my $p = $_;
- #if ($p % 24 == 1 or $p % 24 == 13) {
- #if ($p % 24 == 1) {
- my $z = znorder(2, $p);
- if (2*$z < $limit) {
- foreach my $d (divisors($p - 1)) {
- #if (powmod(2, $d, $p) == 1) {
- if ($d % $z == 0) {
- if ($p > 1e8) {
- if (exists $common_divisors{$d}) {
- push @{$common_divisors{$d}}, $p;
- }
- }
- else {
- push @{$common_divisors{$d}}, $p;
- }
- }
- #}
- }
- }
- } 3, $limit;
- warn ":: Creating combinations...\n";
- #foreach my $arr (values %common_divisors) {
- while (my ($key, $arr) = each %common_divisors) {
- my $nf = 5; # minimum number of prime factors
- #next if @$arr < $nf;
- my $l = scalar(@$arr);
- foreach my $k ($nf .. $l) {
- forcomb {
- my $n = Math::Prime::Util::GMP::vecprod(@{$arr}[@_]);
- $callback->($n);
- } $l, $k;
- }
- }
- }
- open my $fh, '>', 'super_poulet_numbers.txt';
- #$fh->autoflush(1);
- super_poulet_pseudoprimes(
- 2e9, # limit of the largest prime factor
- sub ($n) {
- if ($n > ~0) { # report only numbers greater than 2^64
- #warn "$n\n";
- say $fh $n;
- }
- }
- );
- close $fh;
|