Recover the truncated MRG-with-c driving sequence

'''
This code is used to recover the truncated MRG-with-c's driving sequence.
Suppose that the modulus, the coefficients, the initial states and the constant are unknown, given the first N truncated digits which consist of the alpha high-order bits of the sequence, we need to recover the modulus, the coefficients, the initial states and the constant.
The library used is the Sagemath.
The whole experimental process consists of 5 steps. Searching linear relations, Recovering the modulus, the coefficients, the differences of initial states, the initial states and the constant.
'''

import time

# Search polynomials that annihilate the sequence a'
def Search_linear_relations(y, r, t):
    time_start = time.clock()
    L = Matrix(ZZ, r, r + t)
    C = Matrix(ZZ, 5, r)
    for i in range(0, r):
        for j in range(0, r + t):
            L[i, j] = 0
        L[i, i + t] = 1
        for j in range(0, t):
            L[i, j] = y[i + j]
    reduced_L = L.BKZ(block_size = 20)

    for i in range(0, 5):
        for j in range(0, r):
            C[i, j] = reduced_L[i, j + t]

    time_end = time.clock()
    print('Step 1: time of searching polynomials that annihilate the sequence a`: ', time_end - time_start, 's\n\n')
    return C

# Recover the modulus by the resultant
def Recover_the_modulus(C, r):
    time_start = time.clock()
    R.<x> = ZZ[]
    f1 = C[0, 0]
    f2 = C[1, 0]
    f3 = C[2, 0]
    f4 = C[3, 0]
    f5 = C[4, 0]
    for i in range(1, r):
        f1 = f1 + x^i * C[0, i]
        f2 = f2 + x^i * C[1, i]
        f3 = f3 + x^i * C[2, i]
        f4 = f4 + x^i * C[3, i]
        f5 = f5 + x^i * C[4, i]

    r1 = f1.resultant(f2)
    r2 = f2.resultant(f3)
    r3 = f3.resultant(f1)
    r4 = f4.resultant(f1)
    r5 = f5.resultant(f2)
    r6 = f4.resultant(f3)

    m_power = gcd(gcd(gcd(gcd(gcd(r1, r2), r3), r4), r5), r6)
    m = m_power^(1/16)

    time_end = time.clock()
    print('Step 2: time of recovering the modulus: ', time_end - time_start, 's\n\n')
    return m

# Recover the coefficients by the Chinese Remainder Theorem
def Recover_the_coefficients(C, m, r):
    time_start = time.clock()
    R.<x> = PolynomialRing(Integers(m), implementation='NTL')

    for i in range(0, 2):
        mark = r - 1
        while C[i][mark] == 0:
            mark = mark - 1
        for j in range(0, r):
            C[i, j] = C[i, j] * ((C[i, mark]^(-1)) % m)
    f = C[0, 0]
    g = C[1, 0]
    for i in range(1, r):
        f = f + x^i * C[0, i]
        g = g + x^i * C[1, i]

    time_end = time.clock()
    print('Step 3: time of recovering the coefficients: ', time_end - time_start, 's\n\n')
    return f.gcd(g)

# Recover the differences of initial state by Kannan's embedding technique
def Recover_the_differences_of_initial_state(y, f, m, n, beta):
    time_start = time.clock()
    d = 31
    z = vector(ZZ, n)
    a = vector(ZZ, n)
    Q = Matrix(ZZ, n, n)
    Q_power = identity_matrix(ZZ, n)
    L = Matrix(ZZ, d + 1, d + 1)

    Q[0, n - 1] = (-f.list()[0]) % m
    for i in range(1, n):
        Q[i, i - 1] = 1
        Q[i, n - 1] = (-f.list()[i]) % m

    for i in range(1, n):
        Q_power = (Q_power * Q) % m
    for i in range(n, d):
        Q_power = (Q_power * Q) % m
        b = 0
        for j in range(0, n):
            L[j + 1, i + 1] = Q_power[j, 0]
            b = b + Q_power[j, 0] * y[j]
        L[0, i + 1] = (2^beta * (y[i] - b)) % m

    L[0, 0] = 2^beta
    for i in range(1, n + 1):
        L[i, i] = 1
    for i in range(n + 1, d + 1):
        L[i, i] = m

    reduced_L = L.BKZ(block_size = 20)

    if reduced_L[0, 0] == -2^beta:
        for i in range(0, n):
            z[i] = reduced_L[0, i + 1]
            a[i] = y[i] * 2^(beta) + z[i]
    if reduced_L[0, 0] == 2^beta:
        for i in range(0, n):
            z[i] = -reduced_L[0, i + 1]
            a[i] = y[i] * 2^(beta) + z[i]

    time_end = time.clock()
    print('Step 4: time of recovering the differences of initial state: ', time_end - time_start, 's\n\n')
    return a

# Recover the constant and initial state
def Recover_the_constant_and_initial_state(dif_a, y, beta, n, m, f):
    time_start = time.clock()
    diff_a = vector(ZZ, 9999)
    diff_y = vector(ZZ, 9999)

    for i in range(0, 9999):
        diff_y[i] = y[i + 1] - y[i]
    for i in range(0, n):
        diff_a[i] = dif_a[i]
    for i in range(n, 9999):
        diff_a[i] = 0
        for j in range(1, n + 1):
            diff_a[i] = (diff_a[i] + diff_a[i - j] * (-f.list()[n - j])) % m
        if diff_y[i] < 0:
            diff_a[i] = diff_a[i] - m

    z0_left = 0
    z0_right = 2^beta - 1
    L = 1
    Y = 0
    A = 0
    while z0_left!= z0_right:
        L = L + 1
        Y = Y + diff_y[L - 2]
        A = A + diff_a[L - 2]
        if (A >> beta) == Y:
            right = 2^beta - 1 - A % (2^beta)
            if right < z0_right:
                z0_right = right
        if (A >> beta) + 1 == Y:
            left = 2^beta - A % (2^beta)
            if left > z0_left:
                z0_left = left
    z0 = z0_left

    a = vector(ZZ, n + 1)
    a[0] = z0 + 2^beta * y[0]
    for i in range(1, n + 1):
        a[i] = a[i - 1] + dif_a[i - 1]

    c = 0
    for i in range(0, n + 1):
        c = (c + a[i] * f.list()[i]) % m

    time_end = time.clock()
    print('Step 5: time of recovering the constant and initial state: ', time_end - time_start, 's\n\n''The modulus m = ', m, '\n\n')
f = Recover_the_coefficients(C, m, r)
print('The characteristic polynomial f = ', f, '\n\n')
dif_a = Recover_the_differences_of_initial_state(dif_y, f, m, n, beta)
print('The differences of initial state a` = ', dif_a, '\n\n')
L, a, c = Recover_the_constant_and_initial_state(dif_a, y, beta, n, m, f)
print('The least number of truncated digits L = ', L, '\n\n')
print('The initial state a = ', a, '\n\n')
print('The constant c = ', c, '\n\n')