experimental-projects/pseudoprimes/620-programs/620_strong_fermat_generator.pl

#!/usr/bin/perl

# Author: Daniel "Trizen" Șuteu
# Date: 07 October 2018
# https://github.com/trizen

# A simple algorithm for generating a subset of strong-Lucas pseudoprimes.

# See also:
#   https://oeis.org/A217120 -- Lucas pseudoprimes
#   https://oeis.org/A217255 -- Strong Lucas pseudoprimes
#   https://oeis.org/A177745 -- Semiprimes n such that n divides Fibonacci(n+1).
#   https://oeis.org/A212423 -- Frobenius pseudoprimes == 2,3 (mod 5) with respect to Fibonacci polynomial x^2 - x - 1.

use 5.020;
use warnings;
use experimental qw(signatures);

#use Math::AnyNum qw(prod powmod);
use Math::GMPz;
use Math::Prime::Util::GMP qw(vecprod);
use ntheory
  qw(forcomb forprimes kronecker divisors is_lucas_pseudoprime is_strong_lucas_pseudoprime lucas_sequence random_prime powmod valuation);
use List::Util qw(uniq);

sub fibonacci_pseudoprimes ($limit, $callback) {

    my %common_divisors;

    my $r  = random_prime(1e8);
    my $r2 = random_prime(1e9);

    die 'error' if $r <= 1e7;
    die 'error' if $r2 + 1e7 <= $r;

    my sub generate($p) {

        #~ my $P;
        #~ for (my $k = 1 ; ; ++$k) {
            #~ if (kronecker($p, $k * $k + 4) == -1) {
                #~ $P = $k;
                #~ last;
            #~ }
        #~ }

        foreach my $d (divisors($p - 1)) {
            if (powmod(2, $d, $p) == 1) {
                push @{$common_divisors{$d}}, $p;
            }
        }

        foreach my $d (divisors($p - kronecker($p, 5))) {
            if ((lucas_sequence($p, 1, -1, $d))[1] == 0) {
                push @{$common_divisors{$d}}, $p;
            }
        }
    }

    while (<>) {
        my $p = (split(' ', $_))[-1];
        $p || next;
        $p = Math::GMPz->new($p);
        generate($p);
    }

    forprimes {
        my $p = $_;
        generate($p);
    }
    1e7;

    forprimes {
        my $p = $_;
        generate($p);
    }
    $r, $r + 1e7;

    forprimes {
        my $p = $_;
        generate($p);
    }
    $r2, $r2 + 1e7;

    my %seen;

    foreach my $arr (values %common_divisors) {

        @$arr = uniq(@$arr);
        my $l = $#{$arr} + 1;

        foreach my $k (2 .. $l) {
            forcomb {
                my $nstr = vecprod(@{$arr}[@_]);
                my $n    = Math::GMPz->new($nstr);
                $callback->($n, @{$arr}[@_]) if !$seen{$nstr}++;
            }
            $l, $k;
        }
    }
}

my @pseudoprimes;

sub is_fibonacci_pseudoprime($n) {
    (lucas_sequence($n, 1, -1, $n - kronecker($n, 5)))[0] == 0;
}

fibonacci_pseudoprimes(
    1e3,
    sub ($n, @f) {

        if (is_lucas_pseudoprime($n)) {

            say "Lucas: $n";

            #push @pseudoprimes, $n;

            if (powmod(2, $n - 1, $n) == 1) {
                die "Found a BPSW counter-example: $n = prod(@f)";
            }
        }

        if (powmod(2, $n - 1, $n) == 1) {
            say "Fermat: $n";
            if (kronecker($n, 5) == -1) {
                if (is_fibonacci_pseudoprime($n)) {
                    die "Found a special pseudoprime: $n = prod(@f)";
                }
            }
        }

        #~ if (kronecker($n, 5) == -1) {
        #~ if (powmod(2, $n-1, $n) == 1) {
        #~ die "Found a Fibonacci special number: $n = prod(@f)";
        #~ }
        #~ }
    }
);

#~ @pseudoprimes = sort { $a <=> $b } @pseudoprimes;

#~ say join(', ', @pseudoprimes);

