How do I optimize this XOR sum algorithm?

The solution uses two tricks.

The first trick is to compute the n+1 prefix sums σ[0], …, σ[n] (as tobias_k notes), but with XOR instead of + (e.g., σ[2] = a[0] XOR a[1] XOR a[2]). The sublist from i to j has XOR-sum equal to σ[i−1] XOR σ[j]. If we loop i−1 and j over all possible pairs of values 0, …, n without considering the constraint i ≤ j, we get each sublist twice (once “forward”, once “backward”) and also n+1 extra zeros from when i−1 = j.

The second trick is that a fast Walsh–Hadamard transform can solve the following problem in O(n log n) time: given a list X and a list Y, we want the frequency counts of x XOR y for (x, y) ∈ X × Y. (For this problem, X = Y, but the structure of this trick is clearer if I use separate variables.) Why should we suspect that there is a fast algorithm in the first place? Limits aside, if it were x + y instead of x XOR y, then we would be looking to multiply polynomials quickly.

Let’s replace the list X with its frequency vector f, and Y with its frequency vector g. For example, X = [0, 0, 0, 2, 2, 3] becomes f = [3, 0, 2, 1]. Assuming that f and g have four elements, the desired result is

[f[0] g[0] + f[1] g[1] + f[2] g[2] + f[3] g[3]
,f[0] g[1] + f[1] g[0] + f[2] g[3] + f[3] g[2]
,f[0] g[2] + f[1] g[3] + f[2] g[0] + f[3] g[1]
,f[0] g[3] + f[1] g[2] + f[2] g[1] + f[3] g[0]
].

This is an example of an algebraic object called a symmetric bilinear form, which means that there exists some basis-changing matrix B such that the desired result is B⁻¹ (B f * B g), where * denotes element-wise multiplication. (Spoiler: B is a Walsh matrix.)

For intuition about the fast Walsh–Hadamard transform, which for every vector v can compute B v efficiently, let me show you what happens when I add the first two elements of the result:

f[0] g[0] + f[1] g[1] + f[2] g[2] + f[3] g[3]
+ f[0] g[1] + f[1] g[0] + f[2] g[3] + f[3] g[2]
= f[0] (g[0] + g[1]) + f[1] (g[1] + g[0]) + f[2] (g[2] + g[3]) + f[3] (g[3] + g[2])
= (f[0] + f[1]) (g[0] + g[1]) + (f[2] + f[3]) (g[2] + g[3])

and add the second two elements:

f[0] g[2] + f[1] g[3] + f[2] g[0] + f[3] g[1]
+ f[0] g[3] + f[1] g[2] + f[2] g[1] + f[3] g[0]
= f[0] (g[2] + g[3]) + f[1] (g[3] + g[2]) + f[2] (g[0] + g[1]) + f[3] (g[1] + g[0])
= (f[0] + f[1]) (g[2] + g[3]) + (f[2] + f[3]) (g[0] + g[1])

and subtract the first two elements:

f[0] g[0] + f[1] g[1] + f[2] g[2] + f[3] g[3]
− (f[0] g[1] + f[1] g[0] + f[2] g[3] + f[3] g[2])
= f[0] (g[0] − g[1]) + f[1] (g[1] − g[0]) + f[2] (g[2] − g[3]) + f[3] (g[3] − g[2])
= (f[0] − f[1]) (g[0] − g[1]) + (f[2] − f[3]) (g[2] − g[3])

and subtract the second two elements:

f[0] g[2] + f[1] g[3] + f[2] g[0] + f[3] g[1]
− (f[0] g[3] + f[1] g[2] + f[2] g[1] + f[3] g[0])
= f[0] (g[2] − g[3]) + f[1] (g[3] − g[2]) + f[2] (g[0] − g[1]) + f[3] (g[1] − g[0])
= (f[0] − f[1]) (g[2] − g[3]) + (f[2] − f[3]) (g[0] − g[1]) .

If we let f′ = [f[0] + f[1], f[2] + f[3]] and g′ = [g[0] + g[1], g[2] + g[3]], then the first two quantities are [f′[0] g′[0] + g′[1] g′[1], f′[0] g′[1] + f′[1] g′[0]], which is the same problem but half as large. Ditto the second two quantities. Finally, we can recover

x = f[0] g[0] + f[1] g[1] + f[2] g[2] + f[3] g[3]
y = f[0] g[1] + f[1] g[0] + f[2] g[3] + f[3] g[2]

from x + y and x − y as x = ((x + y) + (x − y))/2 and y = ((x + y) − (x − y))/2 (and likewise for all of the other pairs). (Note that Kache’s code defers the division to the end so that it can reuse the same transform.)

