An Improved Method for Predicting Truncated Multiple Recursive Generators

// This code is used to recover the truncated ZUC's driving sequence.
// The library used is the NTL 11.4.3 library.
// The whole experimental process consists of 4 steps. Searching linear relations, Recovering the modulus, the coefficients and the initial states.

#include <NTL/vec_ZZ.h>
#include <NTL/mat_ZZ.h>
#include <NTL/mat_ZZ_p.h>
#include <NTL/ZZ_pX.h>
#include <NTL/ZZX.h>
#include <NTL/LLL.h>
#include <fstream>
#include <ctime>
#include <math.h>
#include <fstream>
NTL_CLIENT
using namespace std;

// Search polynomials that annihilate the sequence a
// The high-order scenario
// The modulus m is known
mat_ZZ Search_linear_relations_high_m(ZZ m, vec_ZZ y, int beta, int r, int t)
{
    clock_t start, end;
    start = clock();
    mat_ZZ L, C;
    int i, j;

    L.SetDims(r + t, r + t);
    C.SetDims(r, r);

    // Construct the lattice L
    for (i = 0; i < r + t; i++)
    {
        for (j = 0; j < r + t; j++)
        {
            L[i][j] = 0;
        }
    }
    for (i = 0; i < t; i++)
        L[i][i] = m;
    for (i = t; i < r + t; i++)
    {
        L[i][i] = power2_ZZ(beta);
        for (j = 0; j < t; j++)
        {
            L[i][j] = y[i + j - t] * power2_ZZ(beta);
        }
    }

    // Reduce the lattice L
    BKZ_FP(L, 0.99, 20);

    // Collect annihilating polynomials of the sequence a
    for (i = 0; i < r; i++)
    {
        for (j = 0; j < r; j++)
        {
            C[i][j] = L[i][j + t] / power2_ZZ(beta);
        }
    }

    end = clock();
    double time = (double)(end - start) / CLK_TCK;
    cout << "Step 1: time of searching polynomials that annihilate the sequence a: " << time << "s\n\n\n";

    return C;
}

// Search polynomials that annihilate the sequence a
// The high-order scenario
// The modulus m is unknown
mat_ZZ Search_linear_relations_high(vec_ZZ y, int alpha, int r, int t)
{
    clock_t start, end;
    start = clock();
    mat_ZZ L, C;
    int i, j;

    L.SetDims(r + t, r + t);
    C.SetDims(r, r);

    // Construct the lattice L
    for (i = 0; i < r + t; i++)
    {
        for (j = 0; j < r + t; j++)
        {
            L[i][j] = 0;
        }
    }
    for (i = 0; i < t; i++)
        L[i][i] = power2_ZZ(alpha);
    for (i = t; i < r + t; i++)
    {
        L[i][i] = 1;
        for (j = 0; j < t; j++)
        {
            L[i][j] = y[i + j - t];
        }
    }

    // Reduce the lattice L
    BKZ_FP(L, 0.99, 20);

    // Collect annihilating polynomials of the sequence a
    for (i = 0; i < r; i++)
    {
        for (j = 0; j < r; j++)
        {
            C[i][j] = L[i][j + t];
        }
    }

    end = clock();
    double time = (double)(end - start) / CLK_TCK;
    cout << "Step 1: time of searching polynomials that annihilate the sequence a: " << time << "s\n\n\n";

    return C;
}

// Search polynomials that annihilate the sequence a
// The multi-segment scenario, l = 3
// The modulus m is unknown
mat_ZZ Search_linear_relations_multi_segment_3(vec_ZZ a2, int k2, int r, int t)
{
    clock_t start, end;
    start = clock();
    mat_ZZ L, C;
    int i, j;

    L.SetDims(r + t, r + t);
    C.SetDims(r, r);

    // Construct the lattice L
    for (i = 0; i < r + t; i++)
    {
        for (j = 0; j < r + t; j++)
        {
            L[i][j] = 0;
        }
    }
    for (i = 0; i < t; i++)
        L[i][i] = power2_ZZ(k2);
    for (i = t; i < r + t; i++)
    {
        L[i][i] = 1;
        for (j = 0; j < t; j++)
        {
            L[i][j] = a2[i + j - t];
        }
    }

    // Reduce the lattice L
    BKZ_FP(L, 0.99, 20);

    // Collect annihilating polynomials of the sequence a
    for (i = 0; i < r; i++)
    {
        for (j = 0; j < r; j++)
        {
            C[i][j] = L[i][j + t];
        }
    }

    end = clock();
    double time = (double)(end - start) / CLK_TCK;
    cout << "Step 1: time of searching polynomials that annihilate the sequence a: " << time << "s\n\n\n";

    return C;
}

