diff --git a/maths/fast_walsh_hadamard_transform.py b/maths/fast_walsh_hadamard_transform.py new file mode 100644 index 000000000000..1c7d447e3c0a --- /dev/null +++ b/maths/fast_walsh_hadamard_transform.py @@ -0,0 +1,209 @@ +""" +Fast Walsh-Hadamard Transform (FWHT) for Bitwise Convolutions. + +Reference: https://en.wikipedia.org/wiki/Fast_Walsh%E2%80%93Hadamard_transform +Reference: https://cp-algorithms.com/algebra/walsh-hadamard-transform.html + +Computes bitwise XOR, AND, and OR convolutions of two numeric sequences in O(N log N) +time, where N is a positive power of 2. +""" + + +def fwht_xor(sequence: list[int], inverse: bool = False) -> list[int]: + """ + Perform Fast Walsh-Hadamard Transform (or inverse) for XOR operation. + + Time Complexity: O(N log N) + + >>> fwht_xor([1, 2, 3, 4]) + [10, -2, -4, 0] + >>> fwht_xor([10, -2, -4, 0], inverse=True) + [1, 2, 3, 4] + >>> fwht_xor([1, 2, 3]) + Traceback (most recent call last): + ... + ValueError: Length of sequence must be a positive power of 2. + >>> fwht_xor([]) + Traceback (most recent call last): + ... + ValueError: Length of sequence must be a positive power of 2. + >>> fwht_xor([1, 2, 3, 5], inverse=True) + Traceback (most recent call last): + ... + ValueError: Inverse XOR transform requires elements divisible by sequence length. + """ + sequence_length = len(sequence) + if sequence_length == 0 or (sequence_length & (sequence_length - 1)) != 0: + raise ValueError("Length of sequence must be a positive power of 2.") + + result = list(sequence) + half_block = 1 + while half_block < sequence_length: + block_size = half_block * 2 + for block_start in range(0, sequence_length, block_size): + for offset in range(half_block): + left_index = block_start + offset + right_index = left_index + half_block + left_val = result[left_index] + right_val = result[right_index] + result[left_index] = left_val + right_val + result[right_index] = left_val - right_val + half_block *= 2 + + if inverse: + if any(element % sequence_length != 0 for element in result): + raise ValueError( + "Inverse XOR transform requires elements divisible by sequence length." + ) + result = [element // sequence_length for element in result] + return result + + +def xor_convolution(sequence_a: list[int], sequence_b: list[int]) -> list[int]: + """ + Compute bitwise XOR convolution C[k] = sum_{i ^ j = k} (A[i] * B[j]). + + Time Complexity: O(N log N) + + >>> xor_convolution([1, 2], [3, 4]) + [11, 10] + >>> xor_convolution([1, 2], [3]) + Traceback (most recent call last): + ... + ValueError: Input sequences must have equal length. + """ + if len(sequence_a) != len(sequence_b): + raise ValueError("Input sequences must have equal length.") + + transformed_a = fwht_xor(sequence_a) + transformed_b = fwht_xor(sequence_b) + pointwise_product = [ + transformed_a[index] * transformed_b[index] for index in range(len(sequence_a)) + ] + return fwht_xor(pointwise_product, inverse=True) + + +def fwht_or(sequence: list[int], inverse: bool = False) -> list[int]: + """ + Perform Fast Walsh-Hadamard Transform for OR operation. + + Time Complexity: O(N log N) + + >>> fwht_or([1, 2]) + [1, 3] + >>> fwht_or([1, 3], inverse=True) + [1, 2] + >>> fwht_or([1, 2, 3]) + Traceback (most recent call last): + ... + ValueError: Length of sequence must be a positive power of 2. + """ + sequence_length = len(sequence) + if sequence_length == 0 or (sequence_length & (sequence_length - 1)) != 0: + raise ValueError("Length of sequence must be a positive power of 2.") + + result = list(sequence) + half_block = 1 + while half_block < sequence_length: + block_size = half_block * 2 + for block_start in range(0, sequence_length, block_size): + for offset in range(half_block): + left_index = block_start + offset + right_index = left_index + half_block + if not inverse: + result[right_index] += result[left_index] + else: + result[right_index] -= result[left_index] + half_block *= 2 + + return result + + +def or_convolution(sequence_a: list[int], sequence_b: list[int]) -> list[int]: + """ + Compute bitwise OR convolution C[k] = sum_{i | j = k} (A[i] * B[j]). + + Time Complexity: O(N log N) + + >>> or_convolution([1, 2], [3, 4]) + [3, 18] + >>> or_convolution([1, 2], [3]) + Traceback (most recent call last): + ... + ValueError: Input sequences must have equal length. + """ + if len(sequence_a) != len(sequence_b): + raise ValueError("Input sequences must have equal length.") + + transformed_a = fwht_or(sequence_a) + transformed_b = fwht_or(sequence_b) + pointwise_product = [ + transformed_a[index] * transformed_b[index] for index in range(len(sequence_a)) + ] + return fwht_or(pointwise_product, inverse=True) + + +def fwht_and(sequence: list[int], inverse: bool = False) -> list[int]: + """ + Perform Fast Walsh-Hadamard Transform for AND operation. + + Time Complexity: O(N log N) + + >>> fwht_and([1, 2]) + [3, 2] + >>> fwht_and([3, 2], inverse=True) + [1, 2] + >>> fwht_and([1, 2, 3]) + Traceback (most recent call last): + ... + ValueError: Length of sequence must be a positive power of 2. + """ + sequence_length = len(sequence) + if sequence_length == 0 or (sequence_length & (sequence_length - 1)) != 0: + raise ValueError("Length of sequence must be a positive power of 2.") + + result = list(sequence) + half_block = 1 + while half_block < sequence_length: + block_size = half_block * 2 + for block_start in range(0, sequence_length, block_size): + for offset in range(half_block): + left_index = block_start + offset + right_index = left_index + half_block + if not inverse: + result[left_index] += result[right_index] + else: + result[left_index] -= result[right_index] + half_block *= 2 + + return result + + +def and_convolution(sequence_a: list[int], sequence_b: list[int]) -> list[int]: + """ + Compute bitwise AND convolution C[k] = sum_{i & j = k} (A[i] * B[j]). + + Time Complexity: O(N log N) + + >>> and_convolution([1, 2], [3, 4]) + [13, 8] + >>> and_convolution([1, 2], [3]) + Traceback (most recent call last): + ... + ValueError: Input sequences must have equal length. + """ + if len(sequence_a) != len(sequence_b): + raise ValueError("Input sequences must have equal length.") + + transformed_a = fwht_and(sequence_a) + transformed_b = fwht_and(sequence_b) + pointwise_product = [ + transformed_a[index] * transformed_b[index] for index in range(len(sequence_a)) + ] + return fwht_and(pointwise_product, inverse=True) + + +if __name__ == "__main__": + import doctest + + doctest.testmod()