Probability that Secret Santa arrangement will result in perfect pairings So, we had Secret Santa at work.
We're 8 people. 
We each took turns and pulled a small piece of paper from a bowl with a name on it.
The only rule : If you pull your name, you have to put the piece of paper back in the bowl and try again.
Let's call the people A, B, C, D, E, F, G, H, which is also the order in which they picked their piece of paper.
We did the gift exchange last night.
A was F's secret santa.
B was E's secret santa.
C was D's secret santa.
D was C's secret santa.
E was B's secret santa.
F was A's secret santa.
G was H's secret santa.
H was G's secret santa.
See what happened? We made couples.  
A and F were each other's secret santa.
B and E were each other's secret santa.
C and D  were each other's secret santa.
G and H  were each other's secret santa.  
What is the probability of this happening and how do you calculate it?
 A: The total number of assignments among $2n$ people, where nobody is assigned to themselves, is $$d(2n) = (2n)!(1/2 - 1/6 + \cdots + (-1)^k/k! + \cdots + 1/(2n)!).$$  (These are called derangements.)  The value is very close to $(2n)! / e$.
If they correspond to perfect pairings, then they are a product of disjoint transpositions.  This implies their cycle structure is of the form 
$$(a_{11}a_{12})(a_{21}a_{22})\cdots(a_{n1}a_{n2}).$$
The number of distinct such patterns is the order of the group of all permutations of the $2n$ names divided by the order of the stabilizer of the pattern.  A stabilizing element may swap any number of the pairs and it may also permute the $n!$ pairs, whence there are $2^n n!$ stabilizing elements.  Therefore there are
$$p(2n) = \frac{(2n)!}{2^n n!}$$
such pairings.
Since all such perfect pairings are derangements, and all derangements are equally likely, the chance equals
$$\frac{p(2n)}{d(2n)} = \frac{1}{2^n n!(1 - 1/2 + 1/6 - \cdots + (-1)^k/k! + \cdots + 1/(2n)!)} \approx \frac{e}{2^n n!}.$$
For $2n=8$ people the exact answer therefore is $15/2119 \approx 0.00707881$ while the approximation is $e/(2^4\, 4!) \approx 0.00707886$: they agree to five significant figures.

To check, this R simulation draws a million random permutations of eight objects, retains only those that are derangements, and counts those that are perfect pairings.  It outputs its estimate, the standard error of the estimate, and a Z-score to compare it to the theoretical value.  Its output is
       p.hat           se            Z 
 0.006981031  0.000137385 -0.711721705

The small Z-score is consistent with the theoretical value.  (These results would be consistent with any theoretical value between $0.0066$ and $0.0073$.)
paired <- function(x) crossprod(x[x] - 1:length(x))==0
good <- function(x) sum(x==1:length(x)) == 0

n <- 8
set.seed(17)
x <- replicate(1e6, sample(1:n, n))
i.good <- apply(x, 2, good)
i.paired <- apply(x, 2, paired)

n.deranged <- sum(i.good)
k.paired <- sum(i.good & i.paired)
p.hat <- k.paired / n.deranged
se <- sqrt(p.hat * (1-p.hat) / n.deranged)
(c(p.hat=p.hat, se=se, Z=(p.hat - 15/2119)/se))