// Search polynomials that annihilate the sequence a
// The multi-segment scenario, l = 5
// The modulus m is unknown
mat_ZZ Search_linear_relations_multi_segment_5(vec_ZZ a2, vec_ZZ a4, int v2, int v4, int v5, int r, int t)
{
    clock_t start, end;
    start = clock();
    mat_ZZ L, C;
    int i, j;

    L.SetDims(r + t, r + t);
    C.SetDims(r, r);

    // Construct the lattice L
    for (i = 0; i < r + t; i++)
    {
        for (j = 0; j < r + t; j++)
        {
            L[i][j] = 0;
        }
    }
    for (i = 0; i < t; i++)
        L[i][i] = power2_ZZ(v5 - v2);
    for (i = t; i < r + t; i++)
    {
        L[i][i] = power2_ZZ(v4 - v2);
        for (j = 0; j < t; j++)
        {
            L[i][j] = a2[i + j - t] + a4[i + j - t] * power2_ZZ(v4 - v2);
        }
    }

    // Reduce the lattice L
    BKZ_FP(L, 0.99, 20);

    // Collect annihilating polynomials of the sequence a
    for (i = 0; i < r; i++)
    {
        for (j = 0; j < r; j++)
        {
            C[i][j] = L[i][j + t] / power2_ZZ(v4 - v2);
        }
    }

    end = clock();
    double time = (double)(end - start) / CLK_TCK;
    cout << "Step 1: time of searching polynomials that annihilate the sequence a: " << time << "s\n\n\n";

    return C;
}

// Recover the modulus by the resultant
ZZ Recover_the_modulus(mat_ZZ C, int n, int k, int r)
{
    clock_t start, end;
    start = clock();
    ZZ res, m_power, recovered_m;
    ZZX f, g;
    int i, j, flag;

    f.SetLength(r);
    g.SetLength(r);

    flag = 0;
    for (i = 0; i < r; i++)
        if (C[0][i] == 0)
            flag = flag + 1;
    if (flag == r)
        cout << "Fail to recover the modulus m!\n";

    flag = 2;
    m_power = 0;
    while (flag < r)
    {
        for (i = 0; i < r; i++)
            SetCoeff(f, i, C[flag - 1][i]);
        for (j = 0; j < flag - 1; j++)
        {
            for (i = 0; i < r; i++)
            {
                SetCoeff(g, i, C[j][i]);
            }
            res = resultant(f, g);
            m_power = GCD(m_power, res);
            //cout << m_power << "\n";
        }
        if (NumBits(m_power) == k * n)
            break;
        flag = flag + 1;
    }
    recovered_m = SqrRoot(SqrRoot(SqrRoot(SqrRoot(m_power))));

    end = clock();
    double time = (double)(end - start) / CLK_TCK;
    cout << "Step 2: time of recovering the modulus: " << time << "s\n";
    cout << "Number of the annihilating polynomials for recovering the modulus = " << flag << "\n";

    return recovered_m;
}

