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