experimental-projects/highly-composite-numbers/hcn_plus_minus_1_is_prime.pl

#!/usr/bin/perl

use 5.020;
use strict;
use warnings;

use Math::GMPz;
use POSIX qw(ULONG_MAX);
use ntheory qw(:all);
use IO::Uncompress::Bunzip2;
use experimental qw(signatures declared_refs);

local $| = 1;
prime_precalc(1e7);

sub primality_pretest ($n) {

    # Must be positive
    (Math::GMPz::Rmpz_sgn($n) > 0) || return;

    # Check for divisibilty by 2
    if (Math::GMPz::Rmpz_even_p($n)) {
        return (Math::GMPz::Rmpz_cmp_ui($n, 2) == 0);
    }

    # Return early if n is too small
    Math::GMPz::Rmpz_cmp_ui($n, 101) > 0 or return 1;

    # Check for very small factors
    if (ULONG_MAX >= 18446744073709551615) {
        Math::GMPz::Rmpz_gcd_ui($Math::GMPz::NULL, $n, 16294579238595022365) == 1 or return 0;
        Math::GMPz::Rmpz_gcd_ui($Math::GMPz::NULL, $n, 7145393598349078859) == 1  or return 0;
    }
    else {
        Math::GMPz::Rmpz_gcd_ui($Math::GMPz::NULL, $n, 3234846615) == 1 or return 0;
    }

    # Size of n in base-2
    my $size = Math::GMPz::Rmpz_sizeinbase($n, 2);

    # When n is large enough, try to find a small factor (up to 10^8)
    if ($size > 10_000) {

        state %cache;
        state $g = Math::GMPz::Rmpz_init_nobless();

        my @checks = (1e4);

        push(@checks, 1e6) if ($size > 15_000);
        push(@checks, 1e7) if ($size > 20_000);
        push(@checks, 1e8) if ($size > 30_000);

        my $prev;

        foreach my $k (@checks) {

            my $primorial = (
                $cache{$k} //= do {
                    my $z = Math::GMPz::Rmpz_init_nobless();
                    Math::GMPz::Rmpz_primorial_ui($z, $k);
                    Math::GMPz::Rmpz_divexact($z, $z, $prev) if defined($prev);
                    $z;
                }
            );

            Math::GMPz::Rmpz_gcd($g, $primorial, $n);

            if (Math::GMPz::Rmpz_cmp_ui($g, 1) > 0) {
                return 0;
            }

            $prev = $primorial;
        }
    }

    return 1;
}

sub primorial_inflation ($n) {
    state %primorial;

    my $tmp  = Math::GMPz->new(1);
    my $prod = Math::GMPz->new(1);

    foreach my \@pp (factor_exp($n)) {
        my ($p, $e) = @pp;

        my $prim = $primorial{$p} //= do {
            my $z = Math::GMPz::Rmpz_init_nobless();
            Math::GMPz::Rmpz_primorial_ui($z, $p);
            $z;
        };

        if ($e > 1) {
            Math::GMPz::Rmpz_pow_ui($tmp, $prim, $e);
            Math::GMPz::Rmpz_mul($prod, $prod, $tmp);
        }
        else {
            Math::GMPz::Rmpz_mul($prod, $prod, $prim);
        }
    }

    return $prod;
}

# "HCN.bz2" was generated by Achim Flammenkamp, and is available at:
#   http://wwwhomes.uni-bielefeld.de/achim/HCN.bz2
# "HCN_primorial_deflated.txt.bz2" is a primorial deflation of each highly composite number H(n) from "HCN.bz2", such that A108951(A181815(H(n))) = H(n).
my $z = IO::Uncompress::Bunzip2->new("HCN_primorial_deflated.txt.bz2");

while (defined(my $line = $z->getline())) {

    chomp($line);
    my $hcn = primorial_inflation($line);

    if (primality_pretest($hcn - 1) and is_prob_prime($hcn - 1)) {
        print $., ", ";
    }
}