A: I was quite impressed by the elegance in @whuber answer. To be honest I had to do a lot of acquainting myself with new concepts to follow the steps in his solution. After spending a lot of time on it, I've decided to post what I got. So what follows is an exegetical note to his already accepted response. In this way there is no attempt at originality, and my only objective is to provide some additional anchoring points to follow some of the steps involved. 
So here it goes...
1. Why $2n$? Well this one may be way too basic: we need an even number of people to play the game.
2. Can we derive the formula for derangements? 
Following the Wikipedia entry with the example based on matching of people and hats, $n$ hats to be precise, we have that the options for the first person picking up a hat are shown here as follows:
$d(n)=(n-1)[d(n-2) + d(n-1)]=$
$=n\,d(n-2)-d(n-2)+n\,d(n-1)-d(n-1)$, which can be reorganized as:
$d(n) - n\,d(n-1) = -[d(n-1)-(n-1)\,d(n-2)]$. 
Now noticing the parallelism between the LHS of this equation and the part on the RHS within brackets we can continue recursively:
$d(n) - n\,d(n-1)=-[d(n-1)-(n-1)\,d(n-2)] =$ 
$=(-1)^2\,[d(n-2)-(n-2)\,d(n-3)]=\cdots=(-1)^{n-2}\,d(2)-2\,d(1)$
This implies that $d(n) = n\,d(n-1)+(-1)^n$.
Working backwards:
$d(2)=1$
$d(3)= 3\,d(2)-1 =3*1\,-1$
$d(4)= 4\,d(3)+1 =4*3*1\,-4\,\,+1$
$d(5)= 5\,d(4)-1 =5*4*3*1\,-5*4\,+5\,\,-1$
$d(6)= 6\,d(5)+1 = 6*5*4*3*1\,-6*5*4\,+6*5\,-6\,+1=$
$=6!\left(\frac{1}{2}- \frac{1}{3*2}+\frac{1}{4*3*2}-\frac{1}{5*4*3*2}+\frac{1}{6!}\right)=$
$=6!\left(\frac{1}{6!}-\frac{1}{5!}+\frac{1}{4!}-\frac{1}{3!}+\frac{1}{2!}-\frac{1}{1!}+1\right)$
So in general,
$d(n)=n!\left(1-1+\frac{1}{2!}-\frac{1}{3!}+\frac{1}{4!}+\cdots+\frac{1}{n!}\right)$
And remembering the Taylor series of $e^x$ evaluated at $x=-1$:
$d(n)\approx\frac{n!}{e}$
3. Disjoint transposition: The concept of transposition is easy to get from the link provided in the original answer, but "disjoint" was a bit less clear. Looking at an example, from the set ${a,b,c,d,e,f}$, the permutation ${b,d,a,c,f,e}$  can be expressed as a cycle as: a -> b -> d -> c after which it returns to a, but e -> f forms a loop onto itself - a disjoint cycle. Or, putting these two cycles together, the permutation can be expressed as the product $(\text {a b d c})(\text{e f})$.
In the Santa Claus question (eight employees happening to have drawn names in perfectly matched pairs: Anna gives Martha a gift, while Martha has drawn Anna's name) there will be $4$ closed loops.
4. To find the number of two-element loops we need to divide all possible permutations $(2n)!$ of the set of eight ($2n$) people by the total possible number of swaps of two elements $2^n$ and the total number of permutations of these pairs $n!$: $p(2n) = \frac{(2n)!}{2^n n!}$.

For the R simulation:
1. paired <- function(x) crossprod(x[x] - 1:length(x))==0
This function boils down to understanding x[x]: It is meant to evaluate an $8$ element vector representing the present assignments, and determine whether it is composed of 2-element loops, as in the Santa Claus problem. As long as the permutations correspond to transposing elements so that if Paul was supposed to give a present to Maria (Paul -> Maria and vice versa, Maria -> Paul) and Max to John (Max -> John / John -> Max) initially, the resultant transposition just results in new perfect pairing (Max -> Maria / Maria -> Max and Paul -> John / John -> Paul) we are fulfilling the initial condition of perfect pairing:

In other words, if we go back to the example of the hats in Wikipedia person i always takes back hat $1$. 
2. good <- function(x) sum(x==1:length(x)) == 0
This function evaluates whether we are dealing with a derangement by comparing the vector $x$ element wise to the the vector $(1,2,3,4,5,6,7,8)$, and making sure there is no coincidence.
3. k.paired <- sum(i.good & i.paired) is there to exclude paired permutations like the one above in the diagram, which are not derangements:
v <- c(1,2,3,4,5,6,7,8)
w <- c(1,2,3,5,4,6,7,8)

(c("is v paired?" = paired(v), "is w paired?" = paired(w),
   "is v a derang?" = good(w), "is w a derang?" = good(w)))

# not all paired permutations are derangements.

A: So this answer is prompted by a few questions in the comments. The earlier two answers are correct if the implicit assumption that all derangements are equally likely is true. This is sadly not the case with the method actually used in the question which is why quite a lot of secret santas are now using computer based methods as these can afford to discard the full drawing and re-start when the random assignment gives someone themselves.
To check that the difference actually was real, I implemented a set of python functions to systematically generate all possible options for a given set of people.
Firstly I implemented a version using the method that a computer will use which is generate a completely random option and repeat till we get a valid answer. This gives the same result as just counting the number of valid options (which is the same as the number of derangements) and assigning equal probability to each of them.
Second added in the handling of the replace and re-draw if you select your own name. As written a person could keep drawing them selves, putting their name back in and repeating for an infinite number of times, however we can show that this is actually the same as drawing any of the valid options with equal probability as
$$
\frac{1}{n}\sum_{i=0}^\infty\left(\frac{1}{n}\right)^i= \frac{1}{n-1}
$$
This gives us another way of viewing the problem. The probability of selecting each option now depends on if your name is still in the hat or not. We aren't suggesting people should check if they are in the hat, just that the probabilities are the same. The effect this difference in probabilities has is also compounded by the effect of the having to still abandon and re-start if we end up with the last name in the hat being the name of the last person to pick.
Sadly I haven't found an explicit formula for the probabilities of each option but I have implemented python code with this calculation and produced graphs of the resulting probabilities and calculated the probabilities for the pairing to happen.


The calculation done in the other answers gives a probability of 0.0070788107597923545 where as the accurate calculation for the stated method gives a probability of 0.006876714828601176. This isn't a large difference but it definitely is a measurable difference.
My code is included below for reference:
import pylab as pl

def generate_options(count):
    entries = set(range(count))
    return do_generate_options([], entries)

def do_generate_options(selected, available):
    """ Calculate the probabilities of a secret santa draw if we abandom if someone draws their own name. """
    if len(available) == 0:
        yield 1,tuple(selected)
        return
    factor = 1.0/len(available)
    for i in list(available):
        if i == len(selected):
            yield factor, None
        else:
            selected.append(i)
            available.remove(i)
            for p, option in do_generate_options(selected, available):
                yield p * factor, option
            selected.pop()
            available.add(i)

def do_generate_options(selected, available):
    """ Calculate the probabilities of a secret santa draw if you return the name and re-draw from that person if someone draws their own name. """
    if len(available) == 0:
        yield 1,tuple(selected)
        return
    factor = 1.0/len(available)
    if len(selected) in available:
        if len(available) == 1:
            yield 1, None
            return
        # sum 1/n^i for i from 1 to infinity = 1/n + 1/n * 1/(n-1) = 1/(n-1)
        factor = 1.0/(len(available) - 1)
    for i in list(available):
        if i != len(selected):
            selected.append(i)
            available.remove(i)
            for p, option in do_generate_options(selected, available):
                yield p * factor, option
            selected.pop()
            available.add(i)

def plot_probability_counts(count):
  data = {}
  for p, entry in generate_options(count):
    if entry is not None:
      if p in data:
        data[p] = data[p] + 1
      else:
        data[p] = 1
  entries = list(data.items())
  plt = pl.figure().add_subplot(111)
  plt.set_xlabel("probability")
  plt.set_ylabel("count")
  plt.plot([p * total for p, c in entries], [c for p, c in entries], "+")

def plot_probability_counts_scaled(count):
  data = {}
  total = 0
  for p, entry in generate_options(count):
    if entry is not None:
      total = total + 1
      if p in data:
        data[p] = data[p] + 1
      else:
        data[p] = 1
  entries = list(data.items())
  plt = pl.figure().add_subplot(111)
  plt.set_xlabel("probability/uniform probability")
  plt.set_ylabel("fraction")
  plt.plot([p * total for p, c in entries], [c / total for p, c in entries], "+")

def plot_probabilities_per_person(count):
  total = 0
  probabilites = [[0] * count for i in range(count)]
  for p, entry in generate_options(count):
    if entry is not None:
      total = total + p
      for i,j in enumerate(entry):
        probabilites[i][j] = probabilites[i][j] + p
  width = 1/(count + 1.5)
  xticks = range(count)
  plt = pl.figure().add_subplot(111)
  plt.set_xticks(xticks)
  plt.set_xticklabels(xticks)
  plt.set_xlabel("person giving")
  plt.set_ylabel("probability")
  plt.set_xlim(-0.5, count + 1.5)
  for i, a_probabilities in enumerate(probabilites):
    plt.bar([x + (i - (count - 1)/2) * width for x in xticks], [p/total for p in a_probabilities], width, label = i)
  plt.legend(title = "giving to", loc = "upper right")

def probability_paired(count):
  totalProbability = 0
  pairedProbability = 0
  totalCount = 0
  pairedCount = 0
  for p, entry in generate_options(count):
    if entry is not None:
      totalProbability = totalProbability + p
      totalCount = totalCount + 1
      valid = True
      for i in range(count):
        if entry[entry[i]] != i:
          valid = False
          break
      if valid:
        pairedProbability = pairedProbability + p
        pairedCount = pairedCount + 1
  return pairedProbability, totalProbability, pairedCount, totalCount

