#!/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;

# 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));
}

# $val = prime($n) returns the $nth prime.
# prime(0) is 2, prime(1) is 3, prime(2) is 5, and so on.
sub prime {
    my $n = shift;
    # Make sure there are $n+1 primes already found.
    while ( $n >= @primes ) {
        # $n - $#primes is how many more we want, so the
        # highest we will find is at least twice that much
        # larger than the highest we found so far.
        primes( $last_prime + 2*($n - $#primes) );
    }
    return $primes[$n];
}


# $num = prime_encode( $n1, $n2, ... )
sub prime_encode {
    my $index  = 0;
    my $result = 1;
    while ( @_ ) {
        $result *= prime($index++) ** (shift);
    }
}

# ($n1,$n2, ...) = prime_decode( $num );
sub prime_decode {
    my ($num, $index, @result) = (shift, 0);

    while ( $num > 1 ) {
        use integer;
        my $prime = prime($index++);
        my $div = 0;
        until ( $num % $prime ) {
            $num /= $prime;
            ++$div;
        }
        push @result, $div;
        $num     -= $div * $prime;
    }
}

print prime(100), "\n";

# Check whether $n is prime, by trying up to 5 random tests.
sub is_prime {
    use integer;
    my $n  = shift;
    my $n1 = $n - 1;
    my $one = $n - $n1;     # 1, but ensure the right type of number.

    my $witness = $one * 100;

    # find the power of two for the top bit of $n1.
    my $p2 = $one;
    my $p2index = -1;
    ++$p2index, $p2 *= 2
        while $p2 <= $n1;
    $p2 /= 2;

    # number of iterations: 5 for 260-bit number, go up to
    # 25 for much smaller numbers.
    my $last_witness = 5;
    $last_witness += (260 - $p2index)/13 if $p2index < 260;

    for $witness_count ( 1..$last_witness ) {
        $witness *= 1024;
        $witness += rand(1024);
        $witness = $witness % $n if $witness > $n;
        $witness = $one * 100, redo if $witness == 0;

        my $prod = $one;
        my $n1bits = $n1;
        my $p2next = $p2;

        # compute $witness ** ($n - 1).
        while (1) {
            # Is $prod, the power so far, a square root of 1?
            # (plus or minus 1)
            my $rootone = $prod == 1 || $prod == $n1;

            $prod = ($prod * $prod) % $n;

            # An extra root of 1 disproves the primality.
            return 0 if $prod == 1 && !$rootone;

            if ( $n1bits >= $p2next ) {
                $prod = ($prod * $witness) % $n;
                $n1bits -= $p2next;
            }
            last if $p2next == 1;
            $p2next /= 2;
        }
        return 0 unless $prod == 1;
    }
    return 1;
}

foreach (qw(1089 17 32109481 55441)) {
    if (is_prime($_)) {
	print "$_ is prime.\n";
    } else {
	print "$_ is not prime.\n";
    }
}