Several fixes and simplifications in the generation of non-squarefree Fermat pseudoprimes.

trizen · trizen · commit 7d12fc1c3302 · 2023-02-26T16:04:23.000+02:00
diff --git 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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
