th. baruchel's blog

hacks & code

Code snippets for the Binomial modulo 2 transform and its inverse

Posted at — Jan 11, 2020

While working on my latest paper, I encountered the need to compute the Binomial modulo 2 transform (as well as its inverse). Though the initial definition (which can be found in a comment of sequence A100735) obviously is in $O(n^2)$, it is rather easy to figure out some algorithm for computing it in $O(\log n)$.

The code assumes that the length of the list to be computed is a power of 2.

Lisp code

Here is my code for the Binomial modulo 2 transform:

; Binomial Modulo 2 Transform of a sequence (length is a power of 2)
(defun trvbin2 (v)
  (loop
    with a = (coerce v 'vector)
    with n = (car (array-dimensions a))
    for s = 2 then (ash s 1)
      and z = 1 then s
    when (= z n) return (coerce a 'list)
    do (loop
         for x from 0 to (1- n) by s
         do (loop
              for y from x to (+ x z -1)
              do (setf (aref a (+ y z)) (+ (aref a (+ y z)) (aref a y)))))))

The inverse transform is:

; Binomial Modulo 2 Inverse Transform of a sequence (length is a power of 2)
(defun itrvbin2 (v)
  (loop
    with a = (coerce v 'vector)
    with n = (car (array-dimensions a))
    for s = 2 then (ash s 1)
      and z = 1 then s
    when (= z n) return (coerce a 'list)
    do (loop
         for x from 0 to (1- n) by s
         do (loop
              for y from x to (+ x z -1)
              do (setf (aref a (+ y z)) (- (aref a (+ y z)) (aref a y)))))))

Scheme code

Here is my code for the Binomial modulo 2 transform:

; Binomial Modulo 2 Transform of a sequence (length is a power of 2)
(define trvbin2
  (lambda (v)
    (let* ((a (list->vector v))
           (n (vector-length a)))
      (do ((z 1 (do ((s (* 2 z)) (x 0 (+ x s)))
                  ((= x n) s)
                   (do ((y x (+ y 1)) (t (+ x z)))
                     ((= y t))
                     (vector-set! a
                                  (+ y z)
                                  (+ (vector-ref a (+ y z))
                                     (vector-ref a y)))))))
        ((= z n) (vector->list a))))))

The inverse transform is:

; Binomial Modulo 2 Inverse Transform of a sequence (length is a power of 2)
(define itrvbin2
  (lambda (v)
    (let* ((a (list->vector v))
           (n (vector-length a)))
      (do ((z 1 (do ((s (* 2 z)) (x 0 (+ x s)))
                  ((= x n) s)
                   (do ((y x (+ y 1)) (t (+ x z)))
                     ((= y t))
                     (vector-set! a
                                  (+ y z)
                                  (- (vector-ref a (+ y z))
                                     (vector-ref a y)))))))
        ((= z n) (vector->list a))))))

Python code

The main transform is:

import numpy as np

def trbin(v):
    t, n = 1, len(v)
    v = np.array(v, dtype=object).reshape((n, t))
    while n != 1:
        v[1::2,:] += v[::2,:]
        t <<= 1; n >>= 1
        v = v.reshape((n, t))
    return list(v[0])

The inverse transform is:

def trbin_inv(v):
    t, n = 1, len(v)
    v = np.array(v, dtype=object).reshape((n, t))
    while n != 1:
        v[1::2,:] -= v[::2,:]
        t <<= 1; n >>= 1
        v = v.reshape((n, t))
    return list(v[0])