comb_cyclic_2.pl

#!/usr/bin/perl

use 5.020;
use experimental qw(signatures);
use ntheory qw(forcomb carmichael_lambda divisors is_prime binomial factor euler_phi);
use Math::Prime::Util::GMP qw(is_carmichael vecprod gcd is_pseudoprime);
use List::Util qw(uniq);

#~ 1.8173 record: 64075459460541239985
#~ 1.8274 record: 14370921672011352376545
#~ 1.8670 record: 100724996094964745183793345
#~ 1.8937 record: 5981802520294676451270038145
#~ 1.9103 record: 141934960805533535384100259905
#~ 1.9103 record: 2609668563371076823446624111609234945
#~ 1.9176 record: 871005264581548962932418919531990803260747265
#~ 1.9296 record: 74875990602425226938138664024259158432269224416705
#~ 1.9332 record: 1018329989576989704159754753835889677652405111555820545
#~ 1.9483 record: 165502573617193579886078524884554371013787876466850820614145
#~ 1.9498 record: 14701083488299057530174696885922722686956286485260401577115414785
#~ 1.9540 record: 321030150905393790929751720043602006651739765349595158023943724894346115208705

#my @p = (3, 5, 17, 23, 29, 83, 89, 353, 449);
#my @p = (3, 5, 17, 23, 29);
#my @p = (3, 5, 17, 23, 29);
#my @p = (3, 5, 17, 23, 29, 83, 89);

#my @p = (3, 5, 17, 23, 29, 83, 89, 113, 197, 353, 449, 617);
#my @p = (3, 5, 17, 23, 29, 53, 83, 89, 113, 353, 449);
#my @p = (3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 353, 449, 617);
#my @p = (3, 5, 17, 23, 29, 53, 83, 89, 113);
#my @p = (3, 5, 17, 23, 29, 83, 89, 113, 197, 353, 449, 617, 3137);
#my @p = (3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 353, 449, 617, 4019, 656657);
#my @p = (3, 5, 17, 23, 29, 53, 83, 89);
#my @p = (3, 5, 17, 23, 29, 53, 83, 89, 113);
#my @p = (3, 5, 17, 23, 29, 43, 53, 89, 113, 127);
#my @p = (3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 257, 353, 617);
#my @p = (3, 5, 17, 23, 29);
#my @p = (3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 353, 449, 617);
#my @p = (3, 5, 17, 23, 29, 53, 83, 89, 113, 257, 353, 449, 617);
#my @p = (3, 5, 17, 23, 29, 53, 83, 89, 113, 1409, 3137, 10193, 16073, 23297, 88397, 896897, 1500929, 18386369);

#my @p = (3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 353, 617);     # seems promising
#my @p = (3, 5, 17, 23, 29, 53);     # seems promising
#my @p = (3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 353, 617);
#my @p = (3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 257, 353, 617);

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

# Cyclic numbers:
# 5074617663820994745
# 2278503331055626640505
# 7606851878067671122755
# 3415476493252384334116995

sub is_cyclic ($n) {
    gcd(euler_phi($n), $n) == 1;
}

#my @p = factor("3415476493252384334116995");
#push @p, 83;
my @p = (3, 5, 17, 23, 29, 53, 89);
#my @p = (3, 5, 17, 23, 29, 53, 89, 113, 197, 257, 353, 449, 617);
#my @p = (3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 257, 353, 617);
#my @p = (3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 353, 449, 617);
#my @p = (3, 5, 17, 23, 29, 53, 83, 89, 113, 353, 617);
#my @p = (3, 5, 17, 23, 29, 53, 89, 113, 353, 449, 617);
#my @p = (3, 5, 17, 23, 29, 53, 89, 113, 197, 353, 449, 617);
#my @p = (3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 353, 617);

#my @p = (3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 353, 617);
#my @p = (3, 5, 17, 23, 29, 43, 53, 89, 113, 127);

#my @p = (3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 353, 617);
#my @p = (3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 353, 617, 2549);
#my @p = (3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 353, 617);

# 46009993464351973747087455663482348080143005689116838928720068417823745
# 23016037591494880781207254677674405791538920089036061003925025850757168330033146229133479745
# 11741457044150231137356174490544808074004327146682132288353157271768338063142288065

my $P = vecprod(@p);
#my $C = "37839385943068863406967633413004957540054532539686888463944906014566240419460804270776358938980032660929917901837033235462145";
#my $C = "772459017179480479061611372132330246001039753130436193419524315193543873326133868681083905";
#my $C = "11741457044150231137356174490544808074004327146682132288353157271768338063142288065";
#my $C = "46009993464351973747087455663482348080143005689116838928720068417823745";
#my $C = "23016037591494880781207254677674405791538920089036061003925025850757168330033146229133479745";
my $C = "463149379463251167706230703520995680069543921166639491398457345";

