Version: 1

Type: Function

Category: Algorithms

License: GNU General Public License

Description: Two methods: is_prime - Test if a number is prime. get_prime - Find a prime within a given range. Uses the Miller-Rabin randomized primality test.



<?php
/*
 * Miller-Rabin primality tester...
 *
 * Methods:
 * is_prime( $number, $certainty )
 * Returns boolean indicating whether number is prime or not.
 * If it returns false, the number is _definitely_ composite.
 * If it returns true, the probability that it is wrong is 1 / 2^$certainty.
 *
 * get_prime( $min, $max, $certainty )
 * Returns a number between $min and $max that is probably prime, or false if none exist.
 * Again, the probability that a number returned is actually composite is 1/2^$certainty.
 *
 * (In general, a certainty of 10 will be enough - this is a 99% chance of correctness.)
 */

function get_prime( $min, $max, $certainty )
{

  if( $min >= $max )
    return( false );

  while( $min < $max )
    {
      if( is_prime( $min, $certainty ) )
	return( $min );
      else
	$min++;
    }
    
  return( false );

}

function is_prime( $n, $s )
{

  if( $n < 2 )
    return( false );
  
  if( $n == 2 )
    return( true );

  if( ( $n % 2 ) == 0 )
    return( false );

  srand( make_seed() );
  $randval = rand();

  for( $index = 0; $index < $s; $index++ )
    {

      $a = rand( 1, ( $n - 1 ) );
      if( mr_witness( $a, $n ) )
	return( false );
      
    }
  
  return( true );

}

function mr_witness( $a, $n )
{
  if( !is_int( $a ) || !is_int( $n ) )
    return( false );

  $bits = array_reverse( explode( ":", chunk_split( sprintf( "%b", ( $n - 1 ) ), 1, ":" ) ) );
  $d = 1;
  for( $index = count( $bits ); $index  >= 0; $index-- )
    {
      $x = $d;
      $d = ( $d * $d ) % $n;

 
      if( $d == 1 && $x != 1 && $x != ( $n - 1 ) )
	return( true );

      if( $bits[ $index ] == 1 )
	$d = ( $d * $a ) % $n;
    }
  
  if( $d != 1 )
    return( true );
  
  return( false );

}

function make_seed()
{

  list( $usec, $sec ) = explode( ' ', microtime() );
  return( (float)$sec + ( (float)$usec * 100000 ) );

}
?>