but Math::Pari has very fast factorization algorithms.

Indeed it does--blindly fast. And that completely changes the perspective on using the factors in some sort of solution. I'm still finding my way around the huge set of functions in Math::Pari. If POD allowed tables, the looong list of functions could be made more compact and easier to navigate. Even so, I doubt I would have picked out sigma as being useful in this context from the description:

sigma(x,{k = 1})

sum of the k^{th} powers of the positive divisors 
of |x|. x must be of type integer.

The library syntax is sumdiv(x) ( = sigma(x)) 
or gsumdivk(x,k) ( = sigma(x,k)), 
where k is a C long integer.
[download]

Mmm'kay! :)

Do you have the time/inclination to post code for your solution?

Here is a simple version of the code. It just enumerates the divisors and binary-searches them:

#! /usr/bin/perl
use strict;
use Math::Pari qw(:int factorint sqrtint divisors);
for (@ARGV) {
  my $N = $_;
  my $a = factorint($N);
  my $sqrt = sqrtint($N);
  my $div = divisors($a);
  my $len = Math::Pari::matsize($div)->[1];
 
  my $low = 0;
  my $high = $len - 1;
  my $mid;
  while (1) {
     $mid = ($low + $high) >> 1;
     last unless $low < $high;
     if ($div->[$mid] > $sqrt) {
        $high = $mid - 1;
     } else {
        last if $low == $mid;
        $low = $mid;
     }
  }
  print "Result for $N is ",$div->[$mid],"\n";
}
[download]