minimal factors for abundance

[hv at crypt.org]

(I tried submitting this through the webform, but I've not received my
copy of the email after a few days so I assume it failed.)

In searching for odd weird numbers, I found a need to calculate the
minimum number of prime factors required for a number to be abundant,
given that the least prime factor was some specified prime. This
gives rise to the following sequence.

%I A000001
%S A000001 3 5 9 18 31 46 67 91 122 157 194 238 284 334 392 455 522 591 668
%T A000001 748 834 929 1028 1133 1241 1352 1469 1594 1727 1869 2019 2163 2315
%U A000001 2471 2636 2802 2977 3156 3341 3534 3731 3933 4145 4358 4581 4811
%N A000001 a(n) is the least number of prime factors for any abundant number with p_n (the nth prime) as its least factor
%e A000001 a(2) = 5 since 945 = 3^3.5.7 is an abundant number with p_2 = 3 as its smallest prime factor, and no such number exists with fewer than 5 prime factors
%C A000001 If we replace "abundant" in the definition with "non-deficient", we get the same sequence with an initial 2 instead of 3, barring an astronomically unlikely coincidence with some as-yet-undiscovered odd perfect number
%Y A000001 A000040, A005101, A006038, A001222
%K A000001 nonn
%O A000001 1,1
%A A000001 hv at crypt.org (Hugo van der Sanden)

I include the program I used below - it's too long to include in the
database, but may be of interest.

Hugo
---
#!/usr/bin/perl -w
use strict;

=head1

Let p_k be the k'th prime, p_1 = 2. Then a(n) is the least number of prime
factors for an abundant number that has p_n as its smallest prime factor.

We need f(n) = sigma(n)/n > 2; if n = p^k.q with (p, q) = 1 then:
  f(pn) = f(n) * sigma(p^{k+1}) / (p.sigma(p^k))
        = f(n) (1  + 1/(sigma(p^{k+1}) - 1))
To minimise the number of prime factors, at each step we look to include
a permitted new prime that maximises this contribution, and therefore that
minimises sigma(p^{k+1}).

=cut

use Math::Pari qw/ PARI /;

my @p; BEGIN { @p = (PARI(2), PARI(3)) };
sub nextp {
    my $n = $p[$#p];
    N: while ($n += 2) {
        for (0 .. $#p) {
            next N if ($n % $p[$_]) == 0;
            last if $n < $p[$_] * $p[$_];
        }
        push @p, $n;
        return $n;
    }
}

$| = 1;
for (my $k = 0; 1; ++$k) {
    my @power = ($k);
    my $prod = PARI(2);
    my $sigma = PARI(1);
    my $count = 0;
    # ($sigma < $prod) finds non-deficient rather than abundant numbers
    while ($sigma <= $prod) {
        ++$count;
        my($min, $index) = ($p[$power[0]] + 1, 0);
        for (1 .. $#power) {
            my $this = contrib($p[$power[$_]], $_ + 1);
            next unless $this < $min;
            ($min, $index) = ($this, $_);
        }
        $prod *= $p[$power[$index]];
        $sigma *= $min;
        $sigma /= contrib($p[$power[$index]], $index) if $index > 0;
        ++$power[$index];
        nextp() unless $p[$power[$index]];
        push @power, $k if $index == $#power;
    }
    print "$count ";
}

sub contrib {
    my($p, $k) = @_;
    my $n = 1;
    $n = $n * $p + 1 for 1 .. $k;
    $n;
}