// Recover the coefficients by the Chinese Remainder Theorem
ZZ_pX Recover_the_coefficients(mat_ZZ C, ZZ m, int n, int r)
{
    clock_t start, end;
    start = clock();
    ZZ_pX f, g, h;
    mat_ZZ_p C_m;
    ZZ_p c;
    int i, j, mark, flag;
    c.init(m);

    C_m.SetDims(r, r);

    ofstream fout("C.txt");
    for (i = 0; i < r; i++)
    {
        for (j = 0; j < r; j++)
        {
            fout << C[i][j] << " ";
        }
        fout << "\n";
    }
    fout.close();
    freopen("C.txt", "r", stdin);

    // Make the polynomials monic
    for (i = 0; i < r; i++)
    {
        for (j = 0; j < r; j++)
        {
            cin >> c;
            C_m[i][j] = c;
        }
        mark = r - 1;
        while (C_m[i][mark] == 0)
        {
            mark = mark - 1;
        }
        if (mark >= 0)
        {
            for (j = 0; j < r; j++)
            {
                C_m[i][j] = C_m[i][j] * inv(C_m[i][mark]);
            }
        }
    }

    f.SetLength(r);
    g.SetLength(r);

    for (i = 0; i < r; i++)
    {
        SetCoeff(f, i, C_m[0][i]);
        SetCoeff(g, i, C_m[1][i]);
    }

    h = GCD(f, g);
    flag = 2;

    while (deg(h) > n)
    {
        for (i = 0; i < r; i++)
            SetCoeff(f, i, C_m[flag][i]);
        h = GCD(f, h);
        flag = flag + 1;
    }
    end = clock();
    double time = (double)(end - start) / CLK_TCK;
    cout << "Step 3: time of recovering the coefficients: " << time << "s\n";
    cout << "Number of the annihilating polynomials for recovering the coefficients = " << flag << "\n";

    return h;
}

// Recover the initial state by Kannan's embedding technique
// The high-order scenario
vec_ZZ Recover_the_initial_state_high(vec_ZZ y, ZZ_pX f, ZZ m, int n, int d, int beta)
{
    clock_t start, end;
    start = clock();
    vec_ZZ a, z;
    mat_ZZ Q, Q_power, L;
    ZZ b;
    int i, j, l;

    a.SetLength(n);
    z.SetLength(n);
    Q.SetDims(n, n);
    Q_power.SetDims(n, n);
    L.SetDims(d + 1, d + 1);

    Q_power = ident_mat_ZZ(n);

    Q[0][n - 1] = (-rep(f[0])) % m;
    for (i = 1; i < n; i++)
    {
        Q[i][i - 1] = 1;
        Q[i][n - 1] = (-rep(f[i])) % m;
    }

    // Construct the lattice L
    for (i = 1; i < n; i++)
    {
        Q_power = Q_power * Q;
        for (j = 0; j < n; j++)
        {
            for (l = 0; l < n; l++)
            {
                Q_power[j][l] = Q_power[j][l] % m;
            }
        }
    }

    for (i = n; i < d; i++)
    {
        Q_power = Q_power * Q;
        for (j = 0; j < n; j++)
        {
            for (l = 0; l < n; l++)
            {
                Q_power[j][l] = Q_power[j][l] % m;
            }
        }
        b = 0;
        for (j = 0; j < n; j++)
        {
            L[j + 1][i + 1] = Q_power[j][0];
            b += Q_power[j][0] * y[j];
        }
        L[0][i + 1] = (pow(2, beta) * (y[i] - b)) % m + pow(2, beta - 1);
    }

    L[0][0] = pow(2, beta - 1);
    for (i = 1; i <= n; i++)
    {
        L[0][i] = pow(2, beta - 1);
        L[i][i] = 1;
    }
    for (i = n + 1; i <= d; i++)
    {
        L[i][i] = m;
    }

    // Reduce the lattice L
    BKZ_FP(L, 0.99, 20);

    // Recover the initial state
    if (L[0][0] == -pow(2, beta - 1))
    {
        for (i = 0; i < n; i++)
        {
            z[i] = L[0][i + 1] + pow(2, beta - 1);
            a[i] = y[i] * pow(2, beta) + z[i];
        }
    }
    if (L[0][0] == pow(2, beta - 1))
    {
        for (i = 0; i < n; i++)
        {
            z[i] = pow(2, beta - 1) - L[0][i + 1];
            a[i] = y[i] * pow(2, beta) + z[i];
        }
    }

    end = clock();
    double time = (double)(end - start) / CLK_TCK;
    cout << "Step 4: time of recovering the initial state: " << time << "s\n";

    return a;
}

