Beefy Boxes and Bandwidth Generously Provided by pair Networks
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.

Log In?
Username:
Password:

What's my password?
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?
    Notices?