experimental-projects/oeis-research/Amiram Eldar/Abundant Carmichael numbers/y.pl

#!/usr/bin/perl

use 5.014;

use Math::GMPz;
use List::Util qw(uniq);
use ntheory qw(forsemiprimes forprimes factor forsquarefree random_prime euler_phi );
use Math::Prime::Util::GMP qw(mulint is_pseudoprime vecprod divint sqrtint divisors gcd is_carmichael);

# 2435279804766429985, 1114090435183709985, 7626934787034827265, 5092454223678002545, 13294381492836585505, 1886616373665, 5473367446820370945, 180950795673242145, 3193231538989185, 916541202603455265, 11947816523586945, 101817952350880305, 1177908521713261185, 171800042106877185

# 9744543540935384545, 10186034230241370865, 10885189367088417745, 11356964096926592145, 11595472526856600705, 12275497073183645665, 12376763082073211185, 12525679299253814785, 13294381492836585505, 14631273312068834785, 15131182782493999585, 15157968873827062465, 16888413473957767105, 17215556534723568385, 17445622572008652385, 18329749222420733505

#~ 12376763082073211185
#~ 1180617122678545373185
#~ 345319996314995065
#~ 292561243007134465

#~ 12376763082073211185
#~ 1180617122678545373185
#~ 1399093241875807699585
#~ 573852931410299905
#~ 1370294210266345
#~ 9768827424988773505
#~ 29582577574410161665
#~ 305244070146997345
#~ 5486206558145814865
#~ 345319996314995065
#~ 292561243007134465
#~ 63238598993791585

#~ 9768827424988773505
#~ 32379807987407561422105
#~ 4281332081515547439745
#~ 12376763082073211185
#~ 29582577574410161665
#~ 1227189372830785
#~ 8074371355371685

#~ 1399093241875807699585
#~ 219588765971746167938305
#~ 3904542189035075088865
#~ 1031925092561764182145
#~ 1533250834894400523265

#~ 32379807987407561422105
#~ 66497283692022179297665
#~ 828696209007880663751665
#~ 242746392608701345
#~ 104374820521710440065

#~ 4281332081515547439745
#~ 11956408212534751066465
#~ 1370294210266345

# 2333379336546216408131111533710540349903201, 130912961974316767723865201454187955056178415601, 60977817398996785

my %seen;
my $psp = "828696209007880663751665";

my $limit = sqrtint($psp);

#my $gcd = 7313655;
#my $gcd = 4381310850645;
#my $gcd = 174766108521378405;
#my $gcd = 77728835801292945;
#my $gcd = 113375180882140665;
#my $gcd = 495088126122885;
#my $gcd = 34498510635;
#my $gcd = "19976310800932286865";       # record
#my $gcd = "7051637712729097263345";
my $gcd = vecprod(grep{ $_<50 }factor($psp));
#my $gcd = vecprod(5, 7, 13, 17, 19, 23);
#my $gcd = "119241686273722855330215";

my @squarefree;

forsquarefree {
    if ($_ % 2 == 1 and $_ > 1 and gcd($_, $psp) == 1 and gcd(euler_phi($_), $_) == 1) {
        push @squarefree, $_;
    }
} 1e7;

#push @squarefree, grep{ $_>1 } divisors("6283365899669117794965");
#push @squarefree, grep{ $_>1 } divisors("19976310800932286865");
#push @squarefree, grep{$_>1} divisors("480185024160041299368309995874402451");
push @squarefree, grep{$_>1 and gcd($_, $psp) == 1 and gcd(euler_phi($_), $_) == 1 } divisors("119241686273722855330215");

@squarefree = uniq(@squarefree);

if (is_carmichael($psp)) {
    say $psp;
}

foreach my $p(divisors($psp)) {

    $p > 1 or next;
    last if ($p > $limit);
    next if (gcd($p, $gcd) > 1);

    my $n = divint($psp, $p);

    for (@squarefree) {
        if (is_carmichael(mulint($n, $_))) {
            say mulint($n, $_) if !$seen{mulint($n, $_)}++;
        }
    }
}

__END__
foreach my $n(1..1e6) {

    my @copy = @factors;
    #$copy[-1] =  3*$n+1;
    #$copy[-1] =  2*$n+1;
    #$copy[-2] =  4*$n+1;

    #$copy[-2] =  2*$n+1;
    $copy[rand @copy] =  random_prime(1e4);
    $copy[rand @copy] =  random_prime(1e4);

    #$copy[rand @copy] =  random_prime(1e4);
    #$copy[rand @copy] =  2*int(rand(1e3))+1;

    #$copy[-3] =  2*int(rand(1e9))+1;
    #$copy[-2] = 3*$n+1;
    #$copy[-1] = 1;
    #$copy[rand(@copy-9)+8] =

    if (is_pseudoprime(vecprod(@copy), 2)) {
        if (vecprod(@copy) ne $psp) {
            say vecprod(@copy);
        }
    }

}

exit;

my @psp = qw(

4788772759754985
38353281032877674865
2638294879881771254145
24773130788808935997465
39898386437874778336787265
85866492509341408342261785
198408004487570768403697185

);

forsquarefree {
    #if ($_ % 80 == 9) {

    if ($_ % 2 == 1 and $_ > 1) {
        foreach my $n(@psp) {
            if (is_pseudoprime(mulint($n, $_), 2)) {
                say mulint($n, $_);
            }

            #~ my ($p, $q) = factor($_);
            #~ $n = Math::GMPz->new($n);

            #~ if ($n % $p == 0) {
                #~ if (is_pseudoprime(($n/$p)*$q, 2)) {
                    #~ say (($n/$p)*$q);
                #~ }
            #~ }

            #~ if ($n % $q == 0) {
                #~ if (is_pseudoprime(($n/$q)*$p, 2)) {
                    #~ say (($n/$q)*$p);
                #~ }
            #~ }
        }
    }

} 1e7;

__END__