Why not use xor-sum instead of normal sum? The probability is easier to bound this way. Assume that the hashes that you generate are h[1..n]. Suppose that your hashes retain $$$B$$$ bits of information and that they are randomly chosen from an uniform distribution: 0 <= h[i] < 2**B, P(h[i] = X | 0 <= X < 2**B) = 1 / 2**B.

Notice that this means that P(the jth bit of h[i] = 0) = P(the jth bit of h[i] = 1) = 1/2 for any 1 <= i <= n, 0 <= j < B.

Suppose that the array we are checking is v[1..n] and we know that 1 <= v[i] <= n.

The correct answer for a permutation would be h_ok = h[1] ^ h[2] ^ .. ^ h[n]. Suppose that v[1..n] is not a permutation and we want to calculate the probability that h[v[1]] ^ h[v[2]] ^ .. ^ h[v[n]] = h_ok.

By moving all terms on one side, we would be left with an equality like {0, 1} * h[1] ^ {0, 1} * h[2] ^ .. ^ {0, 1} * h[n] = 0. (as in h[1] either once or none ...)

For example, if v[1..7] = [3 4 5 3 3 2 7], we would have the equality h[3] ^ h[3] ^ h[1] ^ h[6] = 0 <=> h[1] ^ h[6] = 0.

We will prove that the number of terms that we have on the left hand side doesn't matter and the probability of equality is $$$1/2^B$$$.

Let's prove first that the probability of the LSB to be equal to $$$0$$$ is $$$1/2$$$ no matter how many terms we have on the LHS.

If we have $$$1 \leq k \leq n$$$ terms on the left hand side, then the result is $$$0$$$ if either $$$0, 2, 4, ... $$$ out of the $$$k$$$ bits are $$$1$$$:

$$$P(LHS = 0) = \frac{{k \choose 0} + {k \choose 2} + {k \choose 4} + ...}{2^k} = \frac{2^{k-1}}{2^k} = \frac{1}{2}$$$.

Now, since the hashes have been uniformly chosen in $$$[0 .. 2^B)$$$, the bits' values are independent, so the probability for the LHS to be $$$0$$$ taking in consideration all of the bits is $$$(1/2) ^ B = \frac{1}{2 ^ B}$$$.

So if you use $$$B = 64$$$, the collision probability is $$$< 10^{-19}$$$.

import functools
import random

n = 1000
B = 10
h = [0] + [random.randint(0, (1<<B) - 1) for _ in range(1, n+1)]
correct_hash = functools.reduce(lambda a, b: a ^ b, h[1:])

nt = 10**5
cntCollisions = 0
for t in range(nt):
    #P(we generate a v that actually is a permutation) = n!/n^n < 1e-432.
    if functools.reduce(lambda a, b: a ^ b, [h[random.randint(1, n)] for i in range(1, n+1)]) == correct_hash:
        cntCollisions += 1

print(f"({cntCollisions}, {nt}) P_collision = {cntCollisions / nt}, expected = {1 / (2 ** B)}.")

#B = 5: P_collision = 0.03056, expected = 0.03125.
#B = 10: P_collision = 0.00087, expected = 0.0009765625.