[WIP] Special Number Field Sieve (SageMath)

special-number-field-sieve.sage

P.<X> = ZZ[]

def select_polynomial(m, d, p):
    f = [0] * (d + 1)
    f[d] = p // m^d
    r = f[d] * m^d
    for i in range(d - 1, -1, -1):
        f[i] = (p - r) // m^i
        r = r + f[i] * m^i
    f = P(f)
    g = X - m
    print("# Polynomials are f = {0} and g = {1}.".format(P(f), g))
    return f, g

def create_rational_factor_bases(f, B, m):
    res = []
    for q in primes(B):
        res.append((m % q, q))
    return res

def create_algebraic_factor_bases(f, B, d):
    res = []
    for q in primes(B):
        roots = []
        for i in range(q):
            # print(f, i, q, f(i), f(i) % q)
            if f(i) % q == 0:
                res.append((i, q))
            if len(roots) >= d:
                break
        if len(roots) >= d:
            res.append(*roots)
    return res

def norm_N(a, b, f):
    res = 0
    if a == 0 and b == 0:
        res = 0
    if b == 0:
        res = a^f.degree()
    elif a == 0:
        res = b^f.degree() * f[0]
    else:
        i0 = 0
        i1 = f.degree()
        while i1 != -1:
            res = res + f[i0] * a^i0 * (-b)^i1
            i0 = i0 + 1
            i1 = i1 - 1
    return abs(res)

def is_Bsmooth(n, B):
    factors = list(factor(int(n)))
    if len(factors) != 0 and factors[-1][0] <= B: 
        return True, dict(factors)
    else:
        return False, dict(factors)

def is_useful_relation(a, b, rational_factor_bases, algebraic_factor_bases):
    isSatisfied1 = False
    isSatisfied2 = False
    res = []
    for root, prime in rational_factor_bases:
        if a == (-b * root) % p:
            isSatisfied1 = True
    if not isSatisfied1:
        return False
    for root, prime in algebraic_factor_bases:
        if a == (-b * root) % p:
            isSatisfied2 = True
    return (isSatisfied1 and isSatisfied2)

def create_vector1(x, B):
    vec = []
    vec.append(1 if x < 0 else 0)
    factors = dict(factor(x))
    for q in primes(B):
        vec.append(factors.get(q, 0))
    return vec

def create_vector2(a, b, nx, B, algebraic_factor_bases):
    vec = []
    factors = dict(factor(nx))

    for root, prime in algebraic_factor_bases:
        elem = 0
        factor_p = factors.get(prime, None)
        if factor_p:
            if (a % prime) == ((-b * root) % prime):
                elem = factor_p
        vec.append(elem)
    return vec

def compute_schirokauermap_exp(f, l):
    Q.<y> = GF(l)[]
    return lcm([l^i[0].degree() - 1 for i in Q(f).factor()])

def create_vector3(a, b, f, l, ε):
    m = matrix(IntegerModRing(l^2), [[a, -f[0]*b], [b, a-b]])
    sm = (m^ε) * vector([1,0]) - vector([1,0])
    res = []
    for i in range(2):
        res.append(int(sm[i]) // l)
    return res


l = 509  # p = 2*l + 1
p = 1019
q = 509
g = 277
y = 487
d = 2
m = int(p^(1/d))
B = 15

f1, f2 = select_polynomial(m, d, p)
rfb = create_rational_factor_bases(f1, B, m)
afb = create_algebraic_factor_bases(f1, B, d)

ε = compute_schirokauermap_exp(f1, l)

print(rfb)
print(afb)

A = []

max_row = len(rfb) + len(afb) + 2
print(max_row)

# for a in range(1, B*2):
#     for b in range(-B//2, B):
for a, b in [(1,-7), (1,1), (1,4), (3,-1), (4,1), (8,7), (9,-1), (9,1), (9,25), (26,-1), (27,-2), (29,2), (35,-1), (37,-17)]:
    if gcd(a, b) != 1:
        continue
    c1 = a + b*m
    c2 = norm_N(a, b, f1)
    if not (is_Bsmooth(c1, B) and is_Bsmooth(c2, B)):
        continue
    v1 = create_vector1(c1, B)
    v2 = create_vector2(a, b, c2, B, afb)
    v3 = create_vector3(a, b, f1, l, ε)
    if (all(e == 0 for e in v1) or sum(e for e in v1) <= 1):
        continue
    if (all(e == 0 for e in v2) or sum(e for e in v2) <= 1):
        continue
    print("useful!", a, b)
    A.append(v1 + v2 + v3)
    print(v1, v2, v3)
    print('---')
    if len(A) >= max_row:
        break

# v_H
print('v_H')
while True:
    R = randint(2, p-1)
    # R = 187
    H = int(pow(g, R, p) * y) % p
    if is_Bsmooth(H, B):
        v1 = create_vector1(H, B)
        v2 = [0] * len(rfb)
        v3 = [0] * 2
        vH = vector(v1 + v2 + v3)
        if (all(e == 0 for e in vH) or sum(e for e in vH) <= 3):
            continue
        break

# v_G
print('v_G')
while True:
    S = randint(2, p-1)
    # S = 299
    G = int(pow(g, S, p))
    if is_Bsmooth(G, B):
        v1 = create_vector1(G, B)
        v2 = [0] * len(rfb)
        v3 = [0] * 2
        vG = vector(v1 + v2 + v3)
        if (all(e == 0 for e in vG) or sum(e for e in vG) <= 3):
            continue
        break

print(R, S)

print(vG)
print('---')

A.insert(0, vG)
A = matrix(GF(l), A).transpose()

vH = vector(GF(l), -vH)

print(A)
print(vH)

X = A.solve_right(vH)
print(X)

x = (-X[0] * S - R) % q
print(x)