-
Notifications
You must be signed in to change notification settings - Fork 34
/
Copy pathsquarefree_fermat_pseudoprimes_in_range_mpz.pl
125 lines (90 loc) · 5.34 KB
/
squarefree_fermat_pseudoprimes_in_range_mpz.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
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
119
120
121
122
123
124
125
#!/usr/bin/perl
# Daniel "Trizen" Șuteu
# Date: 25 February 2023
# https://github.com/trizen
# Generate all the squarefree 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
# PARI/GP program (in range):
# squarefree_fermat(A, B, k, base=2) = A=max(A, vecprod(primes(k))); (f(m, l, lo, k) = my(list=List()); my(hi=sqrtnint(B\m, k)); if(lo > hi, return(list)); if(k==1, forprime(p=max(lo, ceil(A/m)), hi, if(base%p != 0, my(t=m*p); if((t-1)%l == 0 && (t-1)%znorder(Mod(base, p)) == 0, listput(list, t)))), forprime(p=lo, hi, if (base%p != 0, my(z=znorder(Mod(base, p))); if(gcd(m, z) == 1, list=concat(list, f(m*p, lcm(l,z), p+1, k-1)))))); list); vecsort(Vec(f(1, 1, 2, k)));
# PARI/GP program (in range) (faster):
# squarefree_fermat(A, B, k, base=2) = A=max(A, vecprod(primes(k))); (f(m, l, lo, k) = my(list=List()); my(hi=sqrtnint(B\m, k)); if(lo > hi, return(list)); if(k==1, lo=max(lo, ceil(A/m)); my(t=lift(1/Mod(m,l))); while(t < lo, t += l); forstep(p=t, hi, l, if(isprime(p), my(n=m*p); if((n-1)%znorder(Mod(base, p)) == 0, listput(list, n)))), forprime(p=lo, hi, if (base%p != 0, my(z=znorder(Mod(base, p))); if(gcd(m, z) == 1, list=concat(list, f(m*p, lcm(l,z), p+1, k-1)))))); list); vecsort(Vec(f(1, 1, 2, k)));
use 5.036;
use Math::GMPz;
use ntheory qw(:all);
sub divceil ($x, $y) { # ceil(x/y)
(($x % $y == 0) ? 0 : 1) + divint($x, $y);
}
sub squarefree_fermat_pseudoprimes_in_range ($A, $B, $k, $base) {
$A = vecmax($A, pn_primorial($k));
$A = Math::GMPz->new("$A");
$B = Math::GMPz->new("$B");
my $u = Math::GMPz::Rmpz_init();
my $v = Math::GMPz::Rmpz_init();
my @list;
sub ($m, $L, $lo, $k) {
Math::GMPz::Rmpz_tdiv_q($u, $B, $m);
Math::GMPz::Rmpz_root($u, $u, $k);
my $hi = Math::GMPz::Rmpz_get_ui($u);
if ($lo > $hi) {
return;
}
if ($k == 1) {
Math::GMPz::Rmpz_cdiv_q($u, $A, $m);
if (Math::GMPz::Rmpz_fits_ulong_p($u)) {
$lo = vecmax($lo, Math::GMPz::Rmpz_get_ui($u));
}
elsif (Math::GMPz::Rmpz_cmp_ui($u, $lo) > 0) {
if (Math::GMPz::Rmpz_cmp_ui($u, $hi) > 0) {
return;
}
$lo = Math::GMPz::Rmpz_get_ui($u);
}
if ($lo > $hi) {
return;
}
Math::GMPz::Rmpz_invert($v, $m, $L);
if (Math::GMPz::Rmpz_cmp_ui($v, $hi) > 0) {
return;
}
if (Math::GMPz::Rmpz_fits_ulong_p($L)) {
$L = Math::GMPz::Rmpz_get_ui($L);
}
my $t = Math::GMPz::Rmpz_get_ui($v);
$t > $hi && return;
$t += $L * divceil($lo - $t, $L) if ($t < $lo);
for (my $p = $t ; $p <= $hi ; $p += $L) {
if (is_prime($p) and $base % $p != 0) {
Math::GMPz::Rmpz_mul_ui($v, $m, $p);
Math::GMPz::Rmpz_sub_ui($u, $v, 1);
if (Math::GMPz::Rmpz_divisible_ui_p($u, znorder($base, $p))) {
push(@list, Math::GMPz::Rmpz_init_set($v));
}
}
}
return;
}
my $t = Math::GMPz::Rmpz_init();
my $lcm = Math::GMPz::Rmpz_init();
foreach my $p (@{primes($lo, $hi)}) {
$base % $p == 0 and next;
my $z = znorder($base, $p);
Math::GMPz::Rmpz_gcd_ui($Math::GMPz::NULL, $m, $z) == 1 or next;
Math::GMPz::Rmpz_lcm_ui($lcm, $L, $z);
Math::GMPz::Rmpz_mul_ui($t, $m, $p);
__SUB__->($t, $lcm, $p + 1, $k - 1);
}
}
->(Math::GMPz->new(1), Math::GMPz->new(1), 2, $k);
return sort { $a <=> $b } @list;
}
# Generate all the squarefree Fermat pseudoprimes to base 2 with 5 prime factors in the range [100, 10^8]
my $k = 5;
my $base = 3;
my $from = 100;
my $upto = 1e8;
my @arr = squarefree_fermat_pseudoprimes_in_range($from, $upto, $k, $base);
say join(', ', @arr);
__END__
825265, 1050985, 1275681, 2113665, 2503501, 2615977, 2882265, 3370641, 3755521, 4670029, 4698001, 4895065, 5034601, 6242685, 6973057, 7428421, 8322945, 9223401, 9224391, 9890881, 10877581, 12067705, 12945745, 13757653, 13823601, 13992265, 16778881, 17698241, 18007345, 18162001, 18779761, 19092921, 22203181, 22269745, 23386441, 25266745, 25831585, 26553241, 27218269, 27336673, 27736345, 28175001, 28787185, 31146661, 32368609, 32428045, 32756581, 34111441, 34386121, 35428141, 36121345, 36168265, 36507801, 37167361, 37695505, 37938901, 38790753, 40280065, 40886241, 41298985, 41341321, 41424801, 41471521, 42689305, 43136821, 46282405, 47006785, 49084321, 49430305, 51396865, 52018341, 52452905, 53661945, 54177949, 54215161, 54651961, 55035001, 55329985, 58708761, 59586241, 60761701, 61679905, 63337393, 63560685, 64567405, 64685545, 67371265, 67994641, 68830021, 69331969, 71804161, 72135505, 72192021, 72348409, 73346365, 73988641, 74165065, 75151441, 76595761, 77442905, 78397705, 80787421, 83058481, 84028407, 84234745, 85875361, 86968981, 88407361, 88466521, 88689601, 89816545, 89915365, 92027001, 92343745, 92974921, 93614521, 93839201, 93869665, 96259681, 96386865, 96653985, 98124481, 98756281, 99551881