#!perl
use strict;
use warnings;
use integer;


my @primes = ( 2, 3, 5, 7 );

#
# Checks if a number is prime
#
sub isprime($;)
{
    my $number = shift or die;
    die "Not enough primes !" unless $primes[-1] * $primes[-1] >= $number;

    # should get more primes...

    foreach( @primes )
    {
        return 0 unless $number % $_;
    }
    return 1;
}

#
# Returns the prime number at the given index
#
sub prime($;)
{
    my $index = shift or die;
    my $nrprimes = scalar @primes;
    if( $index < $nrprimes )
    {
        # already computed
        die unless defined $primes[$index];
    }
    else
    {
        for( my $idx = $nrprimes-1; $idx < $index; ++$idx )
        {
            my $p = $primes[$idx];
            do
            {
                $p += 2;
            }until( isprime $p );
            push @primes, $p;
        }
    }

    return $primes[$index];
}

#
# returns a hash with prime factors as keys and powers as values
#
sub decompose($;)
{
    my $number = shift or die;
    die "Not enough primes !" unless $primes[-1] * $primes[-1] >= $number;

    # should get more primes instead...
    return ( 1 => 0 ) if $number == 1; # special case for 1
    my %factors;

    foreach my $p ( @primes )
    {
        until( $number % $p )
        {
            ++$factors{$p};
            $number /= $p;
        }
        last if $number == 1;
    }

    warn "$number still remaining" unless $number == 1;
    warn "No factors !" unless scalar keys %factors;
    %factors;
}

#
# Given two hashes returned by decompose above,
# add the powers of the corresponding prime factors
#
# addhash( decompose(x), decompose(y) ) is the same as
# decompose( x*y )
#
sub addhash(\%\%)
{
    my $left  = shift or die;
    my $right = shift or die;
    my %left  = %$left;
    my %right = %$right;

    while( my($k,$v) = each %right )
    {
        $left{$k} += $v; # a.k.a. $right{$k}
    }

    %left;
}

#
# Decompose a list by adding up the factors
# In effect computes decompose( product of list elements )
#
sub decompose_l(@;)
{
    my %r;
    foreach(@_)
    {
        my %h = decompose $_;
        my %s = addhash %r, %h;
        %r = %s;
    }
    %r;
}

#
# Checks whether a hash from decompose corresponds to a square
# (all powers must be even numbers).
#
sub square(%)
{
    my %h = @_;
    foreach( values %h )
    {
        return 0 if $_ & 1;
    }
    return 1;
}


#
# Returns a list of lists, each being a subset of ($low .. $high), of specified length

# Recursive implementation; a list of length $length is obtained
# by prepending $low to all sequences of $length-1,
# starting from $low+1, then $low+2, etc; all ending at $high
#
sub sequence($$$;);
sub sequence($$$;)
{
    my ($low,$high,$length) = @_;
    return undef if $length < 3; # only understands list of at least 3
    my @seq;
    if( 3 == $length )
    {
        foreach( $low+1 .. $high-1 )
        {
            push @seq, [ $low, $_, $high ];
        }
    }
    else
    {
        --$length;
        foreach( $low + 1 .. $high - $length + 1 )
        {
            my @x = sequence( $_, $high, $length);
            foreach( @x )
            {
                unshift @$_, $low;
            }
            push @seq, @x;
        }
    }
    @seq;
}

#
# Calculate Graham's function
#
my $verbose = 0;
sub Graham($;)
{
    my $number = shift or die;

    my %factors = decompose $number;

    my $number2 = $number * 2; # the presumed maximum of G(x)
    if( $number < 4 )
    {
        $number2 = 8;
    }
    my @seq = ( $number ); # The sequence starts with the number itself


    unless( square %factors ) # nothing to do if it's a square
    {
        GRAHAM:
        foreach my $gx ( $number+2 .. $number2 ) # candidate
        {
            print STDERR "$number: Checking $gx\nlength:" if $verbose;
            foreach my $len ( 3 .. $number2-$number ) # length of sequence

            {
                print STDERR "$len," if $verbose;
                my @s = sequence $number, $gx, $len; # get all sequences (ouch)

                foreach my $rseq ( @s )
                {
                    my %dec = decompose_l @$rseq;
                    if ( square %dec )
                    {
                        # The product of the sequence is a square.
                        # We found it !
                        @seq = @$rseq;
                        last GRAHAM;
                    }
                } # next sequence
            } # next length
            print STDERR "\n" if $verbose;
        } # next G(x) candidate
    }

    if( wantarray() )
    {
        @seq; # purely for nicety
    }
    else
    {
        $seq[-1];
    }
}

#####################################################################

# prime (pun intended) the pump
prime(100);

#$,=',';
#print @primes;

#my @grahams = <DATA>;
#chomp @grahams;
#print @grahams;
#print "\n";

my $v = 0;
#print map { $_ / ++$v } @grahams;

# use Win32;

#$v = 0;
# foreach( @grahams )
# {
#     ++$v;
# #    my $before = Win32::GetTickCount;
#     my $before = time;
#     my @graham = Graham $v;
#     my $after = time;
# #    my $after = Win32::GetTickCount;
#     printf( "G(%3d) = %3d ; %10d ms\n", $v, $graham[-1], $after-$before );
# 
#     die "*** " unless $graham[-1]==$_;
# }
# 
__END__
1
6
8
4
10
12
14
15
9
[snip more lines]


- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Running this on a machine with 1GHz Athlon processor and 256M RAM:
G(  1) =   1 ;          0 ms
G(  2) =   6 ;          0 ms
G(  3) =   8 ;          0 ms
G(  4) =   4 ;          0 ms
G(  5) =  10 ;         10 ms
G(  6) =  12 ;          0 ms
G(  7) =  14 ;         10 ms
G(  8) =  15 ;         20 ms
G(  9) =   9 ;          0 ms
G( 10) =  18 ;         40 ms
G( 11) =  22 ;        350 ms
G( 12) =  20 ;         30 ms
G( 13) =  26 ;       1763 ms
G( 14) =  21 ;         20 ms
G( 15) =  24 ;         80 ms
G( 16) =  16 ;          0 ms
G( 17) =  34 ;      38976 ms    39 seconds
G( 18) =  27 ;         80 ms
G( 19) =  38 ;     184085 ms     3 minutes !
G( 20) =  30 ;        180 ms
G( 21) =  28 ;         20 ms
G( 22) =  33 ;        401 ms

After 1 hour and 15 minutes I had to kill it...


--
Csaba Ráduly, Software Engineer, Sophos Anti-Virus
Email: csaba.raduly@sophos.com, Tel: 01235 559933, Web: www.sophos.com
Add live virus info to your website: http://www.sophos.com/link/vfeed



