Factoring Polynomials Mod 2

main.py

import random


class GF2Poly:
    def __init__(self, c):
        self.c = c

    def __add__(self, other):
        return GF2Poly(other.c ^ self.c)

    def __mul__(self, other):
        temp = self.c
        addend = other.c
        result = 0
        while temp != 0:
            if temp % 2 == 1:
                result ^= addend
            temp >>= 1
            addend <<= 1
        return GF2Poly(result)

    def __str__(self):
        return bin(self.c)

    def __mod__(self, other):
        ol = other.deg()
        shift = self.deg() - ol
        result = GF2Poly(self.c)
        while shift >= 0:
            result += GF2Poly(other.c << shift)
            shift = result.deg() - ol
        return result

    def __truediv__(self, other):
        ol = other.deg()
        shift = self.deg() - ol
        remainder = GF2Poly(self.c)
        result = 0
        while shift >= 0:
            result += 1 << shift
            remainder += GF2Poly(other.c << shift)
            shift = remainder.deg() - ol
        return GF2Poly(result)

    def deg(self):
        return self.c.bit_length() - 1

    def derivative(self):
        return GF2Poly((self.c >> 1) & (((1 << (2 * (self.deg() // 2) + 2)) - 1) // 3))


def gcd(a, b):
    if a.c < b.c:
        a, b = b, a
    while b.c:
        shift = a.deg() - b.deg()
        a, b = b, (a + GF2Poly(b.c << shift))
        if a.c < b.c:
            a, b = b, a
    return a


def ddf(f):
    d = f.deg()
    # print(d)
    results = []
    t = GF2Poly(0b10)
    for i in range(1, (d // 2) + 1):
        # compute x^(2^i) + x mod f
        t = (t * t) % f
        other = t + GF2Poly(0b10)
        r = gcd(f, other)
        if r.c != 1:
            # print(f"{f} {r}")
            results.append((r, i))
        f /= r
    if f.c != 1:
        # print(f"{f} {f.deg()}")
        results.append((f, f.deg()))
    return results


def cantor_zassenhaus(f, d):
    if f.deg() == d:
        return
    while True:
        trace = GF2Poly(1)
        b = GF2Poly(random.randint(2, f.c))
        for i in range(0, d):
            trace += b
            b = (b * b) % f
        if f.deg() > gcd(f, trace).deg() > 0:
            return gcd(f, trace)


def cz_recurse(f, d, lst):
    if f.deg() == d:
        # print(f)
        lst.append(f)
        return
    res = cantor_zassenhaus(f, d)
    cz_recurse(f / res, d, lst)
    cz_recurse(res, d, lst)


def sqrt(f):
    old = f.c
    result = 0
    i = 0
    while old != 0:
        result += (old & 1) << i
        i += 1
        old >>= 2
    return GF2Poly(result)


def factor(f, lst):
    g = gcd(f, f.derivative())
    f /= g
    if g.c != 1:
        tmp = []
        factor(sqrt(g), tmp)
        lst += tmp
        lst += tmp
    for t in ddf(f):
        # print(f"Factoring {f} {t[0]} {t[0].deg()} {t[1]}")
        cz_recurse(t[0], t[1], lst)


if __name__ == "__main__":
    f = GF2Poly(0xb0c152f9) * GF2Poly(0xebf2831f)
    lst = []
    factor(f, lst)
    res = GF2Poly(0b1)
    for i in lst:
        res *= i
    print(f"Expect: {f}")
    print(f"Result: {res}")