Here is a variation of that one that deals with log factorials, so it needs only regular arithmetic. It uses the gammln recently posted
here by
ewijaya.
It is extremely fast, solving 1..99 in about 2 seconds. (There may be roundoff error issues. I have not checked all the answers.)
Update: I added an independent confirmation with Math::Pari for every number I insert in the queue. They all come out correct. With checking added, the solution to 1..99 takes a minute and 23 seconds. The hardest to get is 92, which passes through 12985 factorial.
#!/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);
}