__END__
5777, 10877, 75077, 100127, 113573, 161027, 162133, 231703, 430127, 635627, 851927, 1033997, 1106327, 1256293, 1388903, 1697183, 2263127, 2435423, 2662277, 3175883, 3399527, 3452147, 3774377, 3900797, 4109363, 4226777, 4403027, 4828277, 4870847, 5208377, 5942627, 6003923, 7353917, 8518127, 9401893, 9713027, 9793313, 9922337, 10054043, 11637583, 13277423, 13455077, 13695947, 14015843, 14985833, 15754007, 16485493, 16685003, 17497127, 19168477, 20018627, 22361327, 23307377, 24157817, 25948187, 27854147, 29395277, 29604893, 30299333, 31673333, 32702723, 34134407, 34175777, 36061997, 39247393, 39850127, 40928627, 41177993, 42389027, 42525773, 47297543, 49219673, 49476377, 50075027, 51931333, 53697953, 57464207, 59268827, 62133377, 64610027, 67237883, 70894277, 73295777, 73780877, 74580767, 75239513, 75245777, 75983627, 83241013, 83963177, 85015493, 85903277, 86023943, 87471017, 90686777, 91418543, 93400277, 98385377, 104943827, 106728053, 110734667, 116853827, 117772877, 122879063, 124477513, 131017577, 131990627, 136579127, 139904627, 142593827, 144967877, 146278373, 148472347, 153256277, 154308527, 157132127, 158197577, 163578827, 166850777, 168018353, 171579883, 177991277, 179295443, 184135673, 185504633, 186003827, 192227027, 202368143, 207023087, 210089303, 211099877, 213361937, 226525883, 229206347, 231437957, 247030877, 247882963, 253755053, 254194877, 257815277, 259179527, 264250367, 264689963, 276795217, 277932113, 280075277, 284828777, 290256947, 293485877, 306219377, 311387693, 312189697, 316701527, 320234777, 334046627, 344107133, 360783793, 375578683, 376682627, 386628527, 387009737, 400091327, 400657277, 401790377, 403675973, 409245563, 420717527, 432988877, 437118527, 438894377, 439744397, 443146057, 443969063, 448504697, 450825377, 455039027, 456780193, 461700077, 461807147, 464407883, 465964127, 468245207, 469721647, 475167377, 480053573, 480891143, 485326403, 495101777, 500337713, 504097397, 523827527, 540136277, 544339637, 558030527, 562046627, 570122027, 574181327, 577647017, 583031693, 584238563, 598147577, 623709217, 634888253, 638227127, 657665777, 659936423, 664939277, 670042903, 670786877, 686258627, 691455077, 692726473, 704907377, 727615877, 729645563, 731349233, 734498627, 747587777, 768614027, 772719947, 780421277, 788342777, 797102627, 799500077, 811541327, 812957903, 825393997, 839350363, 847053323, 847887823, 856901267, 863097377, 869420473, 873933527, 878330573, 922483693, 923962577, 930039743, 961095923, 969210377, 978920627, 979805777, 985125077, 1011449753, 1015183343, 1032469817, 1034663713, 1055586377, 1085197577, 1113330077, 1171643027, 1173580127, 1194143443, 1203809777, 1218575027, 1226486627, 1230253133, 1280000357, 1295786777, 1296805127, 1308489103, 1309056527, 1326270203, 1345118777, 1364001113, 1371177527, 1387768397, 1435476803, 1437954377, 1477822433, 1524039373, 1538321777, 1541651627, 1546097027, 1561706327, 1598226167, 1599941027, 1620370127, 1663923827, 1689403127, 1749213377, 1757470643, 1770571277, 1783687127, 1826950127, 1839059627, 1885440527, 1897742027, 1909027273, 1966151713, 1986232877, 2021685077, 2031527803, 2044641377, 2056699133, 2087064527, 2093530277, 2132534777, 2152172027, 2179815377, 2189069027, 2207635127, 2231621027, 2271885527, 2273233877, 2336003647, 2382397877, 2407312577, 2411416883, 2444927627, 2509684127, 2525294777, 2535254027, 2564590757, 2630493643, 2641736327, 2660668877, 2767644017, 2774193827, 2775683777, 2807065127, 2817978767, 2832598277, 2834103827, 2837116127, 2865812777, 2887050077, 2915997527, 2985547447, 3045706127, 3059770877, 3092714627, 3174423947, 3181427027, 3226253627, 3234291377, 3279487577, 3316826083, 3331524377, 3333157127, 3342962027, 3400444277, 3424906253, 3435168827, 3501798827, 3525270527, 3558410813, 3560625077, 3582599627, 3585571907, 3601246277, 3643805027, 3690049277, 3744599777, 3752161877, 3782019617, 3797343377, 3860348777, 3881456123, 3923872577, 3966509777, 4041679277, 4068854957, 4092184277, 4097614127, 4131655097, 4141182527, 4188641627, 4223494277, 4249267577, 4256645777, 4265877527, 4484755277, 4601042627, 4622170877, 4639493627, 4681974527, 4840050077, 4887391277, 4893325127, 4938314813, 4968798827, 4998750077, 5294024413, 6039541727, 11851534697, 22200933343, 35646833933, 68055160643, 92402327687, 98831168617, 101590045727, 192900153617, 353348357933, 353833078717, 671092578683, 1118047771487, 2270927963303, 3357827162143, 3601866154427, 3703263099587, 5324864903273, 7973122223753, 8932423389707, 18846129954107, 25022761143923, 29469429987317, 29536817792327, 61561639243505213