hcn_plus1_is_square.pl

#!/usr/bin/perl

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

# Find highly composite numbers H(n) of the form k^2 - 1.

# Known values of k:
#   5, 7, 11, 19, 29, 41, 71, 4159

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

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

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

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

    Math::GMPz::Rmpz_add_ui($hcn, $hcn, 1);

    say "Testing: $." if ($. % 10000 == 0);

    if (Math::GMPz::Rmpz_perfect_square_p($hcn)) {
        my $t = sqrtint($hcn + 1);
        say "Found: $t";

        if ($t > 4159) {
            die "New term: $t";
        }
    }
}