hcn.pl

#!/usr/bin/perl

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

# Print the n-th highly composite number, using data generated by Achim Flammenkamp.

# 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 ntheory qw(:all);
use POSIX qw(ULONG_MAX);
use IO::Uncompress::Bunzip2;
use experimental qw(signatures);

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

# "HCN.bz2" was generated by Achim Flammenkamp, and is available at:
#   http://wwwhomes.uni-bielefeld.de/achim/HCN.bz2
my $z = IO::Uncompress::Bunzip2->new("HCN.bz2");

my $from = shift(@ARGV) // die "usage: perl $0 [from] [to=from]\n";
my $to   = shift(@ARGV) // $from;

if ($from < 0) {
    die "error: $from must be a positive integer\n";
}

if ($to < $from) {
    die "error: $to must be >= $from\n";
}

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

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

    if ($. < $from) {
        next;
    }

    my @fields = split(' ', $line);
    my $omega  = shift(@fields);

    my @primes = (($omega > 0) ? @{primes(nth_prime($omega))} : ());
    my $prod   = Math::GMPz->new(1);

    while (@primes) {
        my $k = shift(@fields) // die "unexpected error: $line\n";

        my $e = 1;
        if ($k =~ /^(\d+)\^(\d+)\z/) {
            $k = $1;
            $e = $2;
        }

        for (1 .. $e) {
            my $p = shift(@primes);

            if ($k == 1) {
                Math::GMPz::Rmpz_mul_ui($prod, $prod, $p);
            }
            elsif ($p**$k < ULONG_MAX) {
                Math::GMPz::Rmpz_mul_ui($prod, $prod, powint($p, $k));
            }
            else {
                Math::GMPz::Rmpz_ui_pow_ui($tmp, $p, $k);
                Math::GMPz::Rmpz_mul($prod, $prod, $tmp);
            }
        }
    }

    if ($. >= $from) {
        say "$. $prod";
        last if $. >= $to;
    }
}