negation & addition of Gaussian integers written in base i-1

def neg(n):
    '''Negate a Khmelnik-Penney-encoded Gaussian integer.

n is interpreted as a binary code with base i-1,
where i is the imaginary unit.
Algorithm from Solomon I. Khmelnik,
"Specialized Digital Computer for Operations with Complex Numbers"
(in Russian), Questions of Radio Electronics, volume 12, number 2, pp. 60-82, 1964,
http://lib.izdatelstwo.com/Papers2/s4.djvu
See table 4.

See also oeis.org/A066321, references therein, and the related SeqFan threads:
http://list.seqfan.eu/pipermail/seqfan/2016-August/016646.html
http://list.seqfan.eu/pipermail/seqfan/2016-August/016650.html
http://list.seqfan.eu/pipermail/seqfan/2016-September/016716.html
http://list.seqfan.eu/pipermail/seqfan/2016-October/016843.html'''
    sigma, shift = 0, 0
    while n >> shift:
        bits = n&(1<<shift)
        shift += 1
        sigma |= bits
        if bits:
            bits = n&(1<<shift)
            shift += 1
            sigma |= bits
            if bits:
                bits = ~n&(1<<shift)
                shift += 1
            else:
                bits = ~n&(7<<shift)
                shift += 3
            sigma |= bits
    return sigma

def neg2(n):
    '''Negate a Khmelnik-Penney-encoded Gaussian integer.

A shorter, more cryptic and slightly slower version of neg(n).'''
    n0, sigma, shift = n, 0, 0
    while n:
        if n&1:
            sigma |= (~n & 31>>(n&2) ^ 3) << shift
        shift += [1, 5, 1, 3][n%4]
        n = n0 >> shift
    return sigma

def add(a, b):
    '''Add up two Khmelnik-Penney-encoded Gaussian integers.

a and b are interpreted as binary codes with base i-1.
Algorithm from Solomon I. Khmelnik,
"Specialized Digital Computer for Operations with Complex Numbers"
(in Russian), Questions of Radio Electronics, volume 12, number 2, pp. 60-82, 1964,
http://lib.izdatelstwo.com/Papers2/s4.djvu
See table 7.'''
    sigma, shift, sk = 0, 0, 0o4444
    # Every octal digit of sk corresponds to one of the "retarded carries" P.
    while a or b or sk not in [0o4444, 0o4322]:
        sk = ((sk>>3) + [0o5000, 0o3670, 0o4000, 0o2670][(sk>>1)%4]
              + (a&1) + (b&1))
        sigma |= (sk&1) << shift
        shift += 1
        a >>= 1
        b >>= 1
    return sigma

def from_khmelnik(n, cx = False):
    '''Obtain a Khmelnik-Penney-encoded Gaussian integer.

Returns a pair of integers, i. e., real and imaginary
parts of the result. Alternatively, if the second argument
is set to True, returns a complex number.'''
    gre, gim, bre, bim = 0, 0, 1, 0
    while n:
        if n&1:
            gre += bre
            gim += bim
        n >>= 1
        bre, bim = -bre-bim, bre-bim
    return gre+1j*gim if cx else (gre, gim)

def to_khmelnik(re, im = None):
    '''Obtain Khmelnik-Penney encoding of a Gaussian integer.

Accepts either a complex number or a pair of integers
(real and imaginary parts).'''
    if im is None:
        if isinstance(re, int):
            im = 0
        else:
            re, im = round(re.real), round(re.imag)
    n = 0 # to_khmelnik(0)
    k = 1 # to_khmelnik(1)
    if re < 0:
        k = neg(k)
        re = -re
    while re:
        if re&1:
            n = add(n, k)
        re >>= 1
        k = add(k, k)
    k = 3 # to_khmelnik(1j)
    if im < 0:
        k = neg(k)
        im = -im
    while im:
        if im&1:
            n = add(n, k)
        im >>= 1
        k = add(k, k)
    return n