[Tutor] A small math puzzle [recreational Python]

[Kirby Urner]

A permutations with no element in its original position is
called a derangement.  Leonhard Euler came up with a formula
for computing the number of derangements, based on N elements,
as a function of N:

drs(n) = n! * [1/0! - 1/1! + 1/2! - 1/3! + ... + ((-1)^n)/n!]

(drs stands for derangements).

However, as Tim Peters pointed out on edu-sig awhile back, the
series in brackets is approaching 1/e (e = math.e) so for
larger values of n, you can just go drs(n) = round(n!/e) to
get a fairly good approximation.

My approach to the above function was to use long integers
for their precision, i.e. to not let floating point numbers
enter into it.  The code is simplified thanks to a fraction
object, which internally keeps numerator and denominator
separate, and which adds fractions by finding a common
denominator.

Many people have written such objects in Python.  Some
languages, like Scheme, have this rational number type built
in, and Python might too, one day -- the type system is
extensible after all.

Let's look at the guts of a simple Fraction object (good
enough to get the job done for this application -- even
a little bit fancier than we need):

    class F:

        def __init__(self,numer,denom):
           # reduce inputs to lowest terms
           gcd = self.gcd(numer,denom)
           self.num = numer//gcd
           self.den = denom//gcd

        def __add__(self,other):
            # find lowest common multiple
            a,b = self.num, other.num
            comden = self.lcm(self.den, other.den)
            if comden != self.den:
               a *= comden//self.den
            if comden != other.den:
               b *= comden//other.den
            return F(a+b,comden)

        def __sub__(self,other):
            return self + (-other)

        def __neg__(self):
            return F(-self.num,self.den)

        def gcd(self,a,b):
            # find greatest common divisor of a,b
            if b==0: return a
            else: return self.gcd(b,a%b)

        def lcm(self,a,b):
            # find lowest common multiple of a,b
            return a*b//self.gcd(a,b)

        def __repr__(self):
            # represent as (a/b)
            return "(%s/%s)" % (self.num, self.den)

  >>> m  = F(1,2)
  >>> n  = F(2,3)
  >>> n+m
  (7/6)
  >>> n-m
  (1/6)	

This isn't a complete Fraction definition, as we don't
have the ability to multiply fraction objects (easier than
__add__), or do some other things we'd like.  But we can
add and subtract, which is enough for computing:

[F(1,0!) - F(1,1!) + F(1,2!) - F(1,3!) + ... + F((-1)^n),n!)]

in the above example.

Except factorial notation isn't legal Python.  We need a
simple factorial function:

   >>> from operator import mul
   >>> def fact(n):
   	  if n<=1: return 1
	  else: return reduce(mul,range(2,n+1))

So now we're ready to implement the above expression in
brackets:

  >>> def brackets(n):
          sign = -1
          sum = F(1,1)
          for i in range(1,n+1):
              if sign == 1:
                  sum = sum + F(1,fact(i))
              if sign == -1:
                  sum = sum - F(1,fact(i))
              sign *= -1
          return sum

  >>> brackets(10)
  (16481/44800)
  >>> brackets(30)
  (3364864615063302680426807870189/9146650338351415815045120000000)
  >>> from __future__ import division
  >>> (3364864615063302680426807870189/9146650338351415815045120000000)
  0.36787944117144233
  >>> from math import e  # e used below as well
  >>> 1/e
  0.36787944117144233

So the brackets expression is indeed approaching 1/e, which is
a check on our algorithm.  The final step is to return the number
of derangements:

  >>> def drs(n):
          b = brackets(n)
          num = b.num
          den = b.den
          return fact(n)* num // den

  >>> drs(1)
  0
  >>> drs(2)
  1
  >>> drs(3)
  2
  >>> drs(4)
  9
  >>> drs(9)
  133496L
  >>> drs(30)
  97581073836835777732377428235481L  # too many to list out!

That's lots of derangements for 30 distinct element sets, but
is still less than the number of permutations, which is 30!:

  >>> fact(30)
  265252859812191058636308480000000L

Indeed, drs(30) is less than fact(30) by a factor of 1/e.
The rounding approach gives drs(30) to a lesser degree of
precision, owing to floating point's limitations:

  >>> round(fact(30)*1/e)
  9.758107383683579e+031

-- which is why I prefer the long integer approach, from a
number theoretic point of view, even though it's slower.

In a math-through-programming curriculum, I'd favor having
students evolve fraction objects, as per above, adding to
their capabilities over time, and then tossing out little
puzzles like this which make use of fraction object
capabilities.

Of course the drs(n) function could have be written without
any reference to fraction objects per se, but having fraction
objects available, like matrix objects, or polynomial objects,
helps keep the code clean.

Kirby