Re^4: Challenge: Chasing Knuth's Conjecture

[Roy Johnson]

My code looks like it agrees with yours; I set the factorial limit to 5000, and it ran for 20 or 30 minutes, I think. 64 is the first skip at that setting.

I did a simple translation of your Inline C sub into Perl, and it is still very fast. Running it for 64 took 33 seconds, and it popped out:

64 => sssssssssfssssssfsssssssssssssfsssssssssssfssssssssssssfssssssss
+ssfssssssf
ssssssssfsssssssssfsssssssssssfsssssssfssssssfsssssssssssfssssssssfsss
+ssssssssfs
sssssfssssssfssssssfssssssssssfssssssssfssssssssfssssssfssssssfssssfss
+ssssssfsfs
sfssssssfsssssfssfsfssff
biggest tried was 5071

Code follows.

#!/usr/bin/perl
use strict;
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;
}

sub gammln {
  my $xx = shift;
  my @cof = (76.18009172947146, -86.50532032941677,
             24.01409824083091, -1.231739572450155,
             0.12086509738661e-2, -0.5395239384953e-5);
  my $y = my $x = $xx;
  my $tmp = $x + 5.5;
  $tmp -= ($x + .5) * log($tmp);
  my $ser = 1.000000000190015;
  for my $j (0..5) {
     $ser += $cof[$j]/++$y;
  }
  -$tmp + log(2.5066282746310005*$ser/$x);
}