123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283 |
- #!/usr/bin/perl
- # Author: Daniel "Trizen" Șuteu
- # Date: 28 January 2019
- # https://github.com/trizen
- # A new algorithm for generating Fermat pseudoprimes to base 2.
- # See also:
- # https://oeis.org/A050217 -- Super-Poulet numbers: Poulet numbers whose divisors d all satisfy d|2^d-2.
- # https://oeis.org/A214305 -- Fermat pseudoprimes to base 2 with two prime factors.
- # See also:
- # https://en.wikipedia.org/wiki/Fermat_pseudoprime
- # https://trizenx.blogspot.com/2018/08/investigating-fibonacci-numbers-modulo-m.html
- use 5.020;
- use warnings;
- use strict;
- use ntheory qw(:all);
- use experimental qw(signatures);
- use List::Util qw(uniqnum);
- # a(n) is the smallest odd composite k such that prime(n)^((k-1)/2) == -1 (mod k) and b^((k-1)/2) == 1 (mod k) for every natural b < prime(n).
- # a(n) : 3277, 5173601, 2329584217, 188985961, ... Is this sequence infinite?
- # I found upper-bounds only for the next 5 terms:
- # a(5) <= 32203213602841
- # a(6) <= 323346556958041
- # a(7) <= 2528509579568281
- # a(8) <= 5189206896360728641
- # a(9) <= 12155831039329417441
- # Notice that the upper-bounds a(5)-a(9) are all semiprimes and are of the form p * (2*(p-1) + 1), for some prime number p.
- # I couldn't find an upper-bound for a(10), which may be larger than 2^64 (assuming that it exists).
- # More generally, we can look for numbers of the form k = p*q, where p and q are odd primes, such that d|((p-1)/2) and d|((q-1)/2), satisfying:
- # prime(n)^d == -1 (mod p)
- # prime(n)^d == -1 (mod q)
- # This can be implemented efficiently by iterating over the odd primes bellow a certain limit and grouping primes together by divisors `d` that satisfy the wanted property.
- # When a group has two or more primes, we iterate over all the combinations of two of those primes, say p and q, and generate k = p*q.
- # Once we have k = p*q, it's very likely that:
- # prime(n)^((k-1)/2) == -1 (mod k)
- # For each such number k, we check if:
- # b^((k-1)/2) == 1 (mod k) for every natural b < prime(n)
- # Several more upper-bounds (example), with possible gaps, for each a(2)-a(7):
- # a(2) <= { 5173601, 6787327, 15978007, 314184487, 581618143, 682528687, 717653129, 1536112001, 1854940231, 1860373241, 2693739751, 3344191241, 4916322919, ... }
- # a(3) <= { 2329584217, 4394140633, 18491988217, 71430713857, 146463308377, 306190965697, 542381064457, 662181331897, 1749760817617, 2600044182457, 3420965597017, ... }
- # a(4) <= { 188985961, 922845241, 41610787921, 55165631041, 85103495641, 130173194161, 178420380841, 451403438041, 577255246921, 669094855201, 842611640041, 1038185022241, ... }
- # a(5) <= { 32203213602841, 92223430598881, 181122243848281, 244362410507881, 323731175307121, 359643366744481, 383326325785201, 503561507122801, ... }
- # a(6) <= { 323346556958041, 449287527233161, 871447814622001, 2419909071073201, ... }
- # a(7) <= { 2528509579568281, 51558565269914641, 251985537187183801, ... }
- # Let b(n) be the smallest odd composite k such that q^((k-1)/2) == -1 (mod k) for every prime q <= prime(n).
- #~ b(1) = 3277
- #~ b(2) = 1530787
- #~ b(3) <= 3697278427
- #~ b(4) <= 118670087467
- #~ b(5) <= 2152302898747
- #~ b(6) <= 614796634515444067
- #~ b(7) <= 614796634515444067
- my $N = 5;
- my $multiplier = 4;
- my $P = nth_prime($N);
- my $MAX = ~0;
- my $LIMIT = 1e5;
- use Math::GMPz;
- sub prod {
- my $prod = 1;
- foreach my $n(@_) {
- $prod *= $n;
- }
- $prod;
- }
- # a(5) = 11530801
- # a(7) = 15656266201
- my @primes = @{primes(100)};
- my @primes_bellow = @{primes($P)};
- sub non_residue {
- my ($n) = @_;
- foreach my $p (@primes) {
- #if ($p > $Q) {
- # return -1;
- # }
- sqrtmod($p, $n) // return $p;
- }
- return -1;
- }
- #~ #139309114031, a(9) <= 7947339136801, a(10) <= 72054898434289.
- #~ say non_residue(139309114031, nth_prime(8));
- #~ say non_residue(7947339136801, nth_prime(9));
- #~ say non_residue(72054898434289, nth_prime(10));
- #~ __END__
- #~ say non_residue(72054898434289, nth_prime(10));
- #~ __END__
- sub isok {
- my ($k) = @_;
- my $t = $k-1;
- my $u = $t>>1;
- vecall { powmod($_, $u, $k) == $t } @primes_bellow;
- }
- #say isok(614796634515444067);
- #say isok(8614572538322761627);
- #392044843, 1568179369
- #~ foreach my $p(factor(2152302898747)) {
- #~ foreach my $d(divisors(($p-1)>>1)) {
- #~ if (vecall {my $t = powmod($_, $d, $p); ($t == $p-1) } @primes_bellow) {
- #~ say "$p : $d -> ", join(', ', divisors($d));
- #~ }
- #~ }
- #~ say '';
- #~ }
- # Let a'(n) be the smallest odd prime p such that b^((p-1)/2) == 1 (mod p) for every natural b < prime(n).
- # a(n) is the least number such that the n-th prime is the least quadratic nonresidue (not necessarily coprime) modulo a(n).
- #~ __END__
- #~ forprimes {
- #~ if (is_prime(($_-1)*$multiplier + 1)) {
- #~ my $k = Math::GMPz->new($_)*(($_-1)*$multiplier + 1);
- #~ if (isok($k)) {
- #~ die "a($N) <= $k\n";
- #~ }
- #~ }
- #~ } 1e11, 1e12;
- #~ __END__
- sub fermat_pseudoprimes ($limit, $callback) {
- say "Collecting primes...";
- my %common_divisors;
- forprimes {
- my $p = $_;
- my $common = 0;
- foreach my $d (divisors(($p-1)>>1)) {
- #if (powmod($P, $d, $p) == $p-1 ) {
- if (vecall { powmod($_, $d, $p) == $p-1 } @primes_bellow) {
- #push @{$common_divisors{$d}}, $p;
- foreach my $g (divisors($d)) {
- if ($g * 1000 >= $d) {
- push @{$common_divisors{$g}}, $p;
- }
- }
- }
- }
- #if (@div) {
- # push @{$common_divisors{"@div"}}, $p;
- #}
- #~ if ($common) {
- #~ foreach my $n(2..100) {
- #~ if (is_prime(($p-1)*$n + 1)) {
- #~ my $q = ($p-1)*$n + 1;
- #~ foreach my $d (divisors($q - 1)) {
- #~ if (powmod($P, $d, $q) == $q-1) {
- #~ push @{$common_divisors{$d}}, $q;
- #~ }
- #~ }
- #~ }
- #~ }
- #~ }
- # }
- } 3, $limit;
- say "Generating combinations...";
- #~ forprimes {
- #~ my $p = $_;
- #~ foreach my $d (divisors($p - 1)) {
- #~ if (powmod($P, $d, $p) == $p-1) {
- #~ push @{$common_divisors{$d}}, $p;
- #~ }
- #~ }
- #~ } $limit>>1, $limit;
- my %seen;
- foreach my $value (values %common_divisors) {
- my $arr = [uniqnum(@$value)];
- my $l = $#{$arr} + 1;
- if ($l > 100) {
- say "Combinations: $l";
- }
- foreach my $k (2..$l) {
- forcomb {
- my $n = prod(@{$arr}[@_]);
- if ($n <= $MAX) {
- $callback->($n) #if !$seen{$n}++;
- }
- } $l, $k;
- }
- }
- }
- #~ sub is_fermat_pseudoprime ($n, $base) {
- #~ powmod($base, $n - 1, $n) == 1;
- #~ }
- #~ sub is_fibonacci_pseudoprime($n) {
- #~ (lucas_sequence($n, 1, -1, $n - kronecker($n, 5)))[0] == 0;
- #~ }
- my @values;
- my %seen;
- fermat_pseudoprimes(
- $LIMIT,
- sub ($k) {
- # is_fermat_pseudoprime($n, 2) || die "error for n=$n";
- #~ if (kronecker($n, 5) == -1) {
- #~ if (is_fibonacci_pseudoprime($n)) {
- #~ die "Found a special pseudoprime: $n = prod(@f)";
- #~ }
- #~ }
- #~ if (powmod($P, ($k-1)>>1, $k) != $k-1) {
- #~ say "false: $k";
- #~ }
- if (isok($k) and !$seen{$k}++) {
- say "Found: $k";
- push @values, $k;
- }
- # push @pseudoprimes, $n;
- }
- );
- say "a($N) <= { ", join(', ', sort { $a <=> $b} @values), " }";
- say "Min: a($N) <= ", vecmin(@values);
|