123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105 |
- #!/usr/bin/perl
- # Daniel "Trizen" Șuteu
- # Date: 24 September 2022
- # https://github.com/trizen
- # Generate all the squarefree strong Fermat pseudoprimes to a given base with n prime factors in a given range [A,B]. (not in sorted order)
- # See also:
- # https://en.wikipedia.org/wiki/Almost_prime
- # https://trizenx.blogspot.com/2020/08/pseudoprimes-construction-methods-and.html
- use 5.020;
- use warnings;
- use ntheory qw(:all);
- use experimental qw(signatures);
- sub divceil ($x, $y) { # ceil(x/y)
- (($x % $y == 0) ? 0 : 1) + divint($x, $y);
- }
- sub squarefree_strong_fermat_pseudoprimes_in_range ($A, $B, $k, $base) {
- $A = vecmax($A, pn_primorial($k));
- if ($A > $B) {
- return;
- }
- my @list;
- my $generator = sub ($m, $L, $lo, $k, $k_exp, $congr) {
- my $hi = rootint(divint($B, $m), $k);
- if ($lo > $hi) {
- return;
- }
- if ($k == 1) {
- $lo = vecmax($lo, divceil($A, $m));
- $lo > $hi && return;
- my $t = invmod($m, $L);
- $t > $hi && return;
- $t += $L * divceil($lo - $t, $L) if ($t < $lo);
- for (my $p = $t ; $p <= $hi ; $p += $L) {
- is_prime($p) || next;
- $base % $p == 0 and next;
- my $val = valuation($p - 1, 2);
- if ($val > $k_exp and powmod($base, ($p - 1) >> ($val - $k_exp), $p) == ($congr % $p)) {
- my $n = $m * $p;
- if (($n - 1) % znorder($base, $p) == 0) {
- push @list, $n;
- }
- }
- }
- return;
- }
- foreach my $p (@{primes($lo, $hi)}) {
- $base % $p == 0 and next;
- my $val = valuation($p - 1, 2);
- $val > $k_exp or next;
- powmod($base, ($p - 1) >> ($val - $k_exp), $p) == ($congr % $p) or next;
- my $z = znorder($base, $p);
- if (gcd($m, $z) == 1) {
- __SUB__->($m * $p, lcm($L, $z), $p + 1, $k - 1, $k_exp, $congr);
- }
- }
- };
- # Case where 2^d == 1 (mod p), where d is the odd part of p-1.
- $generator->(1, 1, 2, $k, 0, 1);
- # Cases where 2^(d * 2^v) == -1 (mod p), for some v >= 0.
- foreach my $v (0 .. logint($B, 2)) {
- $generator->(1, 1, 2, $k, $v, -1);
- }
- return sort { $a <=> $b } @list;
- }
- # Generate all the squarefree strong Fermat pseudoprimes to base 2 with 3 prime factors in the range [1, 10^8]
- my $k = 3;
- my $base = 2;
- my $from = 1;
- my $upto = 1e8;
- my @arr = squarefree_strong_fermat_pseudoprimes_in_range($from, $upto, $k, $base);
- say join(', ', @arr);
- __END__
- 15841, 29341, 52633, 74665, 252601, 314821, 476971, 635401, 1004653, 1023121, 1907851, 1909001, 2419385, 2953711, 3581761, 4335241, 4682833, 5049001, 5444489, 5599765, 5681809, 9069229, 13421773, 15247621, 15510041, 15603391, 17509501, 26254801, 26758057, 27966709, 29111881, 35703361, 36765901, 37769887, 38342071, 44963029, 47349373, 47759041, 53399449, 53711113, 54468001, 60155201, 61377109, 61755751, 66977281, 68154001, 70030501, 71572957, 74329399, 82273201, 91659283, 99036001
|