experimental-projects/pseudoprimes/generators/generate_abundant_carmichael_numbers.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.

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

# Modular product of a list of integers
sub vecprodmod ($arr, $mod) {

    #~ if ($mod > ~0) {
        #~ my $prod = Math::GMPz->new(1);
        #~ foreach my $k(@$arr) {
            #~ $prod = ($prod * $k) % $mod;
        #~ }
        #~ return $prod;
    #~ }

    if ($mod < ~0) {
        my $prod = 1;
        foreach my $k(@$arr) {
            $prod = mulmod($prod, $k, $mod);
        }
        return $prod;
    }

    my $prod = 1;

    foreach my $k (@$arr) {
        $prod = Math::Prime::Util::GMP::mulmod($prod, $k, $mod);
    }

    Math::GMPz->new($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);
     grep { ($_ > 2) and (($L % $_) != 0) and is_prime($_) } map { ($_ >= ~0) ? (Math::GMPz->new($_)+1) : ($_ + 1) } divisors($L);
}

#my @prefix = (3, 5, 17, 23, 29, 53, 83, 89);
#my @prefix = (3,5,17,23, 29);

#my @prefix = (3, 5, 17, 23, 29, 53, 83, 89);
#my @prefix = (3, 5, 17, 23, 29, 53, 83, 89, 113, 353, 617, 2003, 2549);
#my @prefix = (3, 5, 17, 23, 29, 53, 83, 89, 113, 353, 617);
#my @prefix = (3, 5, 17, 23, 29, 53, 83, 89, 113, 419, 449, 617);

#my @prefix = (3, 5, 17, 23, 29, 53, 83, 89, 113, 257, 353, 449, 617, 1409, 2003);
#my @prefix = (3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 353, 449, 617, 1409, 2003);

#~ [3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 257, 353, 617, 1409, 2003, 2549, 3137, 9857, 10193, 16073, 68993, 202049, 275969, 1500929, 18386369]
#~ [3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 353, 449, 617, 2003, 2297, 2549, 3137, 9857, 10193, 68993, 88397, 93809, 5850209, 8044037, 18386369]
#~ [3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 353, 449, 617, 2003, 2549, 3137, 4019, 4289, 10193, 16073, 21977, 38459, 50513, 52529, 76649, 93809, 97553]
#~ [3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 353, 449, 617, 1409, 2003, 2549, 8009, 9857, 23297, 50513, 68993, 275969, 375233, 1500929, 3232769, 18386369]
#~ [3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 353, 449, 617, 1409, 2003, 2297, 2549, 3329, 8009, 10193, 16073, 23297, 50177, 93809, 202049, 275969, 656657, 18386369]
#~ [3, 5, 17, 23, 29, 53, 83, 89, 113, 257, 353, 449, 617, 1409, 2003, 2297, 2549, 3137, 3329, 4019, 7547, 9857, 10193, 16073, 17837, 23297, 68993, 88397, 93809, 202049, 896897, 1500929, 2475089, 8044037, 18386369]