__END__
Indices n such that HCN(n) + 1 is prime (A306587):
1, 2, 3, 4, 5, 7, 9, 11, 12, 18, 20, 22, 23, 26, 28, 30, 34, 35, 44, 49, 54, 57, 60, 63, 74, 78, 84, 91, 97, 102, 104, 108, 111, 112, 114, 116, 118, 126, 134, 143, 145, 149, 159, 162, 167, 173, 177, 179, 188, 204, 208, 214, 230, 236, 238, 247, 249, 280, 294, 298, 299, 307, 311, 312, 318, 329, 330, 337, 350, 354, 356, 361, 362, 388, 415, 419, 427, 458, 462, 464, 478, 487, 518, 551, 555, 584, 603, 633, 640, 652, 667, 705, 713, 732, 751, 817, 826, 850, 852, 884, 889, 904, 914, 935, 944, 948, 952, 953, 957, 968, 987, 995, 1002, 1005, 1007, 1028, 1039, 1047, 1091, 1096, 1133, 1186, 1244, 1245, 1253, 1274, 1320, 1335, 1350, 1371, 1392, 1399, 1420, 1432, 1445, 1451, 1462, 1469, 1480, 1525, 1563, 1577, 1601, 1659, 1662, 1704, 1717, 1733, 1747, 1797, 1811, 1948, 1998, 2035, 2037, 2045, 2061, 2117, 2245, 2259, 2289, 2395, 2410, 2524, 2548, 2566, 2666, 2746, 2785, 2841, 2868, 2925, 2980, 2995, 3022, 3067, 3080, 3101, 3134, 3146, 3160, 3211, 3220, 3276, 3327, 3349, 3369, 3491, 3508, 3642, 3676, 3699, 3748, 3829, 3860, 3905, 3931, 3947, 4095, 4117, 4142, 4146, 4188, 4221, 4242, 4248, 4318, 4338, 4352, 4366, 4452, 4501, 4504, 4557, 4584, 4660, 4878, 4945, 4964, 4995, 5005, 5006, 5037, 5080, 5119, 5177, 5216, 5228, 5308, 5329, 5445, 5463, 5484, 5497, 5507, 5512, 5552, 5593, 5723, 5905, 5975, 5981, 6016, 6026, 6196, 6235, 6309, 6602, 6609, 6775, 6776, 6837, 6853, 6971, 7362, 7365, 7465, 7495, 7506, 7543, 7558, 7592, 7606, 7623, 7707, 7738, 7829, 7881, 7890, 7968, 8121, 8130, 8249, 8382, 8401, 8494, 8868, 8902, 9003, 9100, 9154, 9168, 9263, 9272, 9346, 9421, 9631, 9665, 9782, 9901, 9956, 10261, 10411, 10440, 10535, 10536, 10561, 10631, 10659, 10694, 10725, 10821, 11027, 11143, 11324, 11425, 11428, 11444, 11554, 11605, 11754, 11850, 11863, 12161, 12360, 12408, 12528, 12575, 12710, 12860, 13070, 13111, 13166, 13346, 13463, 13503, 13849, 13951, 14270, 14488, 14764, 14864, 14879, 14947, 15180, 15301, 15406, 15518, 15535, 15606, 15651, 16078, 16137, 16406, 16439, 16574, 16633, 16816, 16826, 16872, 17106, 17293, 17339, 17401, 17472, 17492, 17863, 17985, 18020, 18136, 18169, 18355, 18567, 18598, 18667, 18902, 18930, 19112, 19229, 19348, 19513, 19570, 19928, 20110, 20187, 20350, 20387, 20553, 20721, 20764, 20946, 21059, 21135, 21351, 21366, 21406, 21410, 21510, 21511, 21626, 21711, 21803, 21807, 22074, 22255, 22468, 22543, 22695, 22696, 22776, 22840, 22982, 23004, 23101, 23102, 23223, 23243, 23251, 23403, 23427, 23551, 23813, 24051, 24147, 24367, 24569, 24694, 24724, 24866, 24867, 24950, 24989, 25084, 25287, 25366, 25451, 25590, 25899, 26160, 26454, 26465, 26533, 26774, 26905, 27075, 27101, 27106, 27935, 28073, 28201, 28243, 28299, 28404, 28437, 29117, 29214, 29308, 29420, 29574, 29648, 29690, 30162, 30178, 30382, 30661, 30802, 31066, 31128, 31180, 31324, 31549, 31642, 31719, 31819, 31865, 31870, 31967, 32029, 32066, 32169, 32577, 32876, 32881, 32971, 33094, 33549, 33665, 33666, 33732, 33845, 34029, 34218, 34281, 34422, 34665, 34711, 34883, 35186, 35225, 35317, 35538, 35539, 35697, 35745, 36032, 36308, 36328, 37181, 37350, 37623, 37649, 37652, 37740, 37859, 38136, 38354, 38531, 38619, 38639, 38973, 39067, 39532, 40130, 40632, 40943, 41182, 41583, 42028, 42111, 42118, 42336, 42447, 42649, 42710, 42794, 42889, 42955, 43451, 44011, 44180, 44267, 44832, 44859, 45182, 45303, 45339, 45648, 45906, 46606, 47187, 47354, 47412, 48016, 48200, 48405, 48528, 48530, 48625, 48665, 49118, 49124, 49257, 49261, 49351, 49591, 49851, 49864, 50673, 50678, 50734, 52117, 52145, 52216, 52385, 52437, 52852, 52907, 53178, 53699, 53916, 54035, 54163, 54654, 54656, 55660, 55904, 56967, 57232, 57333, 57362

