 No such thing as a small change PerlMonks

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

by Roy Johnson (Monsignor)
 on Mar 30, 2005 at 01:19 UTC ( #443318=note: print w/replies, xml ) Need Help??

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

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

Caution: Contents may have been coded under pressure.

Create A New User
Node Status?
node history
Node Type: note [id://443318]
help
Chatterbox?
and the web crawler heard nothing...

How do I use this? | Other CB clients
Other Users?
Others meditating upon the Monastery: (5)
As of 2020-01-17 14:23 GMT
Sections?
Information?
Find Nodes?
Leftovers?
Voting Booth?
The worst excuse I have ever heard is:

Results (84 votes). Check out past polls.

Notices?