#my $L = "1233872640";
# , , , , ,

#my $L = 6149448264960;
#my $L = 374902528;

#$L = 4352299952*128;
#$L = vecprod(2, 2, 2, 2, 2, 2, 2, 2, 7, 7, 7, 11, 13, 41);

#$L *= 7;
#$L *= 7;
#$L = 294181888;

foreach my $multiple (
    #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

    216, 420, 1380, 8064, 3960, 442800, 286200, 9666000
) {

    #foreach my $k(7, 11, 13, 19, 41, 77, 91, 133, 143, 209, 247, 287, 451, 533, 779, 1001, 1463, 1729, 2717, 3157, 3731, 5453, 5863, 8569, 10127, 19019, 41041, 59983, 70889, 111397, 779779) {
    #foreach my $k(2, 4, 7, 8, 11, 13, 14, 16, 22, 26, 28, 32, 44, 49, 52, 56, 77, 88, 91, 98, 104, 112, 143, 154, 176, 182, 196, 208, 224, 286, 308, 352, 364, 392, 416, 539, 572, 616, 637, 728, 784, 1001, 1078, 1144, 1232, 1274, 1456, 1568, 2002, 2156, 2288, 2464, 2548, 2912, 4004, 4312, 4576, 5096, 7007, 8008, 8624, 10192, 14014, 16016, 17248, 20384, 28028, 32032, 56056, 112112, 224224) {

    say "# L = $multiple";

    foreach my $k(1..1e5) {

    my $L =  $multiple * $k;

    #say "# [$multiple] lambda = $L";
    my @divisors = divisors($L);

    @divisors = grep { $_ > 1 and is_prime($_+1) } @divisors;
    @divisors = map{ $_ + 1 } @divisors;
    @divisors = grep { "@p" !~ /\b$_\b/ } @divisors;

    @divisors = grep{ is_cyclic(vecprod($_, $P)) } @divisors;

    #@divisors = grep{ $_ < 1e6 } @divisors;
    #@divisors = (@divisors[0..15], @divisors[$#divisors-15..$#divisors]);
    #@divisors = @divisors[0..5];

    @divisors = uniq(@divisors);
    #@divisors = grep {$_ > 617} @divisors;

    #@divisors = grep {  $_ < 1e6} @divisors;
    #while (scalar(@divisors) > 25) {
    #    pop @divisors;
    #}

  #  scalar(@divisors) <= 24 or next;

    #push @divisors, 127;
    #@divisors = (53, 257, 1409, 2003, 2297, 2549, 3137, 3329, 4019, 8009, 9857, 10193, 13313, 16073, 23297, 50177, 50513, 68993, 88397, 93809, 202049, 275969, 375233, 656657, 896897, 1500929, 3232769, 18386369, 22629377);

   # say "Primes = {", join(', ', @divisors), '}';

    my $len = scalar(@divisors);

  #  printf("binomial(%s, %s) = %s\n", $len, $len>>1, binomial($len, $len>>1));

    foreach my $k(1..scalar(@divisors)) {

        binomial($len, $k) > 1e5 and next;

        #say "# Testing: $k";
        forcomb {
            if (is_carmichael(vecprod(@divisors[@_],$P))) {
                say vecprod(@divisors[@_], $P);
            }
        } $len, $k;
    }
 }
}

__END__
166320 -> 3649
15120 -> 2972
110880 -> 2925
55440 -> 2828
277200 -> 2775
25200 -> 2685
196560 -> 2586
332640 -> 2572
75600 -> 2483
65520 -> 2432
131040 -> 2274
720720 -> 2217
327600 -> 2185
100800 -> 2114
45360 -> 1964
221760 -> 1957
30240 -> 1940
831600 -> 1937
257040 -> 1917
60480 -> 1899

7201568 -> 1
4352299952 -> 1
5789168 -> 1
2618528 -> 1
2700544 -> 1
285824 -> 1
1350272 -> 1

2520 -> 226
3960 -> 141
3360 -> 128
3600 -> 127
2160 -> 119
4680 -> 114
6120 -> 106
6336 -> 95
5940 -> 86
4032 -> 82
6048 -> 67
1680 -> 65
9000 -> 63
9072 -> 61
4200 -> 60
3024 -> 53
4320 -> 52
1260 -> 52
4800 -> 52
9720 -> 48

1093200 -> 4
1012480 -> 3
1568400 -> 3
19153848 -> 2
1086840 -> 2
6984864 -> 2
2495532 -> 2
1320840 -> 2
1649640 -> 2
1461564 -> 2
1754970 -> 2
1840230 -> 2
5145300 -> 2
4156860 -> 2
2716452 -> 2
3715860 -> 2
1402092 -> 2
6531960 -> 2
1236300 -> 2
1626576 -> 2