/
generate_fermat_psp_to_multiple_bases_from_prime_factors.pl
93 lines (69 loc) · 2.28 KB
/
generate_fermat_psp_to_multiple_bases_from_prime_factors.pl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
#!/usr/bin/perl
# Author: Daniel "Trizen" Șuteu
# Date: 02 July 2022
# Edit: 12 November 2022
# https://github.com/trizen
# A new algorithm for generating Fermat pseudoprimes to multiple bases.
# See also:
# https://oeis.org/A001567 -- Fermat pseudoprimes to base 2, also called Sarrus numbers or Poulet numbers.
# https://oeis.org/A050217 -- Super-Poulet numbers: Poulet numbers whose divisors d all satisfy d|2^d-2.
# See also:
# https://en.wikipedia.org/wiki/Fermat_pseudoprime
# https://trizenx.blogspot.com/2020/08/pseudoprimes-construction-methods-and.html
use 5.020;
use warnings;
use experimental qw(signatures);
use ntheory qw(:all);
use Math::Prime::Util::GMP qw(vecprod is_pseudoprime);
use Math::GMPz;
use List::Util qw(uniq);
sub fermat_pseudoprimes ($bases, $k_limit, $callback) {
my %common_divisors;
my $bases_lcm = lcm(@$bases);
# for (my $p = 2 ; $p <= $prime_limit ; $p = next_prime($p)) {
my %seen_p;
while (<>) {
my $p = (split(' ', $_))[-1];
# foreach my $p(@plist) {
$p || next;
$p =~ /^[0-9]+\z/ or next;
is_prime($p) || next;
next if $seen_p{$p}++;
if ($p > ~0) {
$p = Math::GMPz->new($p);
}
next if ($bases_lcm % $p == 0);
my @orders = map { znorder($_, $p) } @$bases;
my $lcm_orders = lcm(@orders);
for my $k (1 .. $k_limit) {
my $q = mulint($k, $lcm_orders) + 1;
if ($p != $q and is_prime($q)) {
push @{$common_divisors{$lcm_orders}}, $q;
}
}
}
#my %seen;
foreach my $arr (values %common_divisors) {
my $l = scalar(@$arr);
@$arr = uniq(@$arr);
foreach my $k (2 .. $l) {
forcomb {
my $n = vecprod(@{$arr}[@_]);
$callback->($n) #if !$seen{$n}++;
} $l, $k;
}
}
}
my @pseudoprimes;
my @bases = @{primes(nth_prime(10))}; # generate Fermat pseudoprimes to these bases
my $k_limit = 24; # largest k multiple of the znorder(base, p)
fermat_pseudoprimes(
\@bases, # bases
$k_limit, # k limit
sub ($n) {
if ($n > ~0 and is_pseudoprime($n, @bases)) {
# push @pseudoprimes, $n;
say $n;
}
}
);