Преобразование Уолша-Адамара31.03.2024 19:30

На сайте hackerrank.com есть отличная задача.
По заданному массиву short[] A; найти максимальное количество его подмассивов, xor элементов которых будет одинаковым. Сам этот xor тоже нужно найти.

Максимальная длина массива равна 105, так что квадратичный алгоритм не укладывается в лимит по времени исполнения. Я в своё время с этой задачей не справился и сдался, решив подсмотреть авторское решение. И в этот момент я понял почему не справился — автор предлагал решать задачу через дискретное преобразование Фурье.

Подробная формулировка задачи

Сперва стоит подробно описать задачу, чтобы быть уверенным, что все поняли её одинаково. Задан массив A целых чисел, каждое из которых находится в диапазоне от 1 до 65535 (границы включены). Допустим, что это массив $\{1,2,3\}$ . Рассмотрим все его подмассивы:

2 подмассива имеют xor равный 1, так же 2 других подмассива имеют xor равный 3.
2 — максимальное количество подмассивов с одинаковым xor-ом. Сам же xor выбираем равным 1, как минимальное значение из всех альтернатив (1 или 3). Итого решение задачи — пара (1, 2) .

Нужно лишь научиться делать то же самое для массивов размером вплоть до 100000. Длину массива будем обозначать как N.

Наивное решение

Наивное решение предполагает двойной цикл по вспомогательному массиву B, строящийся следующим образом:

B[0] = 0
B[i] = B[i-1] ^ A[i-1]

Массив B имеет длину N+1. При таком построении каждый B[i] (при i > 0) будет равен xor-у всех элементов A[k] для 0 <= k < i. Данное свойство очень удобно, потому что в этом случае xor подмассива [A[i], ..., A[j]] будет равен B[j+1] ^ B[i], т.е. вычисляем за константное время. Сам же массив строится за линейное время, тоже быстро.

Остаётся лишь посчитать xor всех подмассивов. Псевдокод будет выглядеть следующим образом (65536 — это 216):

int[] frequencies = new int[65536];

for (int i = 0; i < N; i++) {
    for (int j = i; j < N; j++) {
        ++frequencies[B[j + 1] ^ B[i]];
    }
}

int maxFrequency = max(frequencies);
int xor = indexOf(frequencies, maxFrequency);

Тройного цикла мы избежали, а вот двойной остался, и это плохо. В нём будет $\frac{N(N+1)}{2}$ итераций, что в худшем случае равно примерно 5×109.

Учитывая, с какой частотой обычно работают процессоры, можно заключить, что этот цикл будет исполняться несколько секунд. На практике решение не укладывается во временные рамки, настроенные на сайте (4 секунды для Java, например).

Какое решение предлагает автор

Официальное решение из секции Editorial

Неправильное.

То есть как неправильное, последующий код в Editorial приведён верный, а вот текстовое описание со скриншота — полнейшая бессмыслица, с кодом не совпадающая.
Хотя и к коду у меня тоже есть вопросы, но эту часть опустим, скажу просто, что он там странный.

Я буду своими словами объяснять то, что на самом деле имелось в виду, а не то, что там написано. И сделаю это не за 3 строки текста. Идея финального алгоритма такая:

long[] C = new long[65536];

for (int i = 0; i <= N; i++) {
    ++C[B[i]];
}

convolute(C);
C[0] -= N + 1;

long maxFrequency = max(C) / 2;
long xor = indexOf(C, maxFrequency * 2);

Решение состоит из 2-х частей: свёртки вспомогательного массива C
(а не массива B, как было заявлено) и выражению C[0] -= N + 1. К моменту вызова max(C) массив C должен содержать в точности те же данные, что и массив frequencies в наивном решении, но умноженные на 2.

