#!/usr/bin/env python3

R =  0xE1000000000000000000000000000000

# I often use reverseBits() here.
# rationale: identity element would be represented as decimal number 1, not as 0x80....00

# yes, both are XORs:
def gf2_add(x, y):
    return x^y

def gf2_sub(x, y):
    return x^y

def gf128_mul(x, y, R):
    z = 0
    for i in range(128-1, -1, -1):

        if (y >> i) & 1:
            z=gf2_add(z, x)

        # shift and also reduce by R if overflow detected
        # IOW, keep x smaller than R or modulo R
        if x & 1:
            x = gf2_sub(x >> 1, R)
        else:
            x = x >> 1

    return z

def gf2_pow_2(x, R):
    return gf128_mul(x, x, R)

def is_odd(n):
    return n&1==1

# almost as in \url{https://en.wikipedia.org/wiki/Exponentiation_by_squaring}
def gf2_pow(x, n, R):
    if n==1:
         return x
    if is_odd(n):
        return gf128_mul(x, gf2_pow(gf2_pow_2(x, R), (n-1)//2, R), R)
    else:
        return gf2_pow(gf2_pow_2(x, R), n//2, R)

# return $x^{2^{128}-2}$
# AKA reciprocal AKA modulo inverse
def gf2_inv(x, R):
    rslt=reverseBits(1, 128) # init to 1

    for i in range(128-1):
        rslt=gf128_mul(rslt, x, R)
        x=gf2_pow_2(x, R)

    return gf2_pow_2(rslt, R)

# simpler version:
def gf2_inv_v2(x, R):
    return gf2_pow(x, 2**128-2, R)

# https://math.stackexchange.com/questions/943417/square-root-for-galois-fields-gf2m
# https://crypto.stackexchange.com/questions/17988/algorithm-for-computing-square-roots-in-gf2n
# return $x^{2^{127}}$
def gf2_sqrt(x, R):
    for i in range(127):
        x=gf2_pow_2(x, R)

    return x

# simpler version:
def gf2_sqrt_v2(x, R):
    return gf2_pow(x, 2**127, R)

def gf128_div(N, D, R):
    # TODO: check for zero
    # N = numerator (dividend)
    # D = denominator (divisor)
    i1=gf2_inv(D, R)
    # also, test 2nd version:
    i2=gf2_inv_v2(D, R)
    assert i1==i2

    return gf128_mul(N, i1, R)