http://www.enlight.ru/crypto/algorithms/rijndael/rijndael02.htm
http://csrc.nist.gov/CryptoToolkit/aes/rijndael/Rijndael-ammended.pdf

1. Invertibility;
2. Minimisation of the largest non-trivial correlation between linear combinations of
input bits and linear combination of output bits;
3. Minimisation of the largest non-trivial value in the EXOR table;
4. Complexity of its algebraic expression in GF(2^8);
5. Simplicity of description.

Там так и написано, что можно использовать любую таблицу замен с выполненными условиями (5 - необязательно

Нет нужды _явно_ строить все множества таких таблиц, чтобы выбрать из них "самую запутывающую", достаточно чётко сформулировать, чего хочется, подумать и построить для его генерации простой алгоритм.

Rijndael, как алгеброид, поанализироал немножко и предложил для генерации этого замуту с неприводимым многочленом
b(x)=(x^7+x^6+x^2+x)+a(x)(x^7+x^6+x^5+x^4+1) mod x^8 +1
(в русском реферате здесь ошибка - тривиальный многочлен написали , потому что математикам так свойственно думать.

Программеру проще - возьмём генератор случайных чисел и построим случайный цикл по всем 256 байтам (то есть такое заполнение, при котором операцией j=a[j] за 256 итераций обойдём всё), проверив по ходу построения, что не наделали ’opposite fixed points' (я, правда, не понял пока глубокого смысла этого противопоказания) - очевидно, что при идеальном генераторе этот алгоритм будет выдавать одинаково лучшие заполнения:
1. очевидно, т.к. заполнение взаимно однозначно по построению
2 и 3 -  эти критерии выполняются, т.к. 1-цикличное заполнение в этом плане точно лучше 2-и-более цикличных
4 - с непременностью: с большой долей вероятности, кратчайшее выражение будет составлять сумму 256 членов

Можно "улучшить" генератор, проверяя для каждого байта из сгенерированного заполнения все "correlation between linear combinations of input bits and linear combination of output bits". Как делать?
1. Перебираем все линейные комбинации бит в исходном байте - 255; во вложенном цикле - все линейные комбинации бит в замененном байте.
2. Внутри этих циклов проверим независимость линейной комбинации замены, скажем 0-го и 6-го бита выходного байта [06] от, скажем, комбинации 2-го, 3-го и 7-го бита исходного байта [237]. Это выглядит так: 256 байт заполнения произведут 256 отображений [237] -> [06], (например, 000->11, 000->01, 010->01, …), то есть 8 групп (с разными значениями [237]) по 32 отображения в [06] в каждой группе. [06] может иметь всего 4 значения – 00, 01, 10, 11
3. Нужно убедиться, что в каждой группе эти значения «размазаны» равномерно, то есть в группе, где исходные [237]=010 должно быть примерно 8 отображений [06]=00, 8 – 01, и по столько же 10 и 11. Допустим, отобразилось не по 8, а по 6, 8, 11, 7. Посчитаем какое-нибудь усреднение отклонения этих [6,8,11,7] от идеального [8,8,8,8] (тут нужно подумать, какое), сложим по другим группам и найдём максимум этой суммы по комбинациям линейных комбинаций (циклов п.1). Этот максимум можно использовать как соответствие таблицы замены критерию 2.

Но с хорошим генератором случайных чисел все сгенерированные подстановки должны быть «удачными» и перебором 10^70 вариантов врядли можно найти существенно более удачную…

По этим подстановкам, вносящим «нелинейность» в алгоритм шифрования много пишут, и проблема интереснее математикой, чем переборным поиском.