如何生成伪随机对合?

How to generate a pseudo-random involution?

要生成伪随机排列,可以使用 Knuth shuffles。对合是一种自逆排列,我想,我可以通过禁止多次触摸一个元素来调整洗牌。但是,我不确定我是否可以有效地做到这一点以及它是否会生成每个内卷 equiprobably.

恐怕需要举个例子:在一个集合{0,1,2}上,有6个排列,其中4个是对合。我正在寻找一种以相同概率随机生成其中一个的算法。

一个正确但非常低效的算法是:使用 Knuth 洗牌,如果没有对合则重试。

对合是一对一的映射,它是它自己的逆。任何密码都是一对一的映射;必须这样才能明确解密密文。

对于对合,您需要一个与其自身相反的密码。这样的密码是存在的,ROT13 就是一个例子。其他一些见 Reciprocal Cipher

对于你的问题,我建议使用 XOR 密码。选择一个随机密钥,至少与初始数据集中最长的数据一样长。如果您使用的是 32 位数字,则使用 32 位密钥。要置换,依次将密钥与每条数据进行异或。反向排列(相当于解密)是完全相同的异或运算,会得到原始数据。

这将解决数学问题,但绝对不是加密安全的。重复使用相同的密钥将允许攻击者发现密钥。我假设除了需要具有均匀分布的看似随机的内卷外,没有其他安全要求。

ETA:这是 Java 中我在第二条评论中谈论的内容的演示。作为 Java,我为您的 13 元素集使用索引 0..12。

public static void Demo() {

    final int key = 0b1001;

    System.out.println("key = " + key);
    System.out.println();

    for (int i = 0; i < 13; ++i) {

        System.out.print(i + " -> ");
        int ctext = i ^ key;

        while (ctext >= 13) {
            System.out.print(ctext + " -> ");
            ctext = ctext ^ key;
        }
        System.out.println(ctext);
    }

} // end Demo()

演示的输出是:

key = 9

0 -> 9
1 -> 8
2 -> 11
3 -> 10
4 -> 13 -> 4
5 -> 12
6 -> 15 -> 6
7 -> 14 -> 7
8 -> 1
9 -> 0
10 -> 3
11 -> 2
12 -> 5

如果转换后的键从数组的末尾掉下来,它会再次转换,直到它落在数组内。我不确定 while 构造是否符合函数的严格数学定义。

让我们在这里使用 a(n) 作为一组大小 n (as OEIS does) 的对合数。对于给定大小 n 的集合和该集合中的给定元素,该集合的对合总数为 a(n)。该元素必须在对合过程中保持不变或与另一个元素交换。使我们的元素保持不变的对合数是 a(n-1),因为这些是对其他元素的对合。因此,对合的均匀分布必须有 a(n-1)/a(n) 的概率保持该元素固定。如果要修复它,只需保留该元素即可。否则,选择我们的算法尚未检查的另一个元素与我们的元素交换。我们刚刚决定了集合中的一个或两个元素会发生什么:继续前进并一次决定一个或两个元素会发生什么。

为此,我们需要每个i <= n的对合计数列表,但这很容易用递归公式

完成

a(i) = a(i-1) + (i-1) * a(i-2)

(请注意,OEIS 的这个公式也来自我的算法:第一项计算保留第一个元素的对合数,第二项计算与之交换的元素。)如果你是使用对合,这可能很重要,足以分解为另一个函数,预先计算一些较小的值,并缓存函数的结果以提高速度,如以下代码所示:

# Counts of involutions (self-inverse permutations) for each size
_invo_cnts = [1, 1, 2, 4, 10, 26, 76, 232, 764, 2620, 9496, 35696, 140152]

def invo_count(n):
    """Return the number of involutions of size n and cache the result."""
    for i in range(len(_invo_cnts), n+1):
        _invo_cnts.append(_invo_cnts[i-1] + (i-1) * _invo_cnts[i-2])
    return _invo_cnts[n]

我们还需要一种方法来跟踪尚未决定的元素,以便我们可以有效地选择具有均匀概率的元素之一 and/or 将元素标记为已决定。我们可以将它们保存在一个收缩列表中,并在列表的当前末尾添加一个标记。当我们决定一个元素时,我们将当前元素移动到列表的末尾以替换决定的元素然后减少列表。有了这样的效率,这个算法的复杂度是 O(n),每个元素都有一个随机数计算,也许除了最后一个。没有比这更好的订单复杂度了。

这是 Python 3.5.2 中的代码。由于通过未定元素列表涉及的间接寻址,代码有些复杂。

from random import randrange

def randinvolution(n):
    """Return a random (uniform) involution of size n."""

    # Set up main variables:
    # -- the result so far as a list
    involution = list(range(n))
    # -- the list of indices of unseen (not yet decided) elements.
    #    unseen[0:cntunseen] are unseen/undecided elements, in any order.
    unseen = list(range(n))
    cntunseen = n

    # Make an involution, progressing one or two elements at a time
    while cntunseen > 1:  # if only one element remains, it must be fixed
        # Decide whether current element (index cntunseen-1) is fixed
        if randrange(invo_count(cntunseen)) < invo_count(cntunseen - 1):
            # Leave the current element as fixed and mark it as seen
            cntunseen -= 1
        else:
            # In involution, swap current element with another not yet seen
            idxother = randrange(cntunseen - 1)
            other = unseen[idxother]
            current = unseen[cntunseen - 1]
            involution[current], involution[other] = (
                involution[other], involution[current])
            # Mark both elements as seen by removing from start of unseen[]
            unseen[idxother] = unseen[cntunseen - 2]
            cntunseen -= 2

    return involution

我做了几个测试。这是我用来检查有效性和均匀分布的代码:

def isinvolution(p):
    """Flag if a permutation is an involution."""
    return all(p[p[i]] == i for i in range(len(p)))

# test the validity and uniformness of randinvolution()
n = 4
cnt = 10 ** 6
distr = {}
for j in range(cnt):
    inv = tuple(randinvolution(n))
    assert isinvolution(inv)
    distr[inv] = distr.get(inv, 0) + 1
print('In {} attempts, there were {} random involutions produced,'
    ' with the distribution...'.format(cnt, len(distr)))
for x in sorted(distr):
    print(x, str(distr[x]).rjust(2 + len(str(cnt))))

结果是

In 1000000 attempts, there were 10 random involutions produced, with the distribution...
(0, 1, 2, 3)     99874
(0, 1, 3, 2)    100239
(0, 2, 1, 3)    100118
(0, 3, 2, 1)     99192
(1, 0, 2, 3)     99919
(1, 0, 3, 2)    100304
(2, 1, 0, 3)    100098
(2, 3, 0, 1)    100211
(3, 1, 2, 0)    100091
(3, 2, 1, 0)     99954

对我来说这看起来很统一,我检查过的其他结果也是如此。