Приблизимся к пониманию этого кода, подумав, как по изначальному массиву C вообще можно получить ответ. Это и будет являться ключом к решению.
В наивном решении мы делали ++frequencies[k] для всех k == B[i] ^ B[j],
где 0 <= i < j <= N (здесь я произвёл замену переменной j + 1 -> j, чтобы не таскать за собой +1, так меньше путаницы).
С помощью ++ в итоге мы считали количество пар (i, j), для которых B[i] ^ B[j] одинаковы и равны k. Иными словами, для всех i < j и k == B[i] ^ B[j] мы должны просуммировать единицы, чтобы получить frequencies[k].

$F_k = \sum_{i<j,\:b_i \oplus b_j = k}1$

Рассмотрим тот же пример, что был показан при подробной формулировке задачи.
A = [1, 2, 3] B = [0, 1, 3, 0] ([0, 0^1, 0^1^2, 0^1^2^3]) C = [2, 1, 0, 2, ...]Далее в индексах при 1 написаны соответствующие им пары (i, j):

$\begin{align} & F_0 = 1_{0,3} = 1 & (b_0\oplus b_3 = 0)\\ & F_1 = 1_{0,1}+1_{1,3} = 2 & (b_0\oplus b_1 = b_1\oplus b_3 = 1)\\ & F_2 = 1_{1,2} = 1 & (b_1\oplus b_2 = 2)\\ & F_3 = 1_{0,2} + 1_{2,3} = 2 & (b_0\oplus b_2 = b_2\oplus b_3 = 3) \end{align}$

Вроде бы сходится со значениями, которые мы видели ранее.

Начнём преобразовывать эту формулу до тех пор, пока не получится что-нибудь путное. В первую очередь стоит избавиться от i < j: такая сложная индексация в сумме сильно затрудняет формулу. Рассмотрим альтернативное значение:

