Beefy Boxes and Bandwidth Generously Provided by pair Networks
Your skill will accomplish
what the force of many cannot
 
PerlMonks  

Re^2: Challenge: Chasing Knuth's Conjecture

by tall_man (Parson)
on Mar 29, 2005 at 17:44 UTC ( [id://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 (Prior) 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 => ssssssssfsssssssssssfsssssfssssssssssfss ssssssssssfsssssssssssssfssssfsssssssssf sssssssssssfsssssssfssssssssfssssssssssf ssssfsssssssfsssssssssssfssssssssfssssss ssssfsssssssssssfsssssssfsssssssssfsssss sfsssssssssfsssssssssfsssssssfssssssfsss fsssssssssfssssssfssssssfssssfssssssssfs fssfssssssfsssssfssfsfssff

    Hugo

    Update 2023-11-27: broke long line

Log In?
Username:
Password:

What's my password?
Create A New User
Domain Nodelet?
Node Status?
node history
Node Type: note [id://443204]
help
Chatterbox?
and the web crawler heard nothing...

How do I use this?Last hourOther CB clients
Other Users?
Others pondering the Monastery: (5)
As of 2024-03-28 23:16 GMT
Sections?
Information?
Find Nodes?
Leftovers?
    Voting Booth?

    No recent polls found