#!/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);
}