$F_k^{'}=\sum_{b_i\oplus b_j=k}1$

В нём порядок i и j игнорируется, а также допускается их полное равенство. Это можно показать явно:

$F_k^{'}=\sum_{i<j,\:b_i\oplus b_j=k}1 + \sum_{i>j,\: b_i\oplus b_j=k}1 + \sum_{i=j,\: b_i\oplus b_j = k}1» src=«https://habrastorage.org/getpro/habr/upload_files/6e9/f21/133/6e9f21133f6a237b4eaa468d75d23379.svg» /></p> <p>Последнее слагаемое имеет смысл только для <code>k == 0</code>, поскольку для всех других <code>k</code> там будет просто пустая сумма (<code>xor</code> равен <code>0</code> только при равенстве аргументов). А равняться это самое слагаемое будет количеству различных <code>i</code>, по которым мы пробегаемся, т.е. длине массива <code>B</code>, она равна <code>N+1</code>. Далее, суммы для <code>i<j</code> и для <code>i>j</code> численно равны, так как <code>xor</code> коммутативен. Учитывая всё вышесказанное, уравнения можно переписать в следующем виде: </p> <p><img alt=$

А значит, нам достаточно вычислять более простой $F^{'}$ .
Что, если некоторые B[i] совпадают? Кажется, что в общем случае это справедливое допущение. Если бы мы заранее знали, что B[i_1] == B[i_2] == ... == B[i_n], то подсчёт количества всех B[i_k] ^ B[j] был бы тривиален, таких пар было бы ровно n штук.

Процедуру вычисления можно видоизменить, если сделать её следующим образом:
для всех i таких, что все B[i] попарно различны и k == B[i] ^ B[j], мы должны просуммировать count(B[i]), чтобы получить $F_k^{'}$ . Похоже на то, что было раньше, но суммируются не единицы, да и i учитываются не все.

$F_k^{'} = \sum_{b_i \oplus b_j = k,\:b_i\:is\:unique}count(b_i)$

Фраза bi is unique не совсем корректна, но, надеюсь, она даёт интуитивное понимание происходящего. Что, если провести такую же группировку и по индексам j? Получим следующее:
для всех i таких, что все B[i] попарно различны, всех j таких, что все B[j] попарно различны и k == B[i] ^ B[j], мы должны просуммировать count(B[i]) * count(B[j]), чтобы получить $F_k^{'}$ .

$F_k^{'} = \sum_{b_i \oplus b_j = k,\:b_i\:is\:unique,\:b_j\:is\:unique}count(b_i)\cdot count(b_j)$

Если посмотреть внимательно на код решения, то будет видно, что массив C в нём как раз вычисляет count, т.е. C[B[i]] == count(B[i]). Далее, i и j можно будет вообще выбросить из обсуждения, сделав небольшую хитрость.

Если какое-то значение x не встречается в массиве B, то C[x] == 0. Если же какое-то значение x встречается в B, то оно встречается в точности C[x] раз. Это хорошее свойство. Учитывая его, предлагаю посмотреть, что будет при вычислении подобной суммы:

$\sum_{p\oplus q = k}count(p)\cdot count(q)$

которую можно переписать немного иначе, избавившись от буквы q:

$\sum_{p=0}^{65535}count(p)\cdot count(p\oplus k)$

Если p или q отсутствует в массиве B, то соответствующее значение count будет равно нулю и слагаемое будет игнорироваться. В противном случае p и q будут соответствовать уникальным значениям в B. Другими словами, это буквально та же самая сумма.

$F_k^{'} = \sum_{p\oplus q=k}count(p)\cdot count(q)$

Мы можем написать следующий код, который должен быть эквивалентен строке convolute(C):

long[] tmp = new long[65536];

for (int p = 0; p < 65536; p++) {
    for (int q = 0; q < 65536; q++) {
        tmp[p ^ q] += C[p] * C[q];
    }
}

C = tmp;

Именно такая сумма произведений называется свёрткой массива C с самим собой.
Так же стало понятно, откуда взялись C[0] -= N + 1 и последующее деление на 2, из соотношения между $\mathcal{F}$ и $\mathcal{F^{'}}$ .

Вычисляться свёртка будет с помощью быстрого преобразования Уолша-Адамара (FWHT), работающего за O(M*log(M)) (у нас M = 65536) и ещё одного волшебного действия, о котором чуть позже.
Большая часть статьи как раз будет посвящена объяснениям того, что это за преобразование, что такое свёртка в более общем случае, как они связаны да и что такое дискретные преобразования Фурье в целом.

Дискретное преобразование Фурье

Любое дискретное преобразование Фурье (DFT) — вещь достаточно простая. На входе есть вектор длины, на выходе есть вектордлины, посередине находится умножение матрицы на вектор. Или вектора на матрицу, это не принципиально. Начну я с классического дискретного преобразования Фурье, а к преобразованию Уолша-Адамара (WHT) перейду попозже.

Для самых обычных векторов $x=\{x_0, \ldots, x_{n-1}\}$ размерностииспользуется матрица

$\mathcal{F} = \begin{pmatrix} 1 & 1 & 1 & \dots & 1 \\ 1 & \omega & \omega^2 & \dots & \omega^{n-1} \\ 1 & \omega^2 & \omega^{2\cdot 2} & \dots & \omega^{2(n-1)} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \omega^{n-1} & \omega^{(n-1)\cdot 2} & \dots & \omega^{(n-1)^2} \\ \end{pmatrix}$

где $\omega=e^{-\frac{2\pi i}{n}}$ — корень-й степени из 1. Это в буквальном смысле полное определение дискретного преобразования Фурье. Во всяком случае оно эквивалентно всем остальным.
Есть ещё обратное преобразование, с помощью которого из результирующего вектораможно обратно получить исходный вектор, оно выражается обратной матрицей

$\mathcal{F}^{-1} = \frac{1}{n} \begin{pmatrix} 1 & 1 & 1 & \dots & 1 \\ 1 & \omega^{-1} & \omega^{-2} & \dots & \omega^{-(n-1)} \\ 1 & \omega^{-2} & \omega^{-2\cdot 2} & \dots & \omega^{-2(n-1)} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \omega^{-(n-1)} & \omega^{-(n-1)\cdot 2} & \dots & \omega^{-(n-1)^2} \\ \end{pmatrix}$

в которой $\frac{1}{n}$ — это нормирующий коэффициент.

$\begin{align} & X = \mathcal{F}x\\ & x = \mathcal{F}^{-1}X \end{align}$

Тут дело в том, что определитель матрицы $\mathcal{F}$ равен $\sqrt{n}$ .

В общем виде элементы матрицы выглядят следующим образом:

$\mathcal{F}_{p,q} = \omega^{pq} = e^{-\frac{2\pi i p q}{n}}, 0\le p<n, 0\le q < n$

Преобразование системы координат

Про дискретные преобразования Фурье часто пишут, что они раскладывают сигнал на гармоники, или типа того. Все эти интерпретации в рамках данной статьи меня не интересуют. Чисто алгебраически дискретное преобразование Фурье — это просто перевод вектора из одного базиса в другой (частный случай аффинного преобразования). Просто базис подбирается таким образом, что перевод в него начинает обладать некоторыми волшебными свойствами.

Сейчас будет небольшое отступление, которое объясняет как этот перевод работает на деле.

Для DFT базисом является система векторов $f_q=\{\frac{1}{n}, \frac{1}{n}e^{\frac{2\pi i q}{n}}, \frac{1}{n}e^{\frac{2\pi i q \cdot 2}{n}},\dots, \frac{1}{n}e^{\frac{2\pi i q(n-1)}{n}}\}$ :

$x=\sum_q{X_q f_q}$

X_q — это просто коэффициенты в разложении векторапо векторам f_q . Более того, то, что тут написано — это буквально формула обратного DFT в более развёрнутом виде. Теперь увидим, почему X_q , посчитанные через прямое DFT, удовлетворяют данному равенству. Для этого, в первую очередь, стоит заметить, что

$f_q=\{\frac{1}{n}, \frac{1}{n}\overline{\omega}^q, \frac{1}{n}\overline{\omega}^{2q},\dots, \frac{1}{n}\overline{\omega}^{(n-1)q}\}$

Верхняя черта над $\omega$ — это символ комплексного сопряжения, превращающий $\omega=e^{-\frac{2\pi i}{n}}$ в $\overline{\omega}=e^{\frac{2\pi i}{n}}$ . Это сопряжение можно переписать иначе:

$\overline{f}_p=\{\frac{1}{n}, \frac{1}{n}\omega^p, \frac{1}{n}\omega^{2p},\dots, \frac{1}{n}\omega^{(n-1)p}\}$

Данное выражение уже подозрительно похоже на строку матрицы $\mathcal{F}$ , разве что наумножить нужно. Да и вообще, матрица $\mathcal{F}$ — унитарная, с поправкой на нормирующий коэффициент.

$\mathcal{F}^{-1} = \frac{1}{n}\mathcal{F}^* = \frac{1}{n}\overline{\mathcal{F}^\text{T}}$

Если матрицу нормировать, то она станет воистину унитарной. Распишем прямое преобразование:

$\begin{align} & X = \sum_p x_p n \overline{f}_p\\ & X_q = \sum_p x_p n \overline{f}_{p,q} = \sum_p x_p \omega^{pq} = n(x\cdot f_q)\end{align}$

Здесь $(a\cdot b)$ — это скалярное произведение двух комплексных векторов, равное $\sum_i{a_i\overline{b}_i}$ .

Один из множителей в каждом слагаемом обязан быть заменённым на комплексно-сопряжённый, определение такое. Если собрать всё в кучу, то выйдет следующее:

$x=\sum_q{n(x\cdot f_q) f_q}$

А это уже просто разложение вектора по ортогональному базису, которое в общем виде выглядит так:

$x = \sum_q{\frac{(x\cdot f_q) f_q}{(f_q\cdot f_q)}} = \sum_q{\frac{(x\cdot f_q) f_q}{{\lvert f_q\rvert}^2}}$

Проверить ортогональность несложно, нужно лишь показать, что $(f_p \cdot f_q) = 0 \Leftrightarrow p \neq q$ , для этого достаточно знать формулу суммы геометрической прогрессии.

Теорема о свёртке

Свёртка — это операция, позволяющая из 2-х векторов получить 3-й.
Записывается она символом $\ast$ :

$\begin{align} &c = a \ast b\\ &c_k = \sum_{q}{a_q b_{(k-q)\pmod{n}}} \end{align}$

Как и в случае с самим DFT, у меня не стоит цели рассказать, зачем свёртка нужна. Меня просто интересуют её конкретные свойства. Одно из них — это то, как свёртка связана с DFT:

$\mathcal{F}(a \ast b)_p = (\mathcal{F}a)_p \cdot (\mathcal{F}b)_p$

Дискретное преобразование Фурье превращает свёртку в покомпонентное умножение векторов, так же известное как произведение Адамара:

$\mathcal{F}(a \ast b) = (\mathcal{F}a)\circ(\mathcal{F}b)$

На практике это означает, что свёртку можно вычислять с помощью прямого DFT, умножения результатов и применению обратного DFT. И это быстрее прямого вычисления свёртки по формуле из определения, поскольку для DFT есть эффективный алгоритм, называемый быстрым преобразованием Фурье (FFT).

Предлагаю разобрать, почему это свойство работает, потом мы наконец-то перейдём к обещанному преобразованию Уолша-Адамара. В первую очередь перепишем формулу свёртки в немного другой форме:

$c_k=\sum_{p+q\equiv k\pmod{n}}a_p b_q$

Мне такой вариант записи всегда нравился больше, чем классический. Более симметрично, что ли. Подействуем на этот вектор дискретным преобразованием Фурье:

$\mathcal{F}(a \ast b) = \sum_{k}c_k\mathcal{F}_k = \sum_{k}(\sum_{p+q\equiv k\pmod{n}}a_p b_q)\mathcal{F}_k$

Здесь $\mathcal{F}_k$ — это k-я строка матрицы $\mathcal{F}$ (или столбец — это неважно, матрица симметричная).
Каждый компонент результирующего вектора выглядит следующим образом:

$\mathcal{F}(a \ast b)_s = \sum_{k}(\sum_{p+q\equiv k\pmod{n}}a_p b_q)\mathcal{F}_{k,s} = \sum_{p}\sum_{q}{a_p b_q \mathcal{F}_{(p+q)\pmod{n}, s}}$

В выражении выше я произвёл замену индексов суммирования. Надеюсь, больших вопросов она не вызывает.
Как будет выглядеть произведение Адамара, стоящее справа в теореме о свёртке?

$(\mathcal{F}a)_s\cdot(\mathcal{F}b)_s = (\sum_{p}a_p\mathcal{F}_{p,s})(\sum_{q}b_q\mathcal{F}_{q,s}) = \sum_{p}\sum_{q}{a_p b_q \mathcal{F}_{p, s}\mathcal{F}_{q, s}}$

Очень похоже. Вот если бы $\mathcal{F}_{(p+q)\pmod{n}, s}$ равнялось $\mathcal{F}_{p, s}\cdot \mathcal{F}_{q, s}$ …

Хотя погодите, оно и равняется, ведь $\mathcal{F}_{p,s}=\omega^{ps}$ , $\omega^n=1$ и выражение выше следует напрямую из свойств показательной функции.

$\mathcal{F}_{(p+q)\pmod{n}, s}=\omega^{((p+q)\pmod{n})s}=\omega^{(p+q)s}=\omega^{ps}\cdot\omega^{qs}=\mathcal{F}_{p, s}\cdot \mathcal{F}_{q, s}$

Вот оно, то волшебное свойство матрицы $\mathcal{F}$ , которое позволяет теореме о свёртке работать! А это значит, что если взять другую матрицу с тем же свойством, то теорема о свёртке выполнится и для неё.

Отсюда есть несколько интересных выводов:

$\mathcal{F}_{p, s} = \mathcal{F}_{0, s} \cdot \mathcal{F}_{p, s} \Rightarrow \mathcal{F}_{0, s} = 1\\ \mathcal{F}_{p, s} = \mathcal{F}_{1, s} \cdot \mathcal{F}_{p-1, s} \Rightarrow \mathcal{F}_{p, s} = \mathcal{F}_{1, s}^{p}$

Симметричные формулы для переставленных индексов тоже верны. Матрица для DFT полностью определяется элементом $\mathcal{F}_{1, 1}$ : $\mathcal{F}_{p, q} = \mathcal{F}_{1, 1}^{pq}$ , иначе теорема о свёртке просто не выполнится.

Преобразование Уолша-Адамара

WHT является частным случаем преобразования Фурье на группах, которое, в свою очередь, является обобщением дискретного преобразования Фурье.

Преобразование Уолша-Адамара тоже является переводом вектора размерностииз одного базиса в другой и тоже выражается с помощью умножения матрицы на вектор, но при дополнительных ограничениях на:

$n = 2^m, m \in \mathbb{N}$

Вместов формулах принято использовать сам показатель степени.

$X = H_mx\\ x = H_m^{-1}X$

Матрица преобразования, называемая матрицей Адамара, вводится рекуррентно, начиная со случая m=1 (0 я проигнорирую, слишком скучно). При этом, вначале мы получим в точности матрицу обычного DFT. Матрицы бОльших размерностей строятся из матриц меньших размерностей.

$\begin{align}&H_1=\mathcal{F}=\begin{pmatrix}1&1\\1&-1\end{pmatrix}\\ &H_m=\begin{pmatrix}H_{m-1}&H_{m-1}\\H_{m-1}&-H_{m-1}\end{pmatrix}\end{align}$

Например

$H_2 = \begin{pmatrix} \phantom{-}1 & \phantom{-}1 & \phantom{-}1 & \phantom{-}1 \\ \phantom{-}1 & -1 & \phantom{-}1 & -1 \\\phantom{-}1 & \phantom{-}1 & -1 & -1 \\\phantom{-}1 & -1 & -1 & \phantom{-}1\end{pmatrix}\\ H_3 = \begin{pmatrix} \phantom{-}1 & \phantom{-}1 & \phantom{-}1 & \phantom{-}1 & \phantom{-}1 & \phantom{-}1 & \phantom{-}1 & \phantom{-}1 \\ \phantom{-}1 & -1 & \phantom{-}1 & -1 & \phantom{-}1 & -1 & \phantom{-}1 & -1\\ \phantom{-}1 & \phantom{-}1 & -1 & -1 & \phantom{-}1 & \phantom{-}1 & -1 & -1\\ \phantom{-}1 & -1 & -1 & \phantom{-}1 & \phantom{-}1 & -1 & -1 & \phantom{-}1 \\ \phantom{-}1 & \phantom{-}1 & \phantom{-}1 & \phantom{-}1 & -1 & -1 & -1 & -1 \\ \phantom{-}1 & -1 & \phantom{-}1 & -1 & -1 & \phantom{-}1 & -1 & \phantom{-}1\\ \phantom{-}1 & \phantom{-}1 & -1 & -1 & -1 & -1 & \phantom{-}1 & \phantom{-}1\\ \phantom{-}1 & -1 & -1 & \phantom{-}1 & -1 & \phantom{-}1 & \phantom{-}1 & -1 \\ \end{pmatrix}$

Начнём со сходств с дискретным преобразованием Фурье. Матрица— ортогональная, т.е. тоже описывает перевод в ортогональный базис, это легко проверить самостоятельно. Выражение для обратного преобразования идентично выражению для DFT:

$H^{-1} = \frac{1}{2^m}H^*$

Знак * можно вообще убрать из формулы выше, поскольку матрица преобразования симметричная и содержит только вещественные числа.

Теперь различия. Для WHT не выполняется теорема о свёртке, поскольку $H_{(p + q)\pmod{n}, s} \neq H_{p, s} \cdot H_{q, s}$ . Поскольку WHT — это вариант обобщённого DFT, то и его теорема о свёртке должна быть вариантом обобщения теоремы о свёртке для DFT. Попробуем прийти к этому обобщению.

Если смотреть на матрицу «рекурсивно», то видно, что элементы меняют знак, если они находятся в правой нижней подматрице во время «рекурсивного спуска». Если вычислить, сколько раз меняется знак, то можно будет легко получить значение элемента массива как $-1^{count}$ .

Для матрицы 2 на 2 это просто сделать, $H_{p, q} = -1^{pq}, p \in \{0, 1\}, q \in \{0, 1\}$ . Что же касается матриц больших размеров

$H_m = \begin{pmatrix}H_{m-1}&H_{m-1}\\H_{m-1}&-H_{m-1}\end{pmatrix}$

Все индексы в нижней половине матрицы больше либо равны 2m-1, аналогично все индексы в правой половине матрицы больше либо равны 2m-1. Такие индексы имеют единичный бит в позиции m-1 , а значит при рекурсивном спуске нам нужно будет домножать значение элемента на $-1^{bit_{m-1}(p) \cdot bit_{m-1}(q)}$ (по аналогии с матрицей 2 на 2).

Если перемножить эти выражения для всех битов чисели, то получим следующее:

$H_{p, q} = -1^{bit\_count(p \wedge q)}$

Тут главное не запутаться. В программировании & — это конъюнкция, а ^ — это исключающее или. В математике $\wedge$ — это конъюнкция, а $\oplus$ — это исключающее или. Под записью $p \wedge q$ я имею ввиду побитовую конъюнкцию. А выражение $bit\_count(p \wedge q)$ фактически равно скалярному произведениюи, если представлять их себе как битовые вектора.

Мы снова получили показательную функцию в выражении для элементов матрицы, а значит вполне можем рассчитывать на свойства, похожие на свойства матрицы $\mathcal{F}$ .

Свёртка для WHT

Конъюнкция битовых векторов — это их побитовое произведение. А исключающее или битовых векторов — это их побитовое сложение по модулю 2. Эти 2 операции ведут себя очень схоже с обычными умножением и сложением. В частности, выполнено свойство дистрибутивности:

$(p \oplus q)\wedge s = (p \wedge s) \oplus (q \wedge s)$

Так же важно обратить внимание, что

$bit\_count(p\oplus q) \equiv (bit\_count(p)+bit\_count(q))\pmod{2}$

Учитывая это, мы можем снова воспользоваться свойствами показательных функций и заключить следующее:

$H_{(p\oplus q, s)} = -1^{bit\_count((p \oplus q)\wedge s)} = -1^{bit\_count(p\wedge s)} \cdot -1^{bit\_count(q\wedge s)} = H_{p, q}\cdot H_{q, s}$

Как помним, для DFT у нас была похожая формула, только сложение там было другое:

$\mathcal{F}_{(p+q)\pmod{n}, s} = \mathcal{F}_{p, s} \cdot \mathcal{F}_{q, s}$

Именно это свойство использовалось для доказательства теоремы о свёртке, где формулы для свёртки $c = a \ast b$ были следующими:

$\begin{align} & (\mathcal{F}c)_k = (\mathcal{F}a)_k \cdot (\mathcal{F}b)_k\\ &c_k = \sum_{p+q\equiv k\pmod{n}}{a_p b_q}\end{align}$

По аналогии мы можем получить теорему о свёртке для преобразование Уолша-Адамара:

$\begin{align} & (Hc)_k = (Ha)_k \cdot (Hb)_k \\ & c_k = \sum_{p \oplus q = k}{a_p b_q}\end{align}$

Теорема о свёртке выполняется, просто сумма индексов в ней не $+\pmod{n}$ , а $\oplus$ . И это в точности та же самая формула, которую должен вычислить метод convolute.

Используя другие операции вместо $\oplus$ или сложения по модулю можно получать другие варианты дискретных преобразований Фурье. Их, понятное дело, может быть очень много.

Быстрое преобразование Уолша-Адамара

$H(a \ast b) = (Ha)\circ(Hb)$

Чтобы вычислить свёртку двух векторов, пользуясь теоремой о свёртке, нужно произвести WHT над аргументами, вычислить произведение Адамара и выполнить обратное преобразование. В нашем частном случае это приведёт к подобному коду:

static final int M = 65536;

static void convolute(long[] C) {
    wht(C);

    for (int i = 0; i < M; i++) {
        C[i] *= C[i];
    }

    wht(C);

    for (int i = 0; i < M; i++) {
        C[i] /= M;
    }
}

Помним, что $H_m^{-1} = \frac{1}{2^m}H_m$ . Отсюда и деление в конце.

Всё, что нам осталось — это реализовать метод wht:

static void wht(long[] C) {
    for (int h = 1; h < M; h *= 2) {
        for (int i = 0; i < M; i += (h * 2)) {
            for (int j = 0; j < h; j++) {
                long x = C[i + j],
                     y = C[i + j + h];

                C[i + j]     = x + y;
                C[i + j + h] = x - y;
            }
        }
    }
}

Данный алгоритм вовсю эксплуатирует рекурсивную структуру матриц Адамара, буквально повторяя рекуррентное определение матриц при вычислении. Внешний цикл — логарифмический, а 2 внутренних цикла вместе делают один проход по массиву, просто не по порядку. Подробный разбор этой реализации выходит за рамки статьи.

Заключение

Хорошая задача. Сложная. Знаете, что в ней самое интересное? То, что я вас обманул. Не специально, конечно же.

Ближе к началу статьи я говорил, что двойной цикл по массиву B выполнит примерно 5×109 итераций и это слишком много. Это правда. А сколько итераций выполнил бы двойной цикл по массиву C, если бы мы вычисляли свёртку напрямую?
Примерно 2×109, если итерироваться для j >= i, это почти в 2 с половиной раза меньше. Для Java ограничение по времени — 4 секунды. Может прокатит?

public static void main(String[] args) {
    int N = readUnsignedInt();

    int B[] = new int[N + 1];
    for (int i = 0; i < N; i++) {
        B[i + 1] = B[i] ^ readUnsignedInt();
    }

    int C[] = new int[M];
    for (int i = 0; i <= N; i++) {
        ++C[B[i]];
    }

    int tmp[] = C;
    C = new int[M];
    for (int i = 0; i < M; i++) {
        for (int j = i; j < M; j++) {
            C[i ^ j] += tmp[i] * tmp[j];
        }
    }

    C[0] = (C[0] - N - 1) / 2;

    int max = Integer.MIN_VALUE, imax = 0;
    for (int i = 0; i < M; i++) {
        if (C[i] > max) {
            max = C[i];
            imax = i;
        }
    }
    System.out.print(imax);
    System.out.print(' ');
    System.out.println(max);
}

Прокатило!

Получается, что задача решается и без преобразования Уолша-Адамара. Т.е. тот факт, что я изначально не смог решить эту задачу и подглядел ответ, никак не связан с тем, что я не знаком c WHT. Он связан с тем, что я поленился подольше посидеть над задачей, и признавать это сейчас немного стыдно.

А ведь автор задачи мог исключить возможность такого решения, если бы вместо 16-битных чисел использовал, к примеру, 17-битные. Вместе с тем фактом, что в Editorial написана какая-то ерунда, я заключаю, что автор задачи — на самом деле не автор. Почти наверняка он спёр её с другого сайта, может даже не разобравшись до конца в нужной математике и не подогнав условия под тайминги на новом сайте. Искать эту задачу в других источниках мне уже не хочется, так что эти мысли останутся просто гипотезой.

В любом случае я рад, что разобрался в предложенном решении, да и дискретное преобразование Фурье на деле оказалось намного более простым, чем я считал ранее.