experimental-projects/pseudoprimes/generators/erdos_generate_fibonacci_from_lambdas.pl

#!/usr/bin/perl

# Erdos construction method for Carmichael numbers:
#   1. Choose an even integer L with many prime factors.
#   2. Let P be the set of primes d+1, where d|L and d+1 does not divide L.
#   3. Find a subset S of P such that prod(S) == 1 (mod L). Then prod(S) is a Carmichael number.

# Alternatively:
#   3. Find a subset S of P such that prod(S) == prod(P) (mod L). Then prod(P) / prod(S) is a Carmichael number.

# This version generates Fibonacci pseudoprimes and other Lucas pseudoprimes.

use 5.020;
use warnings;
use ntheory qw(:all);
use List::Util qw(shuffle uniq);
use experimental qw(signatures);

# Modular product of a list of integers
sub vecprodmod ($arr, $mod) {
    my $prod = 1;
    foreach my $k (@$arr) {
        $prod = mulmod($prod, $k, $mod);
    }
    $prod;
}

# Primes p such that p-1 divides L and p does not divide L
#~ sub lambda_primes ($L) {
    #~ grep { ($L % $_) != 0 } grep { $_ > 2 and is_prime($_) } map { $_ + 1 } divisors($L);
#~ }

sub lambda_primes ($L) {
    my @D = divisors($L);

    #my $D = 19;              # 5, -7, 9, -11, 13, -15, 17, -19, 21
    my $D = 5;              # 5, 12, 21, 32, 45, 60, 77, 96

    uniq(

    #~ (grep { ($_ > 2) && ($L % $_ != 0) && is_prime($_) && ($L % ($_ - kronecker(5, $_)) == 0) } map { $_ - 1 } @D),
    #~ (grep { ($_ > 2) && ($L % $_ != 0) && is_prime($_) && ($L % ($_ + kronecker(5, $_)) == 0) } map { $_ + 1 } @D)

    #~ (grep { ($_ > 2) and ($L % $_ != 0) and kronecker($D, $_) == -1 and ($L % ($_ - kronecker($D, $_)) == 0) and is_smooth($_-1, 1e2) } map { $_ + 1 } grep { is_prime($_+1) } @D),
    #~ (grep { ($_ > 2) and ($L % $_ != 0) and kronecker($D, $_) == -1 and ($L % ($_ - kronecker($D, $_)) == 0) and is_smooth($_-1, 1e2) } map { $_ - 1 } grep { is_prime($_-1) } @D)

    #~ (grep { ($_ > 2) and ($L % $_ != 0) and kronecker($D, $_) == -1 and ($L % ($_ - kronecker($D, $_)) == 0) and is_smooth($_-1, (factor($_+1))[-1]) } map { $_ + 1 } grep { is_prime($_+1) } @D),
    #~ (grep { ($_ > 2) and ($L % $_ != 0) and kronecker($D, $_) == -1 and ($L % ($_ - kronecker($D, $_)) == 0) and is_smooth($_-1, (factor($_+1))[-1]) } map { $_ - 1 } grep { is_prime($_-1) } @D)

    (grep { ($_ > 2) and ($L % $_ != 0) and (kronecker($D, $_) == -1) and ($L % ($_ - kronecker($D, $_)) == 0) } map { $_ + 1 } grep { is_prime($_+1) } @D),
    (grep { ($_ > 2) and ($L % $_ != 0) and (kronecker($D, $_) == -1) and ($L % ($_ - kronecker($D, $_)) == 0) } map { $_ - 1 } grep { is_prime($_-1) } @D)

    );
}

sub method_1 ($L) {     # smallest numbers first

    my @P = lambda_primes($L);
    my $n = scalar(@P);
    my @orig = @P;

    my $max = 1e4;
    my $max_k = 30;

    foreach my $k (3 .. @P>>1) {
        #next if (binomial($n, $k) > 1e6);

        next if ($k > $max_k);

        @P = @orig;

        my $count = 0;
        forcomb {
            if (vecprodmod([@P[@_]], $L) == 1) {
                say vecprod(@P[@_]);
            }
            lastfor if (++$count > $max);
        } $n, $k;

        next if (binomial($n, $k) < $max);

        @P = reverse(@P);

        $count = 0;
        forcomb {
            if (vecprodmod([@P[@_]], $L) == 1) {
                say vecprod(@P[@_]);
            }
            lastfor if (++$count > $max);
        } $n, $k;

        @P = shuffle(@P);

        $count = 0;
        forcomb {
            if (vecprodmod([@P[@_]], $L) == 1) {
                say vecprod(@P[@_]);
            }
            lastfor if (++$count > $max);
        } $n, $k;
    }

    my $B = $L-vecprodmod(\@P, $L);
    my $T = vecprod(@P);

    foreach my $k (2 .. @P>>1) {
        #next if (binomial($n, $k) > 1e6);

        last if ($k > $max_k);

        @P = @orig;

        my $count = 0;
        forcomb {
            if (vecprodmod([@P[@_]], $L) == $B) {
                my $S = vecprod(@P[@_]);
                say ($T / $S) if ($T != $S);
            }
            lastfor if (++$count > $max);
        } $n, $k;

        next if (binomial($n, $k) < $max);

        @P = reverse(@P);

        $count = 0;
        forcomb {
            if (vecprodmod([@P[@_]], $L) == $B) {
                my $S = vecprod(@P[@_]);
                say ($T / $S) if ($T != $S);
            }
            lastfor if (++$count > $max);
        } $n, $k;

        @P = shuffle(@P);

        $count = 0;
        forcomb {
            if (vecprodmod([@P[@_]], $L) == $B) {
                my $S = vecprod(@P[@_]);
                say ($T / $S) if ($T != $S);
            }
            lastfor if (++$count > $max);
        } $n, $k;
    }
}

