diff --git a/Math/even_fermat_pseudoprimes_in_range.pl b/Math/even_fermat_pseudoprimes_in_range.pl index c8b77c9c..f1fe95e5 100644 --- a/Math/even_fermat_pseudoprimes_in_range.pl +++ b/Math/even_fermat_pseudoprimes_in_range.pl @@ -21,6 +21,8 @@ # PARI/GP program: # even_fermat_psp(A, B, k, base) = A=max(A, vecprod(primes(k))); (f(m, l, p, j) = my(list=List()); forprime(q=p, sqrtnint(B\m, j), if(base%q != 0, my(v=m*q, t=q); while(v <= B, my(L=lcm(l, znorder(Mod(base, t)))); if(gcd(L, v) == 1, if(j==1, if(v>=A && if(k==1, !isprime(v), 1) && (v-1)%L == 0, listput(list, v)), list=concat(list, f(v, L, q+1, j-1))), break); v *= q; t *= q))); list); vecsort(Vec(f(2, 1, 3, k-1))); +# FIXME: it doesn't generate all the terms for bases > 2. + use 5.020; use warnings; @@ -35,6 +37,8 @@ ($A, $B, $k, $base, $callback) $A = vecmax($A, pn_primorial($k)); + my %seen; + sub ($m, $L, $lo, $j) { my $hi = rootint(divint($B, $m), $j); @@ -53,8 +57,9 @@ ($A, $B, $k, $base, $callback) for (my ($q, $v) = ($p, $m * $p) ; $v <= $B ; ($q, $v) = ($q * $p, $v * $p)) { $v >= $A or next; $k == 1 and is_prime($v) and next; - ($v - 1) % znorder($base, $q) == 0 or next; - $callback->($v); + #($v - 1) % znorder($base, $q) == 0 or next; + powmod($base, $v, $v) == $base or next; + $callback->($v) if !$seen{$v}++; } } return; @@ -65,14 +70,13 @@ ($A, $B, $k, $base, $callback) $t += $L while ($t < $lo); for (my $p = $t ; $p <= $hi ; $p += $L) { - if (is_prime($p) and $base % $p != 0) { - for (my ($q, $v) = ($p, $m * $p) ; $v <= $B ; ($q, $v) = ($q * $p, $v * $p)) { - $v >= $A or next; - $k == 1 and is_prime($v) and next; - ($v - 1) % $L == 0 or next; - ($v - 1) % znorder($base, $q) == 0 or next; - $callback->($v); - } + if (is_prime_power($p) and gcd($m, $p) == 1 and gcd($base, $p) == 1) { + my $v = $m * $p; + $v >= $A or next; + $k == 1 and is_prime($v) and next; + #($v - 1) % znorder($base, $p) == 0 or next; + powmod($base, $v, $v) == $base or next; + $callback->($v) if !$seen{$v}++; } } diff --git a/Math/even_squarefree_fermat_pseudoprimes_in_range.pl b/Math/even_squarefree_fermat_pseudoprimes_in_range.pl index c9c37807..99fc8f2e 100644 --- a/Math/even_squarefree_fermat_pseudoprimes_in_range.pl +++ b/Math/even_squarefree_fermat_pseudoprimes_in_range.pl @@ -16,6 +16,8 @@ # PARI/GP program (in range) (version 2): # even_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(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(2, 1, 2, k-1))); +# FIXME: it may not generate all the terms for bases > 2. + use 5.020; use warnings; diff --git a/Math/fermat_pseudoprimes_in_range.pl b/Math/fermat_pseudoprimes_in_range.pl index 5f15f458..b5e4db8b 100644 --- a/Math/fermat_pseudoprimes_in_range.pl +++ b/Math/fermat_pseudoprimes_in_range.pl @@ -30,6 +30,8 @@ ($A, $B, $k, $base, $callback) $A = vecmax($A, pn_primorial($k)); + my %seen; + sub ($m, $L, $lo, $j) { my $hi = rootint(divint($B, $m), $j); @@ -49,7 +51,8 @@ ($A, $B, $k, $base, $callback) $v >= $A or next; $k == 1 and is_prime($v) and next; ($v - 1) % znorder($base, $q) == 0 or next; - $callback->($v); + #powmod($base, $v-1, $v) == 1 or next; + $callback->($v) if !$seen{$v}++; } } return; @@ -60,14 +63,14 @@ ($A, $B, $k, $base, $callback) $t += $L while ($t < $lo); for (my $p = $t ; $p <= $hi ; $p += $L) { - if (is_prime($p) and $base % $p != 0) { - for (my ($q, $v) = ($p, $m * $p) ; $v <= $B ; ($q, $v) = ($q * $p, $v * $p)) { - $v >= $A or next; - $k == 1 and is_prime($v) and next; - ($v - 1) % $L == 0 or next; - ($v - 1) % znorder($base, $q) == 0 or next; - $callback->($v); - } + if (is_prime_power($p) and gcd($m, $p) == 1 and gcd($base, $p) == 1) { + + my $v = $m * $p; + $v >= $A or next; + $k == 1 and is_prime($v) and next; + ($v - 1) % znorder($base, $p) == 0 or next; + #powmod($base, $v-1, $v) == 1 or next; + $callback->($v) if !$seen{$v}++; } } @@ -102,5 +105,25 @@ ($A, $B, $k, $base, $callback) say join(', ', sort { $a <=> $b } @arr); +# Run some tests + +if (0) { # true to run some tests + foreach my $k (1 .. 5) { + + my $lo = pn_primorial($k); + my $hi = mulint($lo, 1000); + my $omega_primes = omega_primes($k, $lo, $hi); + + foreach my $base (2 .. 100) { + my @this = grep { is_pseudoprime($_, $base) and !is_prime($_) } @$omega_primes; + my @that; + fermat_pseudoprimes_in_range($lo, $hi, $k, $base, sub ($n) { push @that, $n }); + @that = sort { $a <=> $b } @that; + join(' ', @this) eq join(' ', @that) + or die "Error for k = $k and base = $base with hi = $hi\n(@this) != (@that)"; + } + } +} + __END__ 91, 121, 286, 671, 703, 949, 1105, 1541, 1729, 1891, 2465, 2665, 2701, 2821, 3281, 3367, 3751, 4961, 5551, 6601, 7381, 8401, 8911, 10585, 11011, 12403, 14383, 15203, 15457, 15841, 16471, 16531, 18721, 19345, 23521, 24046, 24661, 24727, 28009, 29161, 29341, 30857, 31621, 31697, 32791, 38503, 41041, 44287, 46657, 46999, 47197, 49051, 49141, 50881, 52633, 53131, 55261, 55969, 63139, 63973, 65485, 68887, 72041, 74593, 75361, 76627, 79003, 82513, 83333, 83665, 87913, 88561, 88573, 88831, 90751, 93961, 96139, 97567 diff --git a/Math/fermat_pseudoprimes_in_range_mpz.pl b/Math/fermat_pseudoprimes_in_range_mpz.pl index 5eac59ac..d11ce0fd 100644 --- a/Math/fermat_pseudoprimes_in_range_mpz.pl +++ b/Math/fermat_pseudoprimes_in_range_mpz.pl @@ -35,6 +35,8 @@ ($A, $B, $k, $base, $callback) my $v = Math::GMPz::Rmpz_init(); my $w = Math::GMPz::Rmpz_init(); + my %seen; + sub ($m, $L, $lo, $j) { Math::GMPz::Rmpz_tdiv_q($u, $B, $m); @@ -63,24 +65,18 @@ ($A, $B, $k, $base, $callback) $t += $L while ($t < $lo); for (my $p = $t ; $p <= $hi ; $p += $L) { - if (is_prime($p) and $base % $p != 0) { + if (is_prime_power($p) and Math::GMPz::Rmpz_gcd_ui($Math::GMPz::NULL, $m, $p) == 1 and gcd($base, $p) == 1) { - Math::GMPz::Rmpz_set_ui($u, $p); Math::GMPz::Rmpz_mul_ui($v, $m, $p); - while (Math::GMPz::Rmpz_cmp($v, $B) <= 0) { - if ($k == 1 and is_prime($v)) { - ## ok - } - elsif (Math::GMPz::Rmpz_cmp($v, $A) >= 0) { - Math::GMPz::Rmpz_sub_ui($w, $v, 1); - if ((ref($L) ? Math::GMPz::Rmpz_divisible_p($w, $L) : Math::GMPz::Rmpz_divisible_ui_p($w, $L)) - and Math::GMPz::Rmpz_divisible_ui_p($w, znorder($base, $u))) { - $callback->(Math::GMPz::Rmpz_init_set($v)); - } + if ($k == 1 and is_prime($p) and Math::GMPz::Rmpz_cmp_ui($m, 1) == 0) { + ## ok + } + elsif (Math::GMPz::Rmpz_cmp($v, $A) >= 0) { + Math::GMPz::Rmpz_sub_ui($w, $v, 1); + if (Math::GMPz::Rmpz_divisible_ui_p($w, znorder($base, $p))) { + $callback->(Math::GMPz::Rmpz_init_set($v)) if !$seen{Math::GMPz::Rmpz_get_str($v, 10)}++; } - Math::GMPz::Rmpz_mul_ui($u, $u, $p); - Math::GMPz::Rmpz_mul_ui($v, $v, $p); } } } @@ -126,5 +122,25 @@ ($A, $B, $k, $base, $callback) say join(', ', sort { $a <=> $b } @arr); +# Run some tests + +if (0) { # true to run some tests + foreach my $k (1 .. 5) { + + my $lo = pn_primorial($k); + my $hi = mulint($lo, 1000); + my $omega_primes = omega_primes($k, $lo, $hi); + + foreach my $base (2 .. 100) { + my @this = grep { is_pseudoprime($_, $base) and !is_prime($_) } @$omega_primes; + my @that; + fermat_pseudoprimes_in_range($lo, $hi, $k, $base, sub ($n) { push @that, $n }); + @that = sort { $a <=> $b } @that; + join(' ', @this) eq join(' ', @that) + or die "Error for k = $k and base = $base with hi = $hi\n(@this) != (@that)"; + } + } +} + __END__ 91, 121, 286, 671, 703, 949, 1105, 1541, 1729, 1891, 2465, 2665, 2701, 2821, 3281, 3367, 3751, 4961, 5551, 6601, 7381, 8401, 8911, 10585, 11011, 12403, 14383, 15203, 15457, 15841, 16471, 16531, 18721, 19345, 23521, 24046, 24661, 24727, 28009, 29161, 29341, 30857, 31621, 31697, 32791, 38503, 41041, 44287, 46657, 46999, 47197, 49051, 49141, 50881, 52633, 53131, 55261, 55969, 63139, 63973, 65485, 68887, 72041, 74593, 75361, 76627, 79003, 82513, 83333, 83665, 87913, 88561, 88573, 88831, 90751, 93961, 96139, 97567 diff --git a/Math/squarefree_fermat_pseudoprimes_in_range.pl b/Math/squarefree_fermat_pseudoprimes_in_range.pl index e4a7eb90..80b6960f 100644 --- a/Math/squarefree_fermat_pseudoprimes_in_range.pl +++ b/Math/squarefree_fermat_pseudoprimes_in_range.pl @@ -83,5 +83,25 @@ ($A, $B, $k, $base, $callback) say join(', ', sort { $a <=> $b } @arr); +# Run some tests + +if (0) { # true to run some tests + foreach my $k (2 .. 6) { + + my $lo = pn_primorial($k); + my $hi = mulint($lo, 1000); + my @omega_primes = grep { is_square_free($_) } @{omega_primes($k, $lo, $hi)}; + + foreach my $base (2 .. 100) { + my @this = grep { is_pseudoprime($_, $base) } @omega_primes; + my @that; + squarefree_fermat_pseudoprimes_in_range($lo, $hi, $k, $base, sub ($n) { push @that, $n }); + @that = sort { $a <=> $b } @that; + join(' ', @this) eq join(' ', @that) + or die "Error for k = $k and base = $base with hi = $hi\n(@this) != (@that)"; + } + } +} + __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 diff --git a/Math/strong_fermat_pseudoprimes_in_range.pl b/Math/strong_fermat_pseudoprimes_in_range.pl index 59d47e85..21beb2e4 100644 --- a/Math/strong_fermat_pseudoprimes_in_range.pl +++ b/Math/strong_fermat_pseudoprimes_in_range.pl @@ -31,6 +31,8 @@ ($A, $B, $k, $base, $callback) $A = vecmax($A, pn_primorial($k)); $A > $B and return; + my %seen; + my $generator = sub ($m, $L, $lo, $j, $k_exp, $congr) { my $hi = rootint(divint($B, $m), $j); @@ -54,7 +56,7 @@ ($A, $B, $k, $base, $callback) $v >= $A or next; $k == 1 and is_prime($v) and next; ($v - 1) % znorder($base, $q) == 0 or next; - $callback->($v); + $callback->($v) if !$seen{$v}++; } } return; @@ -65,19 +67,17 @@ ($A, $B, $k, $base, $callback) $t += $L while ($t < $lo); for (my $p = $t ; $p <= $hi ; $p += $L) { - if (is_prime($p) and $base % $p != 0) { - my $val = valuation($p - 1, 2); - $val > $k_exp or next; - powmod($base, ($p - 1) >> ($val - $k_exp), $p) == ($congr % $p) or next; + my $val = valuation($p - 1, 2); + $val > $k_exp or next; + powmod($base, ($p - 1) >> ($val - $k_exp), $p) == ($congr % $p) or next; - for (my ($q, $v) = ($p, $m * $p) ; $v <= $B ; ($q, $v) = ($q * $p, $v * $p)) { - $v >= $A or next; - $k == 1 and is_prime($v) and next; - ($v - 1) % $L == 0 or next; - ($v - 1) % znorder($base, $q) == 0 or next; - $callback->($v); - } + if (is_prime_power($p) and gcd($m, $p) == 1 and gcd($base, $p) == 1) { + my $v = $m * $p; + $v >= $A or next; + $k == 1 and is_prime($v) and next; + ($v - 1) % znorder($base, $p) == 0 or next; + $callback->($v) if !$seen{$v}++; } } @@ -123,5 +123,25 @@ ($A, $B, $k, $base, $callback) say join(', ', sort { $a <=> $b } @arr); +# Run some tests + +if (0) { # true to run some tests + foreach my $k (1 .. 5) { + + my $lo = pn_primorial($k); + my $hi = mulint($lo, 10000); + my $omega_primes = omega_primes($k, $lo, $hi); + + foreach my $base (2 .. 100) { + my @this = grep { is_strong_pseudoprime($_, $base) and !is_prime($_) } @$omega_primes; + my @that; + strong_fermat_pseudoprimes_in_range($lo, $hi, $k, $base, sub ($n) { push @that, $n }); + @that = sort { $a <=> $b } @that; + join(' ', @this) eq join(' ', @that) + or die "Error for k = $k and base = $base with hi = $hi\n(@this) != (@that)"; + } + } +} + __END__ 121, 703, 1891, 3281, 8401, 8911, 10585, 12403, 16531, 18721, 19345, 23521, 31621, 44287, 47197, 55969, 63139, 74593, 79003, 82513, 87913, 88573, 97567 diff --git a/Math/strong_fermat_pseudoprimes_in_range_mpz.pl b/Math/strong_fermat_pseudoprimes_in_range_mpz.pl index 5a2f3df6..0363ccad 100644 --- a/Math/strong_fermat_pseudoprimes_in_range_mpz.pl +++ b/Math/strong_fermat_pseudoprimes_in_range_mpz.pl @@ -33,7 +33,8 @@ ($A, $B, $k, $base, $callback) my $u = Math::GMPz::Rmpz_init(); my $v = Math::GMPz::Rmpz_init(); - my $w = Math::GMPz::Rmpz_init(); + + my %seen; my $generator = sub ($m, $L, $lo, $j, $k_exp, $congr) { @@ -48,22 +49,6 @@ ($A, $B, $k, $base, $callback) if ($j == 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) { @@ -79,28 +64,22 @@ ($A, $B, $k, $base, $callback) $t += $L while ($t < $lo); for (my $p = $t ; $p <= $hi ; $p += $L) { - if (is_prime($p) and $base % $p != 0) { + + if (is_prime_power($p) and Math::GMPz::Rmpz_gcd_ui($Math::GMPz::NULL, $m, $p) == 1 and gcd($base, $p) == 1) { my $val = valuation($p - 1, 2); $val > $k_exp or next; powmod($base, ($p - 1) >> ($val - $k_exp), $p) == ($congr % $p) or next; - Math::GMPz::Rmpz_set_ui($u, $p); - Math::GMPz::Rmpz_mul_ui($v, $m, $p); - - while (Math::GMPz::Rmpz_cmp($v, $B) <= 0) { - if ($k == 1 and is_prime($v)) { - ## ok - } - else { - Math::GMPz::Rmpz_sub_ui($w, $v, 1); - if ((ref($L) ? Math::GMPz::Rmpz_divisible_p($w, $L) : Math::GMPz::Rmpz_divisible_ui_p($w, $L)) - and Math::GMPz::Rmpz_divisible_ui_p($w, znorder($base, $u))) { - $callback->(Math::GMPz::Rmpz_init_set($v)); - } + if ($k == 1 and is_prime($p) and Math::GMPz::Rmpz_cmp_ui($m, 1) == 0) { + ## ok + } + else { + 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))) { + $callback->(Math::GMPz::Rmpz_init_set($v)) if !$seen{Math::GMPz::Rmpz_get_str($v, 10)}++; } - Math::GMPz::Rmpz_mul_ui($u, $u, $p); - Math::GMPz::Rmpz_mul_ui($v, $v, $p); } } } @@ -157,5 +136,25 @@ ($A, $B, $k, $base, $callback) say join(', ', sort { $a <=> $b } @arr); +# Run some tests + +if (0) { # true to run some tests + foreach my $k (1 .. 5) { + + my $lo = pn_primorial($k); + my $hi = mulint($lo, 10000); + my $omega_primes = omega_primes($k, $lo, $hi); + + foreach my $base (2 .. 100) { + my @this = grep { is_strong_pseudoprime($_, $base) and !is_prime($_) } @$omega_primes; + my @that; + strong_fermat_pseudoprimes_in_range($lo, $hi, $k, $base, sub ($n) { push @that, $n }); + @that = sort { $a <=> $b } @that; + join(' ', @this) eq join(' ', @that) + or die "Error for k = $k and base = $base with hi = $hi\n(@this) != (@that)"; + } + } +} + __END__ 121, 703, 1891, 3281, 8401, 8911, 10585, 12403, 16531, 18721, 19345, 23521, 31621, 44287, 47197, 55969, 63139, 74593, 79003, 82513, 87913, 88573, 97567