hcn_divisible_by_2^n.pl

#!/usr/bin/perl

# Daniel "Trizen" Șuteu
# Date: 02 July 2019
# https://github.com/trizen

# Find highly composite numbers H(n) divisible by 2^k, for k={1,2,3,...}

# Known values of n:
#  2, 3, 6, 8, 21, 23, 50, 105, 143, 221, 634, 770, 1197, 2996, 5489, 7129, 18673, 25236, 46246, 96437, 179480, 298091, 425614

# See also:
#   https://oeis.org/A072847

use 5.020;
use strict;
use warnings;

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

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

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");

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

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

    if (Math::GMPz::Rmpz_scan1($hcn, 0) == $n) {
        print $., ", ";
        ++$n;
    }
}