Perl: Computation Howto

From AJS.COM

The Perl: Computation Howto is a collection of mathematical and algorithmic tips for Perl programmers.

Primes

Prime numbers are a frequent need in programming. Here are a few times you might need to deal with them:

Generating large primes

Here's a link to the source for a program that can generate large primes. It's pure Perl, so it's a bit slow: [1]

This code uses a two-pass approach. It first tests using quick_prime_test which uses a static list of primes to test against. In Perl, the fastest way to determin if a number is evenly divisible by another number is to use the modulus operator (%), thus the line:

return 0 if (($p % $sp) == 0);

Here $p is our potential prime and $sp is an entry from the list of static primes.

Then, we apply the Rabin-Miller test (subroutine rm_prime_test). It is impossible (at least so far) to determine if a large number is a prime quickly enough to be of use. However, the Rabin-Miller test gets us a very good guess that all but ensures that our number is prime. This is how prime numbers are generated for cryptographic applications. To apply this test, we go through steps relative to our potential prime, p:

• b = the number of 2s in the prime factors of p-1
  • We do this quickly by using the shift operator (>>) which essentially does an integer divide by 2.
• m = (p-1)/2**b
• a = a random integer which is less than p
• z = a**m mod p
• j = 0
• If z is 1 or z is p-1, then p passes
• Begin iterating (this is the loop labeled RM_ITERATION):

• If j is greater than 0 and z is 1, then p is not prime
• Increment j
• If j is less than b and z is not p-1 then:

• z = z**2 mod p
• go to the next iteration (RM_ITERATION)

• If z is p-1, then p passes
• Otherwise, the number is not prime

This whole test determines that a number is likely a prime, but can be repeated with new random values of a to be more certain. 5 repetitions of the test is considered to be a reasonable degree of certainty.

Time series analysis

Imagine that you have a time series "event 1 lasted from time T1 to time T2; event 2 lasted from time T3 to time T4; etc." Now imagine that you want to understand things about this time series like, "how many events were happening at time T5?"

Here's a simple bit of code that you might use:

my $active = 0; # no active events
my @ending;     # active events waiting to end
my $last = 0;   # Last known time
foreach my $event (sort {$a->{start} <=> $b->{start}} @events) {
  my @ended = grep { $_->{end} <= $event->{start} } @ending;
  @ending = grep { $_->{end} > $event->{start} } @ending;
  foreach my $end_event (sort {$a->{end} <=> $b->{end}} @ended) {
    active_between($last,$end_event->{end},$active);
    $last = $end_event->{end};
    $active--;
  }
  active_between($last,$event->{start},$active);
  $last = $event->{start};
  $active++;
}
foreach my $end_event (sort {$a->{end} <=> $b->{end}} @ended) {
  active_between($last,$end_event->{end},$active);
  $last = $end_event->{end};
  $active--;
}

The active_between function is just a subroutine that takes a start time, end time, and number of active events. Write this to do whateven you need to do with that information.