Perl Monk, Perl Meditation PerlMonks

### Re^2: Challenge: Chasing Knuth's Conjecture

by tall_man (Parson)
 on Mar 29, 2005 at 17:44 UTC ( #443204=note: print w/replies, xml ) Need Help??

in reply to Re: Challenge: Chasing Knuth's Conjecture
in thread Challenge: Chasing Knuth's Conjecture

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

Replies are listed 'Best First'.
Re^3: Challenge: Chasing Knuth's Conjecture
by hv (Parson) on Mar 29, 2005 at 19:10 UTC

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 => ssssssssfsssssssssssfsssssfssssssssssfssssssssssssfsssssssssssss
+fssssfsssssssssfsssssssssssfsssssssfssssssssfssssssssssfssssfsssssssf
+sssssssssssfssssssssfssssssssssfsssssssssssfsssssssfsssssssssfssssssf
+sssssssssfsssssssssfsssssssfssssssfsssfsssssssssfssssssfssssssfssssfs
+sssssssfsfssfssssssfsssssfssfsfssff

Hugo

Create A New User
Node Status?
node history
Node Type: note [id://443204]
help
Chatterbox?
and all is quiet...

How do I use this? | Other CB clients
Other Users?
Others browsing the Monastery: (8)
As of 2018-04-22 21:11 GMT
Sections?
Information?
Find Nodes?
Leftovers?
Voting Booth?
My travels bear the most uncanny semblance to ...

Results (83 votes). Check out past polls.

Notices?