while (<>) {
    #next if $_ < 1e4;
    chomp;
    #next if $_ < 49904316;
    method_1($_);
}

__END__
foreach my $L(

#2..5e3

#80, 120, 144, 2520, 5760, 6480, 7920, 15120, 30240, 94248, 110880, 285120, 597168, 604800, 1441440, 1663200, 1738800, 2217600, 5216400, 13305600, 43243200, 64864800, 648648000, 4034016000, 8951342400, 12070749600, 67541947200

#24, 36, 48, 60, 120, 180, 240, 360, 720, 840, 1260, 1680, 2520, 5040, 7560, 10080, 15120, 20160, 25200, 27720, 45360, 50400, 55440, 83160, 110880, 166320, 221760, 277200, 332640, 498960, 554400, 665280, 720720, 1081080, 1441440, 2162160

#10080, 15120, 20160, 25200, 27720, 45360, 50400, 55440, 83160, 110880, 166320, 221760, 277200, 332640, 498960, 554400, 665280, 720720, 1081080, 1441440, 2162160, 2882880, 3603600, 4324320, 6486480, 7207200, 8648640, 10810800, 14414400, 17297280, 21621600, 32432400, 36756720, 43243200, 61261200, 73513440, 110270160, 122522400, 147026880, 183783600, 245044800, 294053760, 367567200, 551350800, 698377680, 735134400, 1102701600, 1396755360, 2095133040, 2205403200, 2327925600, 2793510720, 3491888400, 4655851200, 5587021440, 6983776800, 10475665200, 13967553600, 20951330400, 27935107200, 41902660800, 48886437600, 64250746560, 73329656400, 80313433200, 97772875200, 128501493120, 146659312800, 160626866400, 240940299600, 293318625600, 321253732800, 481880599200, 642507465600, 963761198400, 1124388064800, 1606268664000, 1686582097200, 1927522396800, 2248776129600, 3212537328000, 3373164194400, 4497552259200, 6746328388800, 8995104518400, 9316358251200

#34082048, 64352288, 73545472, 128704576, 147090944, 257409152, 294181888, 374902528, 514818304, 1176727552, 1231886656, 2059273216

#11531520, 80720640, 126846720, 149909760, 329801472, 887927040, 1049368320, 1233872640, 1649007360, 2410087680, 11543051520, 16040344320, 31331139840, 40971490560, 48420852480, 106525875456, 219317978880, 286800433920, 473034481920, 532629377280, 745681128192, 3728405640960, 3979054759680, 5203379301120, 6149448264960, 67643930914560, 473507516401920

#144144, 720720, 1585584, 2306304, 15135120, 25369344, 31279248, 157549392, 1049368320, 1233872640, 1533692160, 1815061248, 2048142096, 6146300160

#5040, 7560, 10080, 21420, 57960, 59616, 109200, 110160, 163800, 178020, 308448, 540540, 811440, 3908520, 7592400, 304045280
#3780, 4032, 5040, 10080, 22176, 37800, 75600, 648000

#73545472, 128704576, 147090944, 202250048, 257409152, 294181888, 374902528, 514818304, 1029636608, 1231886656

#2130128, 2434432, 2818816, 3587584, 4260256, 4868864, 8520512, 10506496, 18386368, 19731712, 34082048, 36644608, 36772736, 53557504, 57209152, 63295232, 73545472, 87030944, 128704576, 147090944, 174061888, 202250048, 294181888, 348123776, 374902528, 618345728

#73545472, 128704576, 147090944, 202250048, 257409152, 294181888, 374902528, 514818304, 1029636608, 1231886656

#5040, 75600, 5040, 3780, 37800, 4032, 10080, 22176, 648000
#1240624, 2626624, 7771456, 10506496, 15542912, 139678448

#map{ $_ * 16 } 1188, 2024, 285824, 5789168, 7201568, 15542912, 4352299952

#map { $_ * 7 } (2130128, 2434432, 2818816, 3587584, 4260256, 4868864, 8520512, 10506496, 18386368, 19731712, 34082048, 36644608, 36772736, 53557504, 57209152, 63295232, 73545472, 87030944, 128704576, 147090944, 174061888, 202250048, 257409152, 294181888, 348123776, 374902528, 514818304, 618345728, 1029636608, 1231886656)

#divisors(2237587312658553856)

#285824, 655424, 1350272, 2618528, 2700544, 5789168, 7201568, 15542912, 34105456, 127359232, 251075968, 1426704224, 4352299952, 83510730688, 1189651277584, 38988702524464

#285824, 655424, 1350272, 2618528, 2700544, 5789168, 7201568, 15542912, 34105456, 127359232, 251075968, 1426704224, 4352299952, 83510730688, 1189651277584, 38988702524464

#(4352299952 * 32)

#655424 * 2

#2024, 1188, 7201568, 15542912, 285824, 5789168, 4352299952

) {
    method_1($L);
    #method_2($L);
}