Here is a variation of that one that deals with log factorials, so it needs only regular arithmetic. It uses the gammln recently posted
It is extremely fast, solving 1..99 in about 2 seconds. (There may be roundoff error issues. I have not checked all the answers.)
#!/usr/bin/perl
use strict;
use Inline qw(C);
use Math::Pari qw(ifact sqrtint PARI);
my $max = shift || 10;
my @try;
my $maxlog = log(1024*1024);
my $log2 = log(2.0);
my $biggest = 0;
# deal with '1' explicitly
my %seen = (1 => 's');
my $waiting = $max - 1;
print "1 => s\n";
insert(3, '');
my $m;
while (@try) {
my $n = shift @try;
$biggest = $n if $n > $biggest;
my $s = 'f' . $seen{$n};
#$n->bfac; # $n = ($n)!
my $logfact = gammln($n + 1.0);
if ($logfact > $maxlog) {
while ($logfact > $maxlog) {
$logfact /= 2.0;
$s = "s$s";
}
# truncate the square root.
$m = int(exp($logfact));
} else {
# round a straight gamma log.
$m = int(exp($logfact) + .5);
}
while (!defined $seen{$m}) {
insert($m, $s);
$m = int(sqrt($m));
$s = "s$s";
}
}
sub insert {
my($n, $s) = @_;
if ($s ne "") {
my $check_n = confirm($s);
if ($check_n != $n) {
print "$n did not confirm correctly -- expression gave $check_
+n\n";
}
}
if ($n <= $max) {
print "$n => $s\n";
if (--$waiting == 0) {
print "biggest tried was $biggest\n";
exit 0;
}
}
$seen{$n} = $s;
my($min, $max) = (0, scalar @try);
while ($min + 1 < $max) {
my $new = ($min + $max) >> 1;
if ($try[$new] > $n) {
$max = $new;
} else {
$min = $new;
}
}
++$min if $min < @try && $try[$min] < $n;
splice @try, $min, 0, $n;
}
sub confirm {
my $s = shift;
my @ops = reverse split //,$s;
my $val = PARI 3;
foreach my $op (@ops) {
if ($op eq "f") {
#print "doing $val factorial\n";
$val = ifact($val);
} else {
$val = sqrtint($val);
#print "doing sqrt\n";
}
}
#"final result is $val\n";
return $val;
}
__END__
__C__
double gammln(double xx)
{
double x,y, tmp,ser;
static double cof[6] = {76.18009172947146, -86.50532032941677,
24.01409824083091, -1.231739572450155,
0.12086509738661e-2, -0.5395239384953e-5};
int j;
y = x = xx;
tmp= x+5.5;
tmp -= (x+0.5)*log(tmp);
ser=1.000000000190015;
for (j=0;j<=5;j++) ser += cof[j]/++y;
return -tmp+log(2.5066282746310005*ser/x);
}