// Recover the initial state by Kannan's embedding technique
// The multi-segment scenario, l = 3
vec_ZZ Recover_the_initial_state_multi_segment_3(vec_ZZ a2, ZZ_pX f, ZZ m, int n, int d, int k1, int k2, int k3)
{
    clock_t start, end;
    start = clock();
    vec_ZZ a, a1, a3;
    mat_ZZ Q, Q_power, L;
    ZZ b;
    int i, j, jj;

    a.SetLength(n);
    a1.SetLength(n);
    a3.SetLength(n);
    Q.SetDims(n, n);
    Q_power.SetDims(n, n);
    L.SetDims(2 * d + 1, 2 * d + 1);

    Q_power = ident_mat_ZZ(n);

    Q[0][n - 1] = (-rep(f[0])) % m;
    for (i = 1; i < n; i++)
    {
        Q[i][i - 1] = 1;
        Q[i][n - 1] = (-rep(f[i])) % m;
    }

    // Construct the lattice L
    for (i = 1; i < n; i++)
    {
        Q_power = Q_power * Q;
        for (j = 0; j < n; j++)
        {
            for (jj = 0; jj < n; jj++)
            {
                Q_power[j][jj] = Q_power[j][jj] % m;
            }
        }
    }

    for (i = n; i < d; i++)
    {
        Q_power = Q_power * Q;
        for (j = 0; j < n; j++)
        {
            for (jj = 0; jj < n; jj++)
            {
                Q_power[j][jj] = Q_power[j][jj] % m;
            }
        }
        b = 0;
        for (j = 0; j < n; j++)
        {
            L[j + 1][i + 1] = Q_power[j][0];
            L[j + d + 1][i + 1] = pow(2, k1 + k2) * Q_power[j][0];
            b += Q_power[j][0] * a2[j];
        }
        L[0][i + 1] = (pow(2, k1) * (a2[i] - b)) % m + pow(2, k1 - 1);
    }

    L[0][0] = pow(2, k1 - 1);
    for (i = 1; i <= n; i++)
    {
        L[0][i] = pow(2, k1 - 1);
        L[i][i] = 1;
    }
    for (i = n + 1; i <= d; i++)
    {
        L[i][i] = m;
        L[d + i][i] = -pow(2, k1 + k2);
    }
    for (i = d + 1; i <= 2 * d; i++)
        L[0][i] = pow(2, k1 - 1);

    if (k1 >= k3)
    {
        for (i = d + 1; i <= 2 * d; i++)
            L[i][i] = pow(2, k1 - k3);
    }
    else
    {
        for (i = 0; i <= 2 * d; i++)
            for (j = 0; j <= 2 * d; j++)
                L[i][j] = L[i][j] * pow(2, k3 - k1);
        for (i = d + 1; i <= 2 * d; i++)
            L[i][i] = 1;
    }

    // Reduce the lattice L
    BKZ_FP(L, 0.99, 20);

    // Recover the initial state
    if (k1 >= k3)
    {
        j = 0;
        while (L[j][0] != pow(2, k1 - 1) && L[j][0] != -pow(2, k1 - 1))
            j = j + 1;
        if (L[j][0] == -pow(2, k1 - 1))
        {
            for (i = 0; i < n; i++)
            {
                a1[i] = L[j][i + 1] + pow(2, k1 - 1);
                a3[i] = (L[j][i + d + 1] + pow(2, k1 - 1)) / pow(2, k1 - k3);
                a[i] = (a3[i] * pow(2, k1 + k2) + a2[i] * pow(2, k1) + a1[i]) % m;
            }
        }
        if (L[j][0] == pow(2, k1 - 1))
        {
            for (i = 0; i < n; i++)
            {
                a1[i] = -L[j][i + 1] + pow(2, k1 - 1);
                a3[i] = (-L[j][i + d + 1] + pow(2, k1 - 1)) / pow(2, k1 - k3);
                a[i] = (a3[i] * pow(2, k1 + k2) + a2[i] * pow(2, k1) + a1[i]) % m;
            }
        }
    }
    if (k1 < k3)
    {
        j = 0;
        while (L[j][0] != pow(2, k3 - 1) && L[j][0] != -pow(2, k3 - 1))
            j = j + 1;
        if (L[j][0] == -pow(2, k3 - 1))
        {
            for (i = 0; i < n; i++)
            {
                a1[i] = (L[j][i + 1] + pow(2, k3 - 1)) / pow(2, k3 - k1);
                a3[i] = L[j][i + d + 1] + pow(2, k3 - 1);
                a[i] = (a3[i] * pow(2, k1 + k2) + a2[i] * pow(2, k1) + a1[i]) % m;
            }
        }
        if (L[j][0] == pow(2, k3 - 1))
        {
            for (i = 0; i < n; i++)
            {
                a1[i] = (-L[j][i + 1] + pow(2, k3 - 1)) / pow(2, k3 - k1);
                a3[i] = -L[j][i + d + 1] + pow(2, k3 - 1);
                a[i] = (a3[i] * pow(2, k1 + k2) + a2[i] * pow(2, k1) + a1[i]) % m;
            }
        }
    }

    end = clock();
    double time = (double)(end - start) / CLK_TCK;
    cout << "Step 4: time of recovering the initial state: " << time << "s\n";

    return a;
}

