[Python] Miller-Rabin Pseudoprime Test

The module file is here.

def modular_exponentiation( base, exponent, modulus ):
        """
        caluclating ( base ** exponent ) % modulus via the Right to Left Binary Algorithm.
        detail of "modular exponentiation", see follow URL.
        http://en.wikipedia.org/wiki/Modular_exponentiation
        """
        remainder = 1
        
        while( exponent > 0 ):
                if( (exponent & 1) > 0 ):
                        remainder = ( remainder * base ) % modulus
                
                exponent >>= 1
                base = ( base * base ) % modulus
                
        return remainder


def miller_rabin( almost_prime, base_list ):
        # base_list should contain some primes
        
        remainder= 1
        
        for base in base_list:
                exponent = almost_prime - 1
                while( exponent > 0 ):
                        # exponent were odd number and ( base ^ exponent ) ≡ 1 mod modulus,
                        # then continue the loop. Because almost_prime has a probability that it is a prime.
                        # The case when ( base ^ exponent ) ≡ -1 mod modulus means passing Miller-Rabin witness with the base then testing with next base.
                        # However ( base ^ exponent ) is not congruent, then immediately almost_prime is identified as pseudoprime: return False.
                        remainder = modular_exponentiation( base, exponent, almost_prime )
                        
                        if remainder == ( almost_prime -1 ):
                                break
                        if( remainder != 1 ) and ( remainder != (almost_prime -1 ) ):
                                return False
                        if( exponent % 2 ):
                                break
                        exponent  >>= 1

        prime = almost_prime
        return prime

Miller-Rabin witness
Define sequence r(n) = b ^ {(m-1)/2n} % m
When an integer m is a prime and also b ( b < m ) is a prime then:

a(1) ≡ 1		a(1) ≡ 1
a(2) ≡ 1		・
・	OR	・
・		a(n-1) ≡ 1
a(n) ≡ 1		a(n) ≡ -1

Miller-Rabin witness is application of this feature. Function modular_exponentiation() calculates value of a(n). Function miller_rabin() evaluates sequence a(n). If a(n) !≡ ( 1| -1 ) then m is not prime. Otherwise m has passed Miller-Rabin pseudoprime test with the base. However, with another bases, it might be failed. Therefore it is necessary to test with more than one primes.

Overview of Miller Rabin Pseudoprime test, see Wikipedia: Miller-Rabin primality test.