perl-scripts/Math/squarefree_fermat_overpseudoprimes_in_range.pl at master · tabulon-ext/perl-scripts · GitHub

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
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
#!/usr/bin/perl

# Daniel "Trizen" Șuteu
# Date: 28 August 2022
# Edit: 04 September 2022
# https://github.com/trizen

# Generate all the squarefree Fermat overpseudoprimes 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.036;
use ntheory 0.74 qw(:all);
use Memoize qw(memoize);

memoize('inverse_znorder_primes');

sub inverse_znorder_primes ($base, $lambda) {
    my %seen;
    grep { !$seen{$_}++ } factor(subint(powint($base, $lambda), 1));
}

sub iterate_over_primes ($x, $y, $base, $lambda, $callback) {

    if ($lambda > 1 and $lambda <= 135) {
        foreach my $p (inverse_znorder_primes($base, $lambda)) {

            next if $p < $x;
            last if $p > $y;

            #znorder($base, $p) == $lambda or next;

            $callback->($p);
        }
        return;
    }

    if ($lambda > 1) {
        for (my $w = $lambda * cdivint($x - 1, $lambda) ; $w <= $y ; $w += $lambda) {
            if (is_prime($w + 1) and powmod($base, $lambda, $w + 1) == 1) {
                $callback->($w + 1);
            }
        }
        return;
    }

    for (my $p = next_prime($x - 1) ; $p <= $y ; $p = next_prime($p)) {
        $callback->($p);
    }
}

sub squarefree_fermat_overpseudoprimes_in_range ($A, $B, $k, $base, $callback) {

    $A = vecmax($A, pn_primorial($k));

    my $F;
    $F = sub ($m, $lambda, $lo, $k) {

        my $hi = rootint(divint($B, $m), $k);

        if ($lo > $hi) {
            return;
        }

        if ($k == 1) {

            $lo = vecmax($lo, cdivint($A, $m));

            if ($lo > $hi) {
                return;
            }

            iterate_over_primes(
                $lo, $hi, $base, $lambda,
                sub ($p) {
                    if (powmod($base, $lambda, $p) == 1) {
                        if (($m * $p - 1) % $lambda == 0 and znorder($base, $p) == $lambda) {
                            $callback->($m * $p);
                        }
                    }
                }
            );

            return;
        }

        iterate_over_primes(
            $lo, $hi, $base, $lambda,
            sub ($p) {
                if ($base % $p != 0) {
                    my $z = znorder($base, $p);
                    if (($z == $lambda or $lambda == 1) and gcd($z, $m) == 1) {
                        $F->($m * $p, $z, $p + 1, $k - 1);
                    }
                }
            }
        );
    };

    $F->(1, 1, 2, $k);
    undef $F;
}

# Generate all the squarefree Fermat overpseudoprimes to base 2 with 3 prime factors in the range [13421773, 412346200100]

my $k    = 3;
my $base = 2;
my $from = 13421773;
my $upto = 412346200100;

my @arr; squarefree_fermat_overpseudoprimes_in_range($from, $upto, $k, $base, sub ($n) { push @arr, $n });

say join(', ', sort { $a <=> $b } @arr);

__END__
13421773, 464955857, 536870911, 1220114377, 1541955409, 2454285751, 3435973837, 5256967999, 5726579371, 7030714813, 8493511669, 8538455017, 8788016089, 10545166433, 13893138041, 17112890881, 18723407341, 19089110641, 21335883193, 23652189937, 37408911097, 43215089153, 47978858771, 50032571509, 50807757529, 54975581389, 59850086533, 65700513721, 68713275457, 78889735961, 85139035489, 90171022049, 99737787437, 105207688757, 125402926477, 149583518641, 161624505241, 168003672409, 175303004581, 206005507811, 219687786701, 252749217641, 262106396551, 265866960649, 276676965109, 280792563977, 294207272761, 306566231341, 355774589609, 381491063773