Note that your algorithms is not O(n²) but O(n³), but you can reduce it to O(n²) by just xor-ing each new number to all the results from the last "slice" (plus that number itself, starting a new subsequence).

from collections import Counter
def xor_sub_sum(arr):
    freq = Counter()
    last = []
    for x in arr:
        last = [x, *(x^y for y in last)]
        freq.update(last)
    # "most_common" does not consider smallest-key constraint...
    return max(freq.items(), key=lambda t: (t[1], -t[0]))

On my machine, this reduces the execution time for 1000 elements from 21.7 seconds to just 0.05 seconds. For 10000, it still takes ~5 seconds, though.

Okay, I still don't fully understand it, but I was able to pass all the tests by implementing an O(n log(n)) fwht() based on Wikipedia's description, plus pulling from another existing solution:

from collections import Counter
from itertools import accumulate
from operator import xor


def fwht(arr):
    # https://en.wikipedia.org/wiki/Fast_Walsh%E2%80%93Hadamard_transform
    if len(arr) == 1:
        return arr

    prefix, suffix = arr[:len(arr) // 2], arr[len(arr) // 2:]

    new_prefix = fwht([p + s for p, s in zip(prefix, suffix)])
    new_suffix = fwht([p - s for p, s in zip(prefix, suffix)])
    return new_prefix + new_suffix


def xorSubsequence(seq):
    next_pow2 = 2**(len(seq) - 1).bit_length()

    histogram = Counter(accumulate([0] + seq, xor))
    histogram = [histogram[value] for value in range(next_pow2)]

    histogram = [x * x for x in fwht(histogram)]
    histogram = [y // next_pow2 for y in fwht(histogram)]

    histogram[0] -= len(seq) + 1             # self combos (diagonal in table)
    histogram = [y // 2 for y in histogram]  # don't count things twice

    max_freq = max(histogram)
    return next((i, freq) for i, freq in enumerate(histogram) if freq == max_freq)

Why Xor sum is Dyadic convolution

Denote the input array as a.

Construct an array b, such that b[i]=a[0]⊕a[1]⊕...⊕a[i]. One can then construct a list M, M[i] stands for the number of element in b which has a value i. Note that some zero-padding is added to make the length of M be a power of 2.

Then consider the Dyadic (XOR) convolution. The difinition is(picture taken form this question)：

Consider conduct this Dyadic convolution between M and M, i.e. N=M*M where * stands for Dyadic convolution. Then N[i] is the sum of M[j]M[k] over all (j,k) that j⊕k=i.

Consider each subsequence xor(a[p:q]), we have xor(a[p:q])=b[p]⊕b[q]. For Every integer i, all the consecutive subsequences the xor results of which are equal to i can be converted to in this form(i=xor(a[p:q])=b[p]⊕b[q]). We further group this family of subsequences by the value of b[p] and the value b[q], for example, if xor(a[p1:q1])=xor(a[p2,q2])=i, and if b[p1]=b[p2],b[q1]=b[q2], these two subsequences will be grouped in the same sub-group. Consider sub-group (j,k), where the subsequences can be represented in the formi=xor(a[p':q'])=b[p']⊕b[q'], b[p']=j, b[q']=k, the number of subsquences in this sub-group (j,k) is M[j]M[k] (Recall that M[i] stands for the number of element in b which has a value i). So N[i] is the number sequence a[p:q] that xor(a[p:q])=i.

Howvever, since a[p:q] and a[q:p] are identical, we count every subsequence twice. So N[i] is twice of "number of consecutive subsequences that xor get i".

Use FWHT to compute the convolution

So now what we need is compute N=M*M, according to The Dyadic (XOR) convolution theorem(see proof here), we can first perform compute H(N)=H(M)×H(M). As H is involutive(see wiki), to get N just apply H on H(N) again.

Analysis of code

In this section, I will analysis the code provided by @Kache.

Actually, b is accumulate([0] + seq, xor). Using histogram = Counter(accumulate([0] + seq, xor)), one can get a dict of {possible_value_in_b: num_of_occurrence}. Then the next line, histogram = [histogram[value] for value in range(next_pow2)], this is the M mentioned above with padding added.

Then in histogram = [x * x for x in fwht(histogram)], so now histogram is H(N). And histogram = [y // next_pow2 for y in fwht(histogram)] serves as inverse transformation.

this is what histogram = [y // next_pow2 for y in fwht(histogram)] do.

The histogram[0] -= len(seq) + 1 eliminate the influence of the fact that a[p:p]=0. And histogram = [y // 2 for y in histogram] avoid counting twice (as stated before, N[i] counts a[p:q] and a[q:p] separately).