Beefy Boxes and Bandwidth Generously Provided by pair Networks
Perl Monk, Perl Meditation

Re: Challenge: Chasing Knuth's Conjecture

by hv (Parson)
on Mar 29, 2005 at 15:07 UTC ( #443148=note: print w/replies, xml ) Need Help??

in reply to Challenge: Chasing Knuth's Conjecture

The code below finds 1..10 in 2 seconds on this machine, having calculated nothing greater than 201!; 1..37 takes 8s, and calculates 366!. Beyond that we start to find some more difficult, but at 37 mins and 24MB it has so far found all but 8 numbers (59, 64, 72, 86, 87, 92, 97, 99) out of 1..99.

Update: still missing (92, 97, 99) after 274 mins (process size 67MB).

The basic ideas are: a) keep a sorted list (@try) of the numbers we've reached but not yet calculated a factorial for; b) at each iteration, take the smallest untried number and find it's factorial, and repeatedly take the square root of the result until we get to a number we've seen before; c) use a binary chop to insert new pending numbers.

#!/usr/bin/perl use strict; use bigint; # optionally with eg C< lib => 'Pari' > my $max = shift || 10; my @try; # deal with '1' explicitly my %seen = (1 => 's'); my $waiting = $max - 1; print "1 => s\n"; insert(3, ''); while (@try) { my $n = shift @try; my $s = 'f' . $seen{$n}; $n->bfac; # $n = ($n)! $s = "f$s"; while (!defined $seen{$n}) { insert($n, $s); $n = sqrt($n); $s = "s$s"; } } sub insert { my($n, $s) = @_; if ($n <= $max) { print "$n => $s\n"; exit 0 unless --$waiting; } $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; }

In the output, eg "5 => ssff" means that 5 = int(sqrt(int(sqrt(fact(fact(3)))))).


Replies are listed 'Best First'.
Re^2: Challenge: Chasing Knuth's Conjecture
by tall_man (Parson) on Mar 29, 2005 at 17:44 UTC
    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.

      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


Re^2: Challenge: Chasing Knuth's Conjecture
by Roy Johnson (Monsignor) on Mar 29, 2005 at 17:55 UTC
    I took your ideas for avoiding re-testing things and revamped my sieve. This generates all the answers that don't require factorial on anything higher than whatever you set (350 gets me all solutions for 1..30). It runs in about 5 seconds on my 1 GHz PC with that setting. Update: pregenerated factorials to save time; also put in a limit on how high to print, so the huge numbers don't show up. Even taking factorials as high as 999 (takes 30-40 seconds), there's no solution for 38.

    Caution: Contents may have been coded under pressure.
      I computed an answer for 38 with my modified version of hv's program, using logs. The result is:
      "sssssssssfsssssssssssfssssssfssssssfssssssfssssssssssfssssssssfssssss +ssfssssssfssssssfssssfssssssssfsfssfssssssfsssssfssfsfssff"

      I confirmed this answer with Math::Pari, and it is correct. The highest factorial required is 1861.

        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.

        Caution: Contents may have been coded under pressure.

Log In?

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

How do I use this? | Other CB clients
Other Users?
Others lurking in the Monastery: (4)
As of 2019-11-19 05:04 GMT
Find Nodes?
    Voting Booth?
    Strict and warnings: which comes first?

    Results (94 votes). Check out past polls.