Re^2: Challenge: Chasing Knuth's Conjecture

[tall_man]

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);
}

[hv]

Interesting, I would have expected rounding errors large enough to skew the results, but your results agree with mine as far as the final value I've discerned so far in 1..99:

87 => 
ssssssssfsssssssssssfsssssfssssssssssfss
ssssssssssfsssssssssssssfssssfsssssssssf
sssssssssssfsssssssfssssssssfssssssssssf
ssssfsssssssfsssssssssssfssssssssfssssss
ssssfsssssssssssfsssssssfsssssssssfsssss
sfsssssssssfsssssssssfsssssssfssssssfsss
fsssssssssfssssssfssssssfssssfssssssssfs
fssfssssssfsssssfssfsfssff

Hugo

Update 2023-11-27: broke long line