Abundant numbers

[hv at crypt.org]

franktaw at netscape.net wrote:
:Based on the recent discussion of "First even term", I have been inspired to submit the following sequence:
: 
:%S A000001 12,945,20,945,12,5391411025,12,945,20,81081,12,5391411025,12,6435,56,945,12,
:5391411025,12,81081,20,945,12,5391411025,12,945,20,6435,12,
:20169691981106018776756331,12,945,20,945,12,5391411025,12,945,20,81081,12,
:366005822969340125,12,945,56,945,12,5391411025,12,81081
:%N A000001 Smallest abundant number relatively prime to n.
[...]
:I would appreciate it if someone would check these before I actually submit
:the sequence.

I can confirm the figures; the code below takes me 22s to calculate
up to a(210) = 49061132957714428902152118459264865645885092682687973.

Hugo
---
#!/usr/bin/perl -w
use strict;
use Math::Pari qw/ PARI prime isprime factor /;

$| = 1;
my $b0 = PARI 0;
my $b1 = PARI 1;
my(@p, %avoid, @map, @factor, @fsigma);
for my $n (1 .. 210) {
    %avoid = map +($_ => 1), @{ factor($n)->[0] };
    print least($b1, $b1, $b0, 0), ",";
}

sub least {
    my($prod, $sigma, $best, $last) = @_;
    return $prod if $prod + $prod < $sigma;
    {
        my $next = @map ? $map[$#map] + 1 : 0;
        ++$next while $avoid{$p[$next] ||= prime($#p + 1)};
        my $p = $p[$next];
        next if $best && $prod * $p >= $best;
        push @map, $next;
        push @factor, 1;
        push @fsigma, $p + $b1;
        $best = least($prod * $p, $sigma * ($p + $b1), $best, $#factor);
        pop @map;
        pop @factor;
        pop @fsigma;
    }
    for (reverse $last .. $#factor) {
        local $factor[$_] = $factor[$_] + 1;
        next if $_ && $factor[$_] > $factor[$_ - 1];
        my $p = $p[$map[$_]];
        next if $best && $prod * $p >= $best;
        my $s2 = $sigma / $fsigma[$_];
        local $fsigma[$_] = $fsigma[$_] * $p + $b1;
        $best = least($prod * $p, $s2 * $fsigma[$_], $best, $_);
    }
    return $best;
}
__END__