// Recover the initial state by Kannan's embedding technique
// The multi-segment scenario, l = 5
vec_ZZ Recover_the_initial_state_multi_segment_5(vec_ZZ a2, vec_ZZ a4, ZZ_pX f, ZZ m, int n, int d, int k1, int k2, int k3, int k4, int k5)
{
    clock_t start, end;
    start = clock();
    vec_ZZ a, a1, a3, a5;
    mat_ZZ Q, Q_power, L;
    ZZ b;
    int i, j, jj;

    a.SetLength(n);
    a1.SetLength(n);
    a3.SetLength(n);
    a5.SetLength(n);
    Q.SetDims(n, n);
    Q_power.SetDims(n, n);
    L.SetDims(3 * d + 1, 3 * d + 1);

    Q_power = ident_mat_ZZ(n);

    Q[0][n - 1] = (-rep(f[0])) % m;
    for (i = 1; i < n; i++)
    {
        Q[i][i - 1] = 1;
        Q[i][n - 1] = (-rep(f[i])) % m;
    }

    // Construct the lattice L
    for (i = 1; i < n; i++)
    {
        Q_power = Q_power * Q;
        for (j = 0; j < n; j++)
        {
            for (jj = 0; jj < n; jj++)
            {
                Q_power[j][jj] = Q_power[j][jj] % m;
            }
        }
    }

    for (i = n; i < d; i++)
    {
        Q_power = Q_power * Q;
        for (j = 0; j < n; j++)
        {
            for (jj = 0; jj < n; jj++)
            {
                Q_power[j][jj] = Q_power[j][jj] % m;
            }
        }
        b = 0;
        for (j = 0; j < n; j++)
        {
            L[j + 1][i + 1] = Q_power[j][0];
            L[j + d + 1][i + 1] = pow(2, k1 + k2) * Q_power[j][0];
            L[j + d + n + 1][i + 1] = pow(2, k1 + k2 + k3 + k4) * Q_power[j][0];
            b += Q_power[j][0] * (pow(2, k1) * a2[j] + pow(2, k1 + k2 + k3) * a4[j]);
        }
        L[0][i + 1] = (pow(2, k1) * a2[i] + pow(2, k1 + k2 + k3) * a4[i] - b) % m + pow(2, k1 - 1);
    }

    L[0][0] = pow(2, k1 - 1);
    for (i = 1; i <= n; i++)
    {
        L[0][i] = pow(2, k1 - 1);
        L[i][i] = 1;
    }
    for (i = n + 1; i <= d; i++)
    {
        L[i][i] = m;
        L[d + n + i][i] = -pow(2, k1 + k2);
        L[2 * d + i][i] = -pow(2, k1 + k2 + k3 + k4);
    }
    for (i = d + 1; i <= 3 * d; i++)
        L[0][i] = pow(2, k1 - 1);

    if (k1 >= k3 && k1 >= k5)
    {
        for (i = d + 1; i <= d + n; i++)
            L[i][i] = pow(2, k1 - k3);
        for (i = d + n + 1; i <= d + 2 * n; i++)
            L[i][i] = pow(2, k1 - k5);
        for (i = d + 2 * n + 1; i <= 2 * d + n; i++)
            L[i][i] = pow(2, k1 - k3);
        for (i = 2 * d + n + 1; i <= 3 * d; i++)
            L[i][i] = pow(2, k1 - k5);
    }
    else if (k3 >= k5)
    {
        for (i = 0; i <= 3 * d; i++)
            for (j = 0; j <= 3 * d; j++)
                L[i][j] = L[i][j] * pow(2, k3 - k1);
        for (i = d + 1; i <= d + n; i++)
            L[i][i] = 1;
        for (i = d + n + 1; i <= d + 2 * n; i++)
            L[i][i] = pow(2, k3 - k5);
        for (i = d + 2 * n + 1; i <= 2 * d + n; i++)
            L[i][i] = 1;
        for (i = 2 * d + n + 1; i <= 3 * d; i++)
            L[i][i] = pow(2, k3 - k5);
    }
    else
    {
        for (i = 0; i <= 3 * d; i++)
            for (j = 0; j <= 3 * d; j++)
                L[i][j] = L[i][j] * pow(2, k5 - k1);
        for (i = d + 1; i <= d + n; i++)
            L[i][i] = pow(2, k5 - k3);
        for (i = d + n + 1; i <= d + 2 * n; i++)
            L[i][i] = 1;
        for (i = d + 2 * n + 1; i <= 2 * d + n; i++)
            L[i][i] = pow(2, k5 - k3);
        for (i = 2 * d + n + 1; i <= 3 * d; i++)
            L[i][i] = 1;
    }

    // Reduce the lattice L
    BKZ_FP(L, 0.99, 20);

    // Recover the initial state
    if (k1 >= k3 && k1 >= k5)
    {
        j = 0;
        while (L[j][0] != pow(2, k1 - 1) && L[j][0] != -pow(2, k1 - 1))
            j = j + 1;
        if (L[j][0] == -pow(2, k1 - 1))
        {
            for (i = 0; i < n; i++)
            {
                a1[i] = L[j][i + 1] + pow(2, k1 - 1);
                a3[i] = (L[j][i + d + 1] + pow(2, k1 - 1)) / pow(2, k1 - k3);
                a5[i] = (L[j][i + d + n + 1] + pow(2, k1 - 1)) / pow(2, k1 - k5);
                a[i] = (a5[i] * pow(2, k1 + k2 + k3 + k4) + a4[i] * pow(2, k1 + k2 + k3) + a3[i] * pow(2, k1 + k2) + a2[i] * pow(2, k1) + a1[i]) % m;
            }
        }
        if (L[j][0] == pow(2, k1 - 1))
        {
            for (i = 0; i < n; i++)
            {
                a1[i] = -L[j][i + 1] + pow(2, k1 - 1);
                a3[i] = (-L[j][i + d + 1] + pow(2, k1 - 1)) / pow(2, k1 - k3);
                a5[i] = (-L[j][i + d + n + 1] + pow(2, k1 - 1)) / pow(2, k1 - k5);
                a[i] = (a5[i] * pow(2, k1 + k2 + k3 + k4) + a4[i] * pow(2, k1 + k2 + k3) + a3[i] * pow(2, k1 + k2) + a2[i] * pow(2, k1) + a1[i]) % m;
            }
        }
    }
    else if (k3 >= k5)
    {
        j = 0;
        while (L[j][0] != pow(2, k3 - 1) && L[j][0] != -pow(2, k3 - 1))
            j = j + 1;
        if (L[j][0] == -pow(2, k3 - 1))
        {
            for (i = 0; i < n; i++)
            {
                a1[i] = (L[j][i + 1] + pow(2, k3 - 1)) / pow(2, k3 - k1);
                a3[i] = L[j][i + d + 1] + pow(2, k3 - 1);
                a5[i] = (L[j][i + d + n + 1] + pow(2, k3 - 1)) / pow(2, k3 - k5);
                a[i] = (a5[i] * pow(2, k1 + k2 + k3 + k4) + a4[i] * pow(2, k1 + k2 + k3) + a3[i] * pow(2, k1 + k2) + a2[i] * pow(2, k1) + a1[i]) % m;
            }
        }
        if (L[j][0] == pow(2, k3 - 1))
        {
            for (i = 0; i < n; i++)
            {
                a1[i] = (-L[j][i + 1] + pow(2, k3 - 1)) / pow(2, k3 - k1);
                a3[i] = -L[j][i + d + 1] + pow(2, k3 - 1);
                a5[i] = (-L[j][i + d + n + 1] + pow(2, k3 - 1)) / pow(2, k3 - k5);
                a[i] = (a5[i] * pow(2, k1 + k2 + k3 + k4) + a4[i] * pow(2, k1 + k2 + k3) + a3[i] * pow(2, k1 + k2) + a2[i] * pow(2, k1) + a1[i]) % m;
            }
        }
    }
    else
    {
        j = 0;
        while (L[j][0] != pow(2, k5 - 1) && L[j][0] != -pow(2, k5 - 1))
            j = j + 1;
        if (L[j][0] == -pow(2, k5 - 1))
        {
            for (i = 0; i < n; i++)
            {
                a1[i] = (L[j][i + 1] + pow(2, k5 - 1)) / pow(2, k5 - k1);
                a3[i] = (L[j][i + d + 1] + pow(2, k5 - 1)) / pow(2, k5 - k3);
                a5[i] = L[j][i + d + n + 1] + pow(2, k5 - 1);
                a[i] = (a5[i] * pow(2, k1 + k2 + k3 + k4) + a4[i] * pow(2, k1 + k2 + k3) + a3[i] * pow(2, k1 + k2) + a2[i] * pow(2, k1) + a1[i]) % m;
            }
        }
        if (L[j][0] == pow(2, k3 - 1))
        {
            for (i = 0; i < n; i++)
            {
                a1[i] = (-L[j][i + 1] + pow(2, k5 - 1)) / pow(2, k5 - k1);
                a3[i] = (-L[j][i + d + 1] + pow(2, k5 - 1)) / pow(2, k5 - k3);
                a5[i] = -L[j][i + d + n + 1] + pow(2, k5 - 1);
                a[i] = (a5[i] * pow(2, k1 + k2 + k3 + k4) + a4[i] * pow(2, k1 + k2 + k3) + a3[i] * pow(2, k1 + k2) + a2[i] * pow(2, k1) + a1[i]) % m;
            }
        }
    }

    end = clock();
    double time = (double)(end - start) / CLK_TCK;
    cout << "Step 4: time of recovering the initial state: " << time << "s\n";

    return a;
}