foreach my $n(

#qw(327976184723917106922696038500316850197097366397296955613721508063986316871814680735 685142249888262836361512024427161900061736398403953340277064230345667415945220868055415 1447705574013899373231874907614593094830449009827553408005436718720395249892251694201091895 3154550445776286734272255423692198353635548392414238876043846610091741249515216441664179239205 7498366409610233567365151142116355486591698528768645808356223392188068950097669481835754051590285 19443264100119335640177836911507709776732274285097098581067687255943662787603256966400110255773609005 6752147266492648336087178462400067734188627847011899231029509812418215 8622492059311111925183326896484886496558877760634195318024684030458060555 11183372200926512166962774984740897786036864455542551327478015187504104539835 15757371431105455643250549953499924980525942017859454820416523399193283296627515 22895460689396227049643049082435390996704193751949787854065208499027840629999779295 4065021504672018226118572207226948264489296099895 2126006246943465532260013264379693942327901860245085 1226705604486379612114027653547083404723199373361414045 837839927864197275073880887372657965425945172005845792735 8400431334513514794031775534168508142396705 3334971239801865373230614887064897732531491885 1444042546834207706608856246099100718186135986205 11304017840824088923079122979486080913570205802842055 8692789719593724381847845571224796222535488262385540295 8075601649502569950736648535667835690735468595756166934055 9020447042494370634972836414340972466551518421459638465339435 13972672468823780113572923605814166350688302034840979982810784815 21895177758646863437968771290310798671528569288595815633064499805105 45739026337813297721916763225459258424823181243876658857471740092864345 108721665604982208684996146186916657275804701816694818104210326200738548065 207093460403925930428246384004848306979999376211770028014102314959732822821308183615 14116844787629511288659115050445171405 6338463309645650568607942657649881960845 )

#qw( 269616807542440095880838902934479765646815 39623553877269279114765582013213715065835439560755)

#qw(80457526688110964090941908021307082343637031298345387811989511961940980765)
#qw(156400339552952614158261462316373572492517490659424507963838717745)
#qw(30300520893268092108675921598650019435254761657479255573181729185)
#qw(23727894199896704861923196240133139730035052198495893166156405)
#qw(5560414251522943074360489383944833466511209118626077678018703909135127547157897173895)
#qw(14807987984119898011183681250338047627675234576733585022035)
#qw(1469693990779095715245182961386962460712186859726612101336192430649765)
#qw(10094456962219244277998564908104312836457170747058383045)
#qw(48976162706891785856850267556951183909509591565)

#qw(232256510885289043333840440634089959783515)
#qw(76901951875930663768872450745768943476147523769184184245)
#qw(93426807299183714704340654128147124952637393355)
qw(100567069213330155763552910794560952586261995)

#qw(2291345950693302094766820752641350945890207963412715)

#qw(327976184723917106922696038500316850197097366397296955613721508063986316871814680735 685142249888262836361512024427161900061736398403953340277064230345667415945220868055415 1447705574013899373231874907614593094830449009827553408005436718720395249892251694201091895 3154550445776286734272255423692198353635548392414238876043846610091741249515216441664179239205 7498366409610233567365151142116355486591698528768645808356223392188068950097669481835754051590285 19443264100119335640177836911507709776732274285097098581067687255943662787603256966400110255773609005 6752147266492648336087178462400067734188627847011899231029509812418215 8622492059311111925183326896484886496558877760634195318024684030458060555 11183372200926512166962774984740897786036864455542551327478015187504104539835 15757371431105455643250549953499924980525942017859454820416523399193283296627515 22895460689396227049643049082435390996704193751949787854065208499027840629999779295 80175683972122511383731338517508443634187270777707283918484112701538301509380205 125635296784315975338307007456935731174771453308667313900264604603310518465198781235 224007734166435384028201394295716408684617501249353820684171790007702654423449426942005 467952156673683517234912712683751577742165960109900131409234869326090845090585852881848445 988782907051493271917370561900767083769196673712218977667713278886029955676407907139345764285 2154557954465203839507950454381771475533079552018925152337947234692659273418892829656634420377015 5121384257763789526510398230065470797342130095148985087107300576864451092916708256093820017236164655 13279749380381506242241462610559765777508143336721318330869230395809521683933024508051275304693374950415 4065021504672018226118572207226948264489296099895 2126006246943465532260013264379693942327901860245085 1226705604486379612114027653547083404723199373361414045 837839927864197275073880887372657965425945172005845792735 19560340012826425693128516437578862635 8400431334513514794031775534168508142396705 3334971239801865373230614887064897732531491885 1444042546834207706608856246099100718186135986205 11304017840824088923079122979486080913570205802842055 8692789719593724381847845571224796222535488262385540295 8075601649502569950736648535667835690735468595756166934055 9020447042494370634972836414340972466551518421459638465339435 13972672468823780113572923605814166350688302034840979982810784815 21895177758646863437968771290310798671528569288595815633064499805105 45739026337813297721916763225459258424823181243876658857471740092864345 108721665604982208684996146186916657275804701816694818104210326200738548065 17277031840122951876810012573270045985 99436691737018840272248469415851811397274754265692619259832323874449185818279122553922438597548102645740113148322901071530942568237363765081022282605920356808140142566166169841228255378791958207751846269805134897319182266879091065192147432933981465795852869698389629543929 895029662324906581290508473212082154386870063145499265957750747193917121550330382107855869816530471914306758448054432544850014056704511249494281565735889131630069423238061694740895526664506415827974368274516019210769959584178698677794519043838767173628471680155205055524904929 8306770296037457980957209139881334474864541056053378687353884684706744805108616276343010327767219309836681025156393188448752980460274568906556427211594787030658674317072450588890251382973284045299430111955783174295155994900762502428610931245867598138445845663520458120326642646049 78058720471863992647054894287464900060302092303733599525064454382189280933605667148795268050028559854535291593394626791852931757385200124014910746507356213727099562557529818183801692245799950173678744762048494488851580884082465235321656920917417819706975611700101744956709460944922453 749441775250366193404374040053950505478960388208146289040143826523399286243548010295583368548324203163393334588181811828579997802655306390667158077217127007993882900114843784382680047251925321617489628460427595587464028068075748724323228097728128487006672847932676853329367534532200471253 207093460403925930428246384004848306979999376211770028014102314959732822821308183615 14116844787629511288659115050445171405 6338463309645650568607942657649881960845 83562900542610185478464992069940404422726948660826302385114906951342076455803610529344640789 50221303226108721472557460234034183058058896145156607733454059077756587949937969928136129114189 32191855367935690463909332010015911340215752429045385557144051868841972875910238723935258762195149 73447150516193077029324563940188527242280274750433560293457441337357999076671648487217809)
#qw(1568442344968558041435703706985298315 622671610952517542449974371673163431055 269616807542440095880838902934479765646815)

#qw(80457526688110964090941908021307082343637031298345387811989511961940980765)
#qw(183946261763273755985808210039437380929387193716936037921324457691523453096923545)
#qw(80457526688110964090941908021307082343637031298345387811989511961940980765)

#qw(68671670497867034860945549416314930790009427315 39623553877269279114765582013213715065835439560755 28806323668774765916434578123606370852862364560668885 22152062901287794989738190577053299185851158347154372565 20579266435296361545466779046082514943655726104506412112885 23727894199896704861923196240133139730035052198495893166156405 30300520893268092108675921598650019435254761657479255573181729185 39299775598568715464952670313449075207525425869750594478416702752945 55373383818383320090118312471649746967403325050478587620089134178899505 80457526688110964090941908021307082343637031298345387811989511961940980765 117387531437953896608684243803087033139366428664285920817692697952471890936135 183946261763273755985808210039437380929387193716936037921324457691523453096923545)
#qw(100567069213330155763552910794560952586261995)

#qw(39623553877269279114765582013213715065835439560755)
#qw(23727894199896704861923196240133139730035052198495893166156405 30300520893268092108675921598650019435254761657479255573181729185 39299775598568715464952670313449075207525425869750594478416702752945 55373383818383320090118312471649746967403325050478587620089134178899505 80457526688110964090941908021307082343637031298345387811989511961940980765)

) {

    my @prefix = factor($n);
#my @prefix = (3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 257, 353, 449, 617, 1409, 2003, 2297, 2549, 3137, 3329, 4019);
#my @prefix = (5, 7, 13, 17, 19, 23, 37, 59, 67, 73, 89, 97, 109, 163, 193, 199, 233, 257, 349, 353, 397, 433, 487, 523, 577, 727, 769, 929, 1153, 1277, 1297, 1409, 1453, 1459, 1567, 1783, 2089, 2113, 2179, 2377, 2593);
#~ [3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 353, 449, 617, 2003, 2297, 2549, 3137, 9857, 10193, 68993, 88397, 93809, 5850209, 8044037, 18386369]
#~ [3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 353, 449, 617, 2003, 2549, 3137, 4019, 4289, 10193, 16073, 21977, 38459, 50513, 52529, 76649, 93809, 97553]
#~ [3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 353, 449, 617, 1409, 2003, 2549, 8009, 9857, 23297, 50513, 68993, 275969, 375233, 1500929, 3232769, 18386369]
#~ [3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 353, 449, 617, 1409, 2003, 2297, 2549, 3329, 8009, 10193, 16073, 23297, 50177, 93809, 202049, 275969, 656657, 18386369]
#~ [3, 5, 17, 23, 29, 53, 83, 89, 113, 257, 353, 449, 617, 1409, 2003, 2297, 2549, 3137, 3329, 4019, 7547, 9857, 10193, 16073, 17837, 23297, 68993, 88397, 93809, 202049, 896897, 1500929, 2475089, 8044037, 18386369]

#my @prefix = (3, 5, 17, 23, 29);
my $prefix_prod = Math::GMPz->new(vecprod(@prefix));

my sub method_1 ($L) {     # smallest numbers first

    (vecall { ($L % ($_-1)) == 0 } @prefix) or return;

    my @P = lambda_primes($L);

    #@P = grep{$_ < 1e5 } @P;

    #~ @P = grep {
        #~ (($prefix_prod % $_) == 0)
            #~ ? 1
            #~ : Math::Prime::Util::GMP::gcd(Math::Prime::Util::GMP::totient(Math::Prime::Util::GMP::mulint($prefix_prod, $_)), Math::Prime::Util::GMP::mulint($prefix_prod, $_)) eq '1'
    #~ } @P;

    #vecprodmod(@P, 3*5*17*23) == 0 or return;
    #vecprodmod(\@P, 3*5*17*23*29) == 0 or return;
    if (@prefix) {
        vecprodmod(\@P, $prefix_prod) == 0 or return;
    }
    #return if (vecprod(@P) < ~0);

    #@P = grep { $_ > $prefix[-1] } @P;
    @P = grep { gcd($prefix_prod, $_) == 1 } @P;

    say "# Testing: $L -- ", scalar(@P);

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

    my $max = 1e6;
    my $max_k = scalar(@prefix)>>1;

    my $L_rem = invmod($prefix_prod, $L);

    foreach my $k (1 .. @P) {

        #next if (binomial($n, $k) > 1e6);

        next if ($k > $max_k);

        @P = @orig;

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

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

        @P = reverse(@P);

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

        @P = shuffle(@P);

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

    my $B = vecprodmod([@P, $prefix_prod], $L);
    my $T = Math::GMPz->new(Math::Prime::Util::GMP::vecprod(@P));

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

        last if ($k > $max_k);

        @P = @orig;

        my $count = 0;
        forcomb {
            if (vecprodmod([@P[@_]], $L) == $B) {
                my $S = Math::GMPz->new(Math::Prime::Util::GMP::vecprod(@P[@_]));
                say vecprod($prefix_prod, $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 = Math::GMPz->new(Math::Prime::Util::GMP::vecprod(@P[@_]));
                say vecprod($prefix_prod, $T / $S) if ($T != $S);
            }
            lastfor if (++$count > $max);
        } $n, $k;

        @P = shuffle(@P);

        $count = 0;
        forcomb {
            if (vecprodmod([@P[@_]], $L) == $B) {
                my $S = Math::GMPz->new(Math::Prime::Util::GMP::vecprod(@P[@_]));
                say vecprod($prefix_prod, $T / $S) if ($T != $S);
            }
            lastfor if (++$count > $max);
        } $n, $k;
    }
}

use Math::GMPz;

#my %seen;
#my $count = 0;

#method_1(205058304)
method_1(carmichael_lambda($prefix_prod));

}

__END__
while (<>) {
    chomp(my $n = $_);

    #$n > 1e5 or next;

    #$n *= 656656;
    #$n *= 78848;
    #$n *= 1232;

    #$n = lcm($n, 29586292736);
    #$n = lcm($n, 1232);
    #$n = lcm($n, 5789168, 147090944);
    #$n = lcm($n, 78848);
    #$n = lcm($n, 656656);
    #$n = lcm($n, 9193184);

    #$n = lcm($n, 10506496);
    #$n = lcm($n, 36772736);

    #~ next if ($n < 707981814540);
    #~ next if ($n > 44351949725003712);

    #next if ($n < 1e6);

    #next if ($n < 1e8);
    #next if ($n < 7813080);        # for 2^64

    #next if ($n < ~0);

    #next if ($n < 17125441200);
    #next if (length($n) > 45);

   # if (++$count % 1000 == 0) {
        #say "Testing: $n";
       # $count = 0 ;
   # }

    #say "Testing: $n";

    if ($n > ~0) {
        $n = Math::GMPz->new($n);
    }

    next if $seen{$n}++;
    method_1($n);
}