#!/usr/bin/perl

use integer;

# Create some private permanent variables.
my @primes     = ( 2, 3, 5 );
my @multiples  = ( 4, 9, 25 );
my $last_prime = 5;

# @heap is a heap of indices into @primes and @multiples.  The
# indexed @multiples value is used to sort the heap.  The increment
# used for successive tests always skips even numbers and multiples
# of three, so we start with only 2, the index of the prime number 5,
# on the heap.
my @heap = ( 2 );

# $last_incr is how we skip even and triple numbers.
# We always increment by 2 or 4.  Since $last_prime starts
# with an odd number (5), that always keeps us odd, and
# since we alternate adding 2 and 4, we never hit the odd
# multiples of 3 within any group of 6.  For instance,
# the first numbers are:  5+2 -> 7, 7+4 -> 11
my $last_incr = 4;

sub goldbach {
    use integer;
    my ($n) = shift;
    my ($low, $high, $primes);
    ($primes, $high) = primes($n); # Shown earlier in chapter.
    $low = 0;
#   return 1 unless $n > 2e10;     # Already tested; don't bother.
#                                  # (But primes() will cause problems
#                                  # if you go far beyond this point.)
    return if $n % 2;              # Return if the number is odd.
    while( $low <= $high ) {
        my $total = $primes->[$low] + $primes->[$high];
        if ($total == $n) {
            return ($primes->[$low], $primes->[$high]);
        } elsif ($total < $n) {
            ++$low;
        } else {
            --$high;
        }
    }

    print "GOLDBACH CONJECTURE REFUTED: $n\n";
    print "\a" while 1;
}

print "992 is ", join(' + ', goldbach(992)), "\n";

# Binary search @primes for the index of the largest prime
# which is <= the argument.
sub _bin_prime {
    my ($n, $low, $high) = (shift, 0, $#primes);
    while ( $low < $high ) {
	my $mid = ($low+$high+1) / 2;
	if ( $prime[$mid] > $n ) {
	    $high = $mid -1
            } else {
                $low = $mid;
            }
    }
    return $low;
}

# Maintain @heap as a heap of indices into @primes and @multiples.
# @heap is ordered by the indexed value of @multiples.  The first
# element has just been incremented and might be out of order.
sub _heap_down {
    my $hi = 0;                     # heap index
    my $pi = $heap[0];              # prime index
    my $pv = $multiples[$pi];       # prime value
    my $ci;                         # child index
    while ( ($ci = $hi*2 + 1) < @heap ) {
	++$ci
	    if ($ci < $#primes)
                && ($multiples[$heap[$ci+1]] < $multiples[$heap[$ci]]);
	last unless $multiples[$heap[$ci]] < $pv;
	$heap[$hi] = $heap[$ci];
	$hi = $ci;
    }
    $heap[$hi] = $pi;
}

# Return all of the primes up to at least $n
#   and the index of the highest prime which is <= $n.
sub primes {
    my $n = shift;

    # Tiny numbers are neither prime nor composite.
    return (\@primes, undef) if $n < 2;

    # If we've previously generated enough primes, return
    # the current array and the right index.
    return (\@primes, _bin_prime($n) ) if $n < $last_prime;

    # Otherwise, we need to extend the list of primes.
    while ( $n > $last_prime ) {
	# Generate another prime, go up by 2 or 4, skipping
	# all even numbers and all odd multiples of 3 as well.
	$last_prime += ($last_incr = 6-$last_incr);
	my $heap;
	while ( $last_prime > $multiples[$heap=$heap[0]] ) {
	    # Advance all primes whose next multiple has been passed.
	    $multiples[$heap] += 2 * $primes[$heap];
	    _heap_down();
	}
	# Skip this candidate if it's a multiple.
	next if $last_prime == $multiples[$heap];
	# $last_prime is the next prime.
	push ( @primes, $last_prime );
	# The first multiple we'll check is its square
	# (anything smaller is also a multiple of some
	# smaller prime).
	push ( @multiples, $last_prime*$last_prime );
	# Add it to the heap.  It must be the larger than any
	# element already on the heap, so there needn't be a
	# _heap_up function.

	push ( @heap, $#primes );
    }

    # No more searching needed - we stopped when we got high enough.
    # The value to return is either the index of the last element
    # (if it matches the target, and then the target is prime) or
    # the index of the second last element (when the last
    # element is greater than the target, which happens when the
    # target is composite); so we determine whether to subtract
    # 1 from the index of the last element:
    return (\@primes, $#primes - ($n != $last_prime));
}