void Experiment_1_1(vec_ZZ a, ZZ m, int n, int k)
{
    vec_ZZ y, recovered_a;
    mat_ZZ C;
    ZZ_pX f;
    int r, t, N, d, alpha, beta;
    int i;

    r = 71;
    t = 71;
    N = r + t - 1;
    alpha = 17;
    beta = 14;
    d = 30;

    y.SetLength(N);
    for (i = 0; i < N; i++)
    {
        y[i] = a[i] >> beta;
    }

    C = Search_linear_relations_high_m(m, y, beta, r, t);
    f = Recover_the_coefficients(C, m, n, r);
    cout << "f = " << f << "\n\n\n";
    recovered_a = Recover_the_initial_state_high(y, f, m, n, d, beta);
    cout << "Number of truncated digits for recovering the initial state = " << d << "\n";
    cout << "a = " << recovered_a << "\n\n\n\n\n";
}

void Experiment_2_1(vec_ZZ a, int n, int k)
{
    vec_ZZ y, recovered_a;
    mat_ZZ C;
    ZZ m;
    ZZ_pX f;
    int r, t, N, d, alpha, beta;
    int i;

    r = 72;
    t = 70;
    N = r + t - 1;
    alpha = 17;
    beta = 14;
    d = 30;

    y.SetLength(N);
    for (i = 0; i < N; i++)
    {
        y[i] = a[i] >> beta;
    }

    C = Search_linear_relations_high(y, alpha, r, t);
    m = Recover_the_modulus(C, n, k, r);
    cout << "m = " << m << "\n\n\n";
    f = Recover_the_coefficients(C, m, n, r);
    cout << "f = " << f << "\n\n\n";
    recovered_a = Recover_the_initial_state_high(y, f, m, n, d, beta);
    cout << "Number of truncated digits for recovering the initial state = " << d << "\n";
    cout << "a = " << recovered_a << "\n\n\n\n\n";
}