Indices n such that HCN(n) - 1 is prime (A306588):
3, 4, 5, 6, 8, 9, 11, 12, 13, 14, 15, 16, 19, 20, 21, 28, 30, 31, 37, 39, 40, 45, 52, 55, 65, 66, 67, 68, 70, 79, 81, 84, 101, 108, 118, 131, 132, 136, 143, 148, 149, 151, 163, 170, 174, 185, 191, 200, 203, 208, 212, 231, 259, 261, 286, 289, 297, 317, 326, 327, 328, 330, 336, 348, 353, 358, 359, 362, 366, 401, 429, 435, 444, 454, 465, 481, 494, 517, 522, 527, 545, 595, 614, 622, 645, 658, 668, 699, 737, 755, 787, 830, 838, 843, 845, 847, 868, 872, 917, 918, 938, 973, 990, 1002, 1004, 1016, 1075, 1088, 1105, 1125, 1155, 1168, 1177, 1197, 1211, 1232, 1236, 1260, 1294, 1306, 1312, 1334, 1341, 1373, 1394, 1407, 1473, 1506, 1515, 1518, 1539, 1540, 1544, 1547, 1565, 1580, 1689, 1705, 1742, 1767, 1779, 1789, 1823, 1833, 1835, 1855, 1893, 1897, 1911, 1963, 1990, 1993, 2015, 2038, 2039, 2070, 2075, 2095, 2116, 2130, 2156, 2267, 2329, 2335, 2395, 2396, 2433, 2454, 2455, 2466, 2499, 2500, 2507, 2510, 2561, 2570, 2582, 2624, 2634, 2644, 2647, 2667, 2730, 2738, 2754, 2765, 2768, 2772, 2787, 2800, 2808, 2905, 2917, 2937, 2947, 3023, 3084, 3111, 3137, 3160, 3176, 3227, 3237, 3252, 3314, 3384, 3426, 3438, 3466, 3499, 3547, 3593, 3627, 3720, 3733, 3796, 3812, 3929, 4004, 4015, 4020, 4060, 4072, 4077, 4127, 4252, 4260, 4313, 4375, 4389, 4408, 4454, 4479, 4595, 4688, 4700, 4718, 4736, 4789, 4800, 4828, 5060, 5091, 5245, 5377, 5431, 5442, 5469, 5500, 5510, 5517, 5645, 5649, 5893, 5989, 6171, 6181, 6192, 6225, 6241, 6433, 6530, 6545, 6552, 6642, 6651, 6691, 6729, 6755, 6808, 6931, 7056, 7100, 7102, 7107, 7136, 7162, 7269, 7516, 7621, 7661, 7686, 7711, 7717, 7773, 7798, 7864, 7915, 8010, 8034, 8258, 8322, 8324, 8338, 8416, 8644, 8678, 8725, 9033, 9103, 9132, 9198, 9247, 9287, 9288, 9429, 9456, 9487, 9558, 9787, 9792, 9815, 9922, 10049, 10535, 10569, 10949, 10958, 10971, 11197, 11342, 11460, 11568, 11647, 11753, 11898, 11947, 11972, 12118, 12173, 12174, 12236, 12277, 12495, 12517, 12629, 12756, 13167, 13202, 13309, 13368, 13374, 13383, 13431, 13705, 13839, 14051, 14161, 14188, 14195, 14267, 14367, 14464, 14525, 14541, 14859, 14893, 14930, 15258, 15316, 15502, 15637, 15676, 15713, 15777, 15823, 15897, 16009, 16088, 16272, 16302, 16519, 16591, 16617, 16646, 16819, 16829, 17187, 17287, 17297, 17597, 17741, 17955, 18013, 18061, 18137, 18364, 18492, 18799, 18971, 19020, 19102, 19123, 19439, 19575, 19814, 19817, 19841, 20029, 20136, 20251, 20270, 20427, 20444, 20467, 20550, 20632, 20741, 20788, 21122, 21302, 21334, 21415, 21478, 21591, 21931, 22554, 22620, 22646, 22841, 22846, 22850, 23198, 23379, 23384, 23393, 23778, 24119, 24278, 24538, 24764, 24916, 24946, 25204, 25253, 25673, 25747, 25857, 26056, 26111, 26321, 26851, 26870, 26909, 27065, 27300, 27695, 27756, 27813, 27836, 28003, 28372, 28449, 28798, 29044, 29224, 29346, 29424, 29551, 29745, 29747, 29914, 30011, 30350, 30357, 30654, 30984, 31109, 31158, 31268, 31284, 31320, 31359, 31879, 32082, 32198, 32362, 32377, 32381, 32407, 32805, 33500, 33600, 34243, 34867, 34874, 34904, 35465, 35615, 35795, 35897, 36967, 36995, 37040, 37067, 37120, 37286, 37315, 37560, 37703, 38258, 38291, 38442, 38465, 38469, 38686, 38985, 39012, 39191, 39579, 39839, 39906, 40978, 41393, 41567, 41951, 42055, 42104, 42234, 42623, 42701, 42712, 42921, 43124, 43720, 43841, 44167, 44481, 44487, 44907, 45000, 45078, 45445, 46188, 46267, 46292, 46396, 46550, 46649, 46658, 46926, 47210, 47395, 47989, 48042, 48060, 48234, 48534, 48546, 48560, 48849, 49308, 49336, 49500, 50074, 50499, 50575, 50888, 50953, 51163, 51296, 52056, 52401, 52630, 52988, 53057, 53173, 53254, 53300, 53533, 53669