Recover the truncated ZUC's driving sequence

// This code is used to recover the truncated ZUC's driving sequence.
// Suppose that the modulus, the coefficients and the initial states 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 and the initial states.
// 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
mat_ZZ Search_linear_relations(vec_ZZ y, int r, int t)
{
    clock_t start, end;
    start = clock();
    mat_ZZ L, C;
    int i, j;

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

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

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

    // Choose 3 polynomials that annihilate the sequence a
    for (i = 0; i < 3; 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";

    return C;
}

// Recover the modulus by the resultant
ZZ Recover_the_modulus(mat_ZZ C, int r)
{
    clock_t start, end;
    start = clock();
    ZZ m, res1, res2, res3, m_power;
    ZZX k1, k2, k3;
    int i;

    k1.SetLength(r);
    k2.SetLength(r);
    k3.SetLength(r);

    // Transform the vectors into polynomials
    for (i = 0; i < r; i++)
    {
        SetCoeff(k1, i, C[0][i]);
        SetCoeff(k2, i, C[1][i]);
        SetCoeff(k3, i, C[2][i]);
    }

    // Compute the resultants of the polynomials
    res1 = resultant(k1, k2);
    res2 = resultant(k1, k3);
    res3 = resultant(k2, k3);

    // Compute the greatest common factor of the resultants to obtain m to the power of n
    m_power = GCD(GCD(res1, res2), res3);
    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\n";

    return m;
}

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

    C_m.SetDims(2, r);

    ofstream fout("C.txt");
    for (i = 0; i < 2; 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 < 2; 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;
        }
        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]);
    }

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

    return GCD(f, g);
}

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

    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
    G_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\n";

    return a;
}

void main()
{
    vec_ZZ y, a;
    mat_ZZ C;
    ZZ m;
    ZZ_pX f;
    int n, r, t, N, alpha, beta;
    int i, j;

    n = 16;
    r = 175;
    t = 65;
    N = r + t - 1;
    alpha = 17;
    beta = 14;

    y.SetLength(N);
    a.SetLength(n);

    long Y[239] = { 312,6195,20905,36760,33275,128773,5120,82629,3434,13450,115969,78595,46080,21032,126472,116368,39147,118169,53608,16848,30654,130814,92427,66994,64940,95674,83420,12319,51686,121628,103101,4578,55561,114165,99105,66029,11183,33637,126822,45172,71211,75549,31249,100570,57872,40625,53029,98594,35382,114241,86331,100286,74451,71452,113487,22824,119007,42576,9420,18312,62834,3099,103626,124788,28318,119993,120470,1417,66987,121877,100551,36538,107257,95791,90747,11856,78683,27967,33439,112227,42492,76803,73986,127640,127988,113849,62073,121003,67257,62990,28618,14066,37196,88819,91188,119086,47509,14040,46365,117363,96302,61501,69537,19261,5302,29632,112053,60563,73212,123325,7485,52451,92665,108139,31895,6630,28136,97771,1874,10433,123318,98757,45897,34866,118234,75543,95739,41629,82990,25310,28330,21183,102442,39164,13744,125171,78275,42927,27280,113230,115192,108433,91178,101247,100622,63782,105110,40546,20787,97620,33627,97857,54439,74890,65682,117061,34213,28237,72077,40838,44105,41674,48232,10965,79577,106188,27672,79314,114402,64028,44159,130067,47639,20721,16188,1301,53216,29228,52467,45721,17559,28506,122718,84634,109576,34973,129964,10245,48392,100922,5358,128436,37415,74936,89446,101529,75414,39696,7480,20690,113561,65658,28894,48013,100734,80683,119101,50176,38679,50948,30958,67432,120501,62535,103920,69944,73629,7114,19565,10250,685,37711,17024,38274,53345,89835,93532,61993,35816,55362,11054,112389,34513,44560,123360,26732,69001,67919,67745 };
    for (i = 0; i < N; i++)
        y[i] = ZZ(Y[i]);

    C = Search_linear_relations(y, r, t);
    m = Recover_the_modulus(C, r);
    cout << "m = " << m << "\n\n\n";
    f = Recover_the_coefficients(C, m, r);
    cout << "f = " << f << "\n\n\n";
    a = Recover_the_initial_state(y, f, m, n, beta);
    cout << "a = " << a << "\n\n\n";
}