void Experiment_3_1(vec_ZZ a, int n, int k)
{
    vec_ZZ y, recovered_a;
    mat_ZZ A, C;
    ZZ m;
    ZZ_pX f;
    int r, t, N, d, l;
    int i;

    // Experiment_3_1
    int K[3] = { 9, 17, 5 };
    int V[3] = { 0, 9, 26 };
    r = 70;
    t = 70;
    N = r + t - 1;
    d = 30;
    l = 3;

    // Experiment_3_2
    /*int K[3] = { 12, 14, 5 };
    int V[3] = { 0, 12, 26 };
    r = 105;
    t = 105;
    N = r + t - 1;
    d = 36;
    l = 3;*/

    A.SetDims(l, N);
    for (i = 0; i < N; i++)
    {
        A[0][i] = a[i] % power2_ZZ(K[0]);
        A[1][i] = (a[i] >> V[1]) % power2_ZZ(K[1]);
        A[2][i] = a[i] >> V[2];
    }

    C = Search_linear_relations_multi_segment_3(A[1], K[1], r, t);
    m = Recover_the_modulus(C, n, k, r);
    cout << "m = " << m << "\n\n\n";
    f = Recover_the_coefficients(C, m, n, r);
    cout << "f = " << f << "\n\n\n";
    recovered_a = Recover_the_initial_state_multi_segment_3(A[1], f, m, n, d, K[0], K[1], K[2]);
    cout << "Number of truncated digits for recovering the initial state = " << d << "\n";
    cout << "a = " << recovered_a << "\n\n\n\n\n";
}

