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

#!/usr/bin/perl

# Daniel "Trizen" Șuteu
# Date: 27 June 2019
# https://github.com/trizen

# Find highly composite numbers H(n) such that H(n)-1 and H(n)+1 are both primes.
# In other words, find highly composite numbers sandwiched between twin primes.

# The sequence of indices of highly composite numbers A002182 which are between a twin prime pair, begins as:
#   3, 4, 5, 9, 11, 12, 20, 28, 30, 84, 108, 118, 143, 149, 208, 330, 362, 1002, 2395, 3160, 10535

# Example:
#   H(10535) = A108951(52900585920)
#            = 14742217487368791965347653720647452690286549052234444179664342042930370966727413549068727214664401976854238590421417268673037399536054005777393104248210539172848500736334237168727231561710827753972114334247396552090671649834020135652920430241738510495400044737265204738821393451152066370913670083496651044937158497896720493198891148968218874744806522767468280764179516341996273430700779982929787918221844760577694188288275419541410142336911631623319041967633591283303769044016192030492715535641753600000

# where H(10535)-1 and H(10535)+1 are both prime numbers.

# See also:
#   https://oeis.org/A321995 -- Indices of highly composite numbers A002182 which are between a twin prime pair.
#   https://oeis.org/A108951 -- Completely multiplicative with a(p) = p# for prime p, where x# is the primorial A034386(x).
#   https://oeis.org/A002182 -- Highly composite numbers, definition (1): where d(n), the number of divisors of n (A000005), increases to a record.

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 primality_pretest($hcn + 1)
        and is_prob_prime($hcn + 1)
        and is_prob_prime($hcn - 1)) {
        print $., ", ";
    }
}