Пусть у нас есть 

 брачных пар и 

 лодок. Будем считать, что лодки различимы (пронумерованы), а места в лодках нет. Тогда количество всех рассаживаний людей по лодкам (исходов) равно 

.
Подсчитаем количество благоприятных исходов, используя принцип включений-исключений. Пусть свойство 

, где 

 и 

, состоит в том, что в лодку номер 

 попадает брачная пара номер 

. При этом рассаживание будет благоприятным тогда и только тогда, когда оно не удовлетворяет ни одному из свойств 

. Пусть 

 - это множество всех свойств, и для любого подмножества 

 число рассаживаний удовлетворяющих свойствам из 

 равно 

. Тогда по принципу включений-исключений число рассаживаний, не удовлетворяющий никакому свойству из 

 равно:
 
Попробуем представить эту сумму в удобоваримом виде.
Прежде всего заметим, что если в 

 больше двух свойств вида 

 для фиксированного 

, то 

 (в лодку не могут сесть больше 2-х пар). Аналогично, если в 

 больше одного 

 для фиксированного 

, то 

 (каждая пара может сесть только в одну лодку). Для остальных 

 определим сигнатуру как двухкомпонентный вектор 

, где 

 и 

, то есть 

 - это число лодок куда согласно 

 должно сесть как минимум 

 брачных пар. Понятно, что если 

 имеет сигнатуру 

, то 

.
Вычислим 

 для подмножества 

 заданной сигнатуры 

. Понятно, что элементы 

 однозначно определяют 

 лодок и 

 брачных пар, которые в них сидят. Для оставшихся 

 пар остается 

 пустых лодок и 

 лодок, заполненных наполовину (т.е. с одной уже сидящей парой). Поэтому
 
Чтобы дать окончательный ответ, осталось подсчитать, а сколько у нас есть подмножеств 

 заданной сигнатуры 

. Так как 

, то чтобы определить 

, можно сначала выбрать брачные пары, которым 

 предписывает сидеть в лодках, потом выбрать 

 лодок, и наконец распределить выбранные пары по выбранным лодкам. Отсюда мы немедленно получаем, что количество различных 

 заданной сигнатуры 

 равно
 
Собирая все вместе, получаем количество благоприятных исходов:
 
Численные значения для 

:
Код:
? f(n)=sum(m1=0,n,sum(m2=0,n-m1,(-1)^m1*(2*n)!*n!*(2*n-m1-2*m2)!/(n-m1-m2)!/m1!/m2!/2^(2*n-2*m1-m2)))
? vector(10,n,f(n))
%1 = [0, 6, 900, 748440, 1559930400, 6928346502000, 58160619655538400, 845986566719614320000, 19957466912796971445888000, 724891264860942581350908960000]
Соответственно, для 

 получаем искомую вероятность:
