How do I optimize this XOR sum algorithm?

Question

I'm trying to solve this hackerrank problem https://www.hackerrank.com/challenges/xor-subsequence/problem

from functools import reduce

def xor_sum(arr):
    return reduce(lambda x,y: x^y, arr)    

def xorSubsequence(arr):
    freq = {}
    max_c = float("-inf") # init val
    min_n = float("inf") # init val
    
    for slice_size in range(1, len(arr)+1):
        for step in range(0, len(arr)+1-slice_size):
            n = xor_sum(arr[i] for i in range(step,step+slice_size))
            
            freq[n] = freq.get(n,0)+1
            if freq[n] >= max_c and (n < min_n or freq[n]> max_c):
                min_n = n
                max_c = freq[n]

    return  min_n, freq[min_n]

But it times out since it's ~O(n^3). I feel like there is some math trick, can someone explain the solution to me? I tried to read some solutions in the discussion but I didn't quite get them.

Problem copy:

Consider an array, A, of n integers (A=a₀,a₁,...,a_n-1). We take all consecutive subsequences of integers from the array that satisfy the following:
{a_i,a_i+1,...,a_j-1,a_j}, where 0≤i≤j≤n

For each subsequence, we apply the bitwise XOR (⊕) operation on all the integers and record the resultant value.

Given array A, find the XOR sum of every subsequence of A and determine the frequency at which each number occurs. Then print the number and its respective frequency as two space-separated values on a single line.

Output Format

Print 2 space-separated integers on a single line. The first integer should be the number having the highest frequency, and the second integer should be the number's frequency (i.e., the number of times it appeared). If there are multiple numbers having maximal frequency, choose the smallest one.

Constraints

• 1≤n≤10⁵
• 1≤a_i<2¹⁶

Answer 1

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 definition 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).

Answer 2

I propose this as a functional bitwise xor sum of subsequence between indexes l and r.

# Sequence (An)n is defined by A(n)=A(n-1)^n with ^ bitwise xor operation
# Compute bitwise xor of subsequence of consecutive sequence elements
# ...between a left and right index, A(k), k=l..r, i.e., A(l)^A(l+1)^...^A(r-1)^A(r)

#!/bin/python3

import math
import os
import random
import re
import sys

# Hints ##
# Notice how to bitwise xor of consecutive integers in range 1..n simplifies
# e.g., printing out bitwise xor of consecutive integers in 1..13
# 0 1 3 0 4 1 7 0 8 1 11 0 12
# ...someting happens, some pattern wrt n%4
# Method comp_xor_in_range(n) shows a corresponding fast implementation
# for bitwise xor of consecutive integers in range 1..n

# Then, notice how bitwise xor of (An)n subsequence elements simplifies, depending on 
# ...parity of left index and parity of number of elements in subsequence 

# compute bitwise xor of elements in range 1..n
def comp_xor_in_range(l):
    res=l%4
    if res==0:
        return l
    elif res==1:
        return 1
    elif res==2:
        return l+1
    else:
        return 0

# compute bitwise xor of elements in range l(eft)..r
def comp_xor_in_range_lr(l,r):
    return comp_xor_in_range(l-1)^comp_xor_in_range(r)

#compute bitwise xor of even elements in range 1..n 
def comp_xor_even_in_range(l):
    res=l%2
    if res==0:
        return 2*comp_xor_in_range(l//2)
    else:
        return 2*comp_xor_in_range((l-1)//2)
    
# compute bitwise xor of even elements in range l(eft)..r
def comp_xor_even_in_range_lr(l,r):
    return comp_xor_even_in_range(l-1)^comp_xor_even_in_range(r)

#compute bitwise xor of odd elements in range 1..n 
def comp_xor_odd_in_range(l):
    return comp_xor_even_in_range(l)^comp_xor_in_range(l)

# compute bitwise xor of even elements in range l(eft)..r
def comp_xor_odd_in_range_lr(l,r):
    return comp_xor_odd_in_range(l-1)^comp_xor_odd_in_range(r)

def xorSequence(l, r):
    n=r-l+1
    
    if n%2==0:
        if l%2==0:
            return comp_xor_odd_in_range_lr(l,r)
        else:
            return comp_xor_even_in_range_lr(l,r)
    else:
        if l%2==0:
            return comp_xor_even_in_range_lr(l+1,r)^comp_xor_in_range(l)
        else:
            return comp_xor_odd_in_range_lr(l+1,r)^comp_xor_in_range(l)


    
if __name__ == '__main__':
    fptr = open(os.environ['OUTPUT_PATH'], 'w')

    q = int(input())

    for q_itr in range(q):
        lr = input().split()

        l = int(lr[0])

        r = int(lr[1])

        result = xorSequence(l, r)

        fptr.write(str(result) + '\n')

    fptr.close()

Answer 3

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.)

Answer 4

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)

Answer 5

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.