[Tutor] primality testing

[Danny Yoo]

> def prime(n):
> 	for in in range(2, sqrt(n)+1):
> 		if n % i == 0:
> 			return 'Composite'
>
> 	return 'Prime'


Small typo; Dana meant to use 'i', not 'in', for the statement,

    for i in range(2, sqrt(n) + 1):






There's an alternative way of checking for primality that's very cool;
it uses randomness and number theory!  There's some code here that shows
how to do it:

    http://www-mitpress.mit.edu/sicp/chapter1/node17.html

The code in the link above is written in a language called Scheme, but
it's actually not too bad to do the conversion from Scheme to Python.
For example, if we see something like this:


;;; Scheme code
(define (expmod base exp m)
  (cond ((= exp 0) 1)
        ((even? exp)
         (remainder (square (expmod base (/ exp 2) m))
                    m))
         (else
          (remainder (* base (expmod base (- exp 1) m))
                     m))))
;;;


This code in Scheme does something akin to:

    (base ** exp) % m

It does it in a cute way so that we can avoid doing huge products --- it
applies the modulo '%' function at the same time that it's calculating the
powers, to keep the numbers small and easy to calculate.


That expmod function can be translated into Python like this:

###
def expmod(base, exp, m):
    """Calcuates base**exp modulo m."""
    if exp == 0:
        return 1
    elif is_even(exp):
        return square(expmod(base, exp / 2, m)) % m
    else:
        return (base * (expmod(base, exp-1, m))) % m
###




For those who are interested, here's the complete translation:

###
"""Primality testing using Fermat's little theorem."""

import random


def is_fast_prime(n, times):
    """Perform the fermat primality test 'times' times on the number n."""
    for i in range(times):
        if not fermat_test(n): return False
    return True


def fermat_test(n):
    """Tests to see if n is possibly prime.  It can give false
    positives, although it never gives false negatives."""
    def try_it(a):
        return expmod(a, n, n) == a
    return try_it(1 + random.randrange(n-1))


def is_even(x):
    """Returns true if x is even."""
    return x % 2 == 0


def square(x):
    """Returns the square of x."""
    return x * x


def expmod(base, exp, m):
    """Calcuates base**exp modulo m."""
    if exp == 0:
        return 1
    elif is_even(exp):
        return square(expmod(base, exp / 2, m)) % m
    else:
        return (base * (expmod(base, exp-1, m))) % m
###



Here's more information on primality tests from MathWorld:

    http://mathworld.wolfram.com/PrimalityTest.html


Hope this helps!