void Experiment_3_3(vec_ZZ a, int n, int k)
{
    vec_ZZ y, recovered_a;
    mat_ZZ A, C;
    ZZ m;
    ZZ_pX f;
    int r, t, N, d, l;
    int i;

    // Experiment_3_3
    int K[5] = { 6, 3, 3, 14, 5 };
    int V[5] = { 0, 6, 9, 12, 26 };
    r = 104;
    t = 104;
    N = r + t - 1;
    d = 30;
    l = 5;

    // Experiment_3_4
    /*int K[5] = { 6, 2, 3, 15, 5 };
    int V[5] = { 0, 6, 8, 11, 26 };
    r = 90;
    t = 90;
    N = r + t - 1;
    d = 31;
    l = 5;*/

    A.SetDims(l, N);
    for (i = 0; i < N; i++)
    {
        A[0][i] = a[i] % power2_ZZ(K[0]);
        A[1][i] = (a[i] >> V[1]) % power2_ZZ(K[1]);
        A[2][i] = (a[i] >> V[2]) % power2_ZZ(K[2]);
        A[3][i] = (a[i] >> V[3]) % power2_ZZ(K[3]);
        A[4][i] = a[i] >> V[4];
    }

    C = Search_linear_relations_multi_segment_5(A[1], A[3], V[1], V[3], V[4], r, t);
    m = Recover_the_modulus(C, n, k, r);
    cout << "m = " << m << "\n\n\n";
    f = Recover_the_coefficients(C, m, n, r);
    cout << "f = " << f << "\n\n\n";
    recovered_a = Recover_the_initial_state_multi_segment_5(A[1], A[3], f, m, n, d, K[0], K[1], K[2], K[3], K[4]);
    cout << "Number of truncated digits for recovering the initial state = " << d << "\n";
    cout << "a = " << recovered_a << "\n\n\n\n\n";
}


void main()
{
    vec_ZZ a;
    ZZ m;
    int n, k, i;

    n = 16;
    k = 31;
    m = 2147483647;

    a.SetLength(300);
    a[0] = 5113033;
    a[1] = 101513454;
    a[2] = 342517303;
    a[3] = 602278014;
    a[4] = 545182365;
    a[5] = 2109832480;
    a[6] = 83899935;
    a[7] = 1353806595;
    a[8] = 56274808;
    a[9] = 220365776;
    a[10] = 1900047307;
    a[11] = 1287703311;
    a[12] = 754976298;
    a[13] = 344592060;
    a[14] = 2072123851;
    a[15] = 1906585150;
    for (i = n; i < 300; i++)
    {
        a[i] = (power2_ZZ(15) * a[i - 1] + power2_ZZ(17) * a[i - 3] + power2_ZZ(21) * a[i - 6] + power2_ZZ(20) * a[i - 12] + 257 * a[i - 16]) % m;
    }

    cout << "***Experiment 1.1***\n";
    Experiment_1_1(a, m, n, k);

    cout << "***Experiment 2.1***\n";
    Experiment_2_1(a, n, k);

    cout << "***Experiment 3.1***\n";
    Experiment_3_1(a, n, k);

    cout << "***Experiment 3.3***\n";
    Experiment_3_3(a, n, k);
}