Beefy Boxes and Bandwidth Generously Provided by pair Networks
Perl-Sensitive Sunglasses
 
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); }


Comment on Re^2: Challenge: Chasing Knuth's Conjecture
Download Code
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

Log In?
Username:
Password:

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

How do I use this? | Other CB clients
Other Users?
Others having an uproarious good time at the Monastery: (12)
As of 2014-10-23 12:45 GMT
Sections?
Information?
Find Nodes?
Leftovers?
    Voting Booth?

    For retirement, I am banking on:










